Add carra2vobs

data_files/carra2vobs.f90

+PROGRAM carra2vobs
+
+  INTEGER nmax,nsmax
+
+  PARAMETER (nmax = 3000, nsmax = 12000)
+
+  REAL    rval_lat, rval_lon, rval_hei, rval_snow, rval_prec24
+  REAL    rval_t2m, rval_ff, rval_pmsl, rval_td2m, rval_dd
+  REAL    rval_baro, rval_ps, rval_rh, rval_n, rval_vi, rval_prec1, rval_prec12
+  REAL    rval_year, rval_month, rval_day
+  REAL    rval_hour, rval_minute
+  REAL    rval_block, rval_station, rval_clim
+  REAL    lat(nmax), lon(nmax), hei(nmax), prec24(nmax)
+  REAL    t2m(nmax), ff(nmax), pmsl(nmax), td2m(nmax), dd(nmax)
+  REAL    ps(nmax), rh(nmax), n8(nmax), vis(nmax), prec1(nmax), prec12(nmax)
+  REAL    BLOCK(nmax), station(nmax), clim(nmax)
+  REAL    slat(nsmax),slon(nsmax),ilat(nsmax),ilon(nsmax)
+  REAL    qsat,pp,lat1,lon1,lat2,lon2,elev
+  INTEGER iy,im,id,ih,idate,syn_num(nsmax),ice1_num(nsmax),ice2_num(nsmax),nice
+  INTEGER anum
+  CHARACTER*10 cval_date,cdate
+  CHARACTER*8 cval_day
+  CHARACTER*40 name
+  CHARACTER*90 line
+  CHARACTER*7 indic, cstn(nmax)
+  CHARACTER*21 filename_ic/'IC/obs_yyyymmddhh.txt'/
+  CHARACTER*21 filename_se/'SE/obs_yyyymmddhh.txt'/
+  CHARACTER*21 filename_fi/'FI/obs_yyyymmddhh.txt'/
+  CHARACTER*21 filename_dk/'DK/obs_yyyymmddhh.txt'/
+  CHARACTER*36 filename_no1/'NO/carra_Norway_synop_yyyymmddhh.txt'/
+  CHARACTER*38 filename_no2/'NO/carra_Norway_SA_RR24_yyyymmddhh.txt'/
+
+  CHARACTER*3 country(6)/'DK ','FI ','IC ','NO1','NO2','SE '/
+
+  LOGICAL found
+
+  PARAMETER (rmissval=-99.)
+  PARAMETER (iversion = 4, ntemp = 0, nvar = 12)
+
+  CHARACTER*256 CINP,COUT,CARG(4),INPATH 
+                                                                      
+  CALL getenv('INPATH',INPATH)                              
+
+  NARG=IARGC()
+  pi   = 4.*ATAN(1.)
+  rad = pi/180.
+  rrad = 1./rad
+
+  DO J=1,NARG
+     CALL GETARG(J,CARG(J))
+  ENDDO
+
+  IF(NARG<1) THEN
+     PRINT*,'USAGE: exe date'
+     STOP
+  END IF
+
+  CINP = TRIM(INPATH)//'allsynop.list'
+  WRITE(*,*)TRIM(CINP)
+
+  OPEN(37,FILE=CINP,ACTION='read', STATUS='OLD',IOSTAT=IEOF)
+  
+  IF (IEOF /= 0) STOP 'OPEN FAILED ON INPUT FILE allsynop.list'
+
+  icount = 0
+
+  DO WHILE (IEOF == 0)
+
+     READ(37,*, IOSTAT=IEOF)num,lat1,lon1,elev,name
+    
+     IF  (IEOF .EQ. 0) THEN  
+
+        icount = icount + 1
+
+        syn_num(icount) = num
+        slat(icount) = lat1
+        slon(icount) = lon1
+
+     ENDIF
+
+  ENDDO
+
+  CLOSE(37)
+
+  nsyn = icount
+
+  CINP = TRIM(INPATH)//'iceland_remap.txt'
+  WRITE(*,*)TRIM(CINP)
+
+  OPEN(37,FILE=CINP,ACTION='read', STATUS='OLD',IOSTAT=IEOF)
+  
+  IF (IEOF /= 0) STOP 'OPEN FAILED ON INPUT FILE iceland_remap.txt'
+
+  icount = 0
+
+  DO WHILE (IEOF == 0)
+
+     READ(37,'(a)', IOSTAT=IEOF)line
+    
+     IF  (IEOF .EQ. 0) THEN  
+
+        icount = icount + 1
+
+        READ(line(1:7),'(i7)')ice1_num(icount)
+        READ(line(9:15),'(i7)')ice2_num(icount)
+        READ(line(17:25),'(f9.5)')ilon(icount)
+        READ(line(27:34),'(f8.5)')ilat(icount)
+
+     ENDIF
+
+  ENDDO
+
+  CLOSE(37)
+
+  nice = icount
+
+  cdate = CARG(1)
+  READ(cdate,'(i10)')idate
+  COUT = 'vobs.txt'
+
+  CALL split(idate,iy,im,id,ih)
+  
+  icount = 0
+
+  DO k = 1,6
+
+     IF(country(k) == 'NO2' .AND. ih .NE. 06)CYCLE
+
+     SELECT CASE (TRIM(country(k)))
+
+     CASE ('NO1')
+
+        WRITE(filename_no1(23:32),'(i4,3i2.2)')iy,im,id,ih
+        CINP=TRIM(INPATH)//filename_no1
+        WRITE(*,*)TRIM(CINP)
+ 
+     CASE ('NO2')
+
+        WRITE(filename_no2(25:34),'(i4,3i2.2)')iy,im,id,ih
+        CINP=TRIM(INPATH)//filename_no2
+        WRITE(*,*)TRIM(CINP)
+ 
+     CASE ('FI ')
+
+        WRITE(filename_fi(8:17),'(i4,3i2.2)')iy,im,id,ih
+        CINP=TRIM(INPATH)//filename_fi
+        WRITE(*,*)TRIM(CINP)
+ 
+     CASE ('IC ')
+
+        WRITE(filename_ic(8:17),'(i4,3i2.2)')iy,im,id,ih
+        CINP=TRIM(INPATH)//filename_ic
+        WRITE(*,*)TRIM(CINP)
+ 
+     CASE ('SE ')
+
+        WRITE(filename_se(8:17),'(i4,3i2.2)')iy,im,id,ih
+        CINP=TRIM(INPATH)//filename_se
+        WRITE(*,*)TRIM(CINP)
+ 
+     CASE ('DK ')
+
+        WRITE(filename_dk(8:17),'(i4,3i2.2)')iy,im,id,ih
+        CINP=TRIM(INPATH)//filename_dk
+        WRITE(*,*)TRIM(CINP)
+ 
+     CASE default        
+      
+        WRITE(*,*)'not allowed'
+
+     END SELECT
+
+     OPEN(37,FILE=CINP,ACTION='read', STATUS='OLD',IOSTAT=IEOF)
+  
+     IF (IEOF /= 0) STOP 'OPEN FAILED ON INPUT FILE'
+
+     DO WHILE (IEOF == 0)
+
+        SELECT CASE (country(k))
+
+        CASE ('FI ','DK ','SE ')
+      
+           READ(37,*, IOSTAT=IEOF) cval_date, rval_block, rval_station, rval_lat, rval_lon, &
+                rval_hei, rval_baro, rval_t2m, rval_rh, rval_td2m, rval_snow, rval_ps, &
+                rval_pmsl, rval_ff, rval_dd, rval_prec24, rval_prec12, rval_prec1, &
+                rval_n, rval_vi
+
+        CASE('NO2')
+
+           rval_baro = -9.
+           rval_block = -9.
+           rval_t2m = -9.
+           rval_td2m = -9.
+           rval_rh = -9.
+           rval_pmsl = -9.
+           rval_ps = -9.
+           rval_ff = -9.
+           rval_dd = -9.
+           rval_n = -9.
+           rval_vi = -9.
+           rval_prec1 = -9.
+           rval_prec12 = -9.
+
+           READ(37,*, IOSTAT=IEOF) cval_date, rval_station, rval_lat, rval_lon, &
+                rval_hei, rval_snow, rval_prec24
+
+           IF(rval_station > 0.)THEN
+
+              iiIII = NINT(rval_station)
+              ii = iiIII/1000
+              if(ii == 1)then
+                 rval_block = float(ii)
+              endif
+              iii = iiIII - ii*1000
+              rval_station = float(iii)
+
+           ENDIF
+
+        CASE('NO1')
+
+           rval_block = -9.
+           rval_prec24 = -9.
+           rval_snow = -9.
+
+           READ(37,*, IOSTAT=IEOF) cval_date, rval_station, rval_lat, rval_lon, &
+                rval_hei, rval_baro, rval_t2m, rval_rh, rval_td2m, rval_pmsl, rval_ps, &
+                rval_ff, rval_dd, rval_prec12, rval_prec1, rval_n, rval_vi
+
+           IF(rval_station > 0.)THEN
+
+              iiIII = NINT(rval_station)
+              ii = iiIII/1000
+              if(ii == 1)then
+                 rval_block = float(ii)
+              endif
+              iii = iiIII - ii*1000
+              rval_station = float(iii)
+
+           ENDIF
+
+        CASE('IC ')
+
+           READ(37,*, IOSTAT=IEOF) cval_date, rval_block, rval_station, rval_lon, rval_lat, &
+                rval_hei, rval_baro, rval_t2m, rval_rh, rval_td2m, rval_snow, rval_ps, &
+                rval_pmsl, rval_ff, rval_dd, rval_prec24, rval_prec12, rval_prec1, &
+                rval_n, rval_vi
+
+        CASE default
+           
+           WRITE(*,*)carg(2),' not allowed'
+           STOP
+
+        END SELECT
+
+        IF  (IEOF .EQ. 0) THEN  
+
+           found = .TRUE.
+
+           SELECT CASE (country(k))
+
+           CASE ('DK ')
+
+              IF (rval_station < 1000.) THEN
+              
+                 WRITE(indic(1:4),'(i4.4)')INT(rval_block)
+                 WRITE(indic(5:7),'(i3.3)')INT(rval_station)
+
+              ELSE
+           
+                 dist = 99999.
+
+                 DO i = 1,nsyn
+
+                    lat1 = rval_lat*rad
+                    lon1 = rval_lon*rad
+
+                    lat2 = slat(i)*rad
+                    lon2 = slon(i)*rad
+
+                    d = 2.*ASIN(SQRT((SIN((lat1-lat2)*0.5))**2 + &
+                         COS(lat1)*COS(lat2)*(SIN((lon1-lon2)*0.5))**2))
+
+                    d = 1852.*60*d*rrad
+                    
+                    IF( d < dist)THEN
+                       dist = d
+                       nval = i
+                    ENDIF
+
+                 ENDDO
+
+                 IF(dist < 2000.)THEN
+                    WRITE(indic,'(i7.7)')syn_num(nval)
+                 ELSE
+                    WRITE(*,*)'no number found: ',rval_lat,rval_lon
+                    found = .FALSE.
+                 ENDIF
+              
+              ENDIF
+            
+           CASE ('NO1','NO2')
+
+              IF (rval_block == 1.) THEN
+              
+                 WRITE(indic(1:4),'(i4.4)')INT(rval_block)
+                 WRITE(indic(5:7),'(i3.3)')INT(rval_station)
+
+              ELSE
+           
+                 dist = 99999.
+
+                 DO i = 1,nsyn
+
+                    lat1 = rval_lat*rad
+                    lon1 = rval_lon*rad
+
+                    lat2 = slat(i)*rad
+                    lon2 = slon(i)*rad
+
+                    d = 2.*ASIN(SQRT((SIN((lat1-lat2)*0.5))**2 + &
+                         COS(lat1)*COS(lat2)*(SIN((lon1-lon2)*0.5))**2))
+                    
+                    d = 1852.*60*d*rrad
+                    
+                    IF( d < dist)THEN
+                       dist = d
+                       nval = i
+                    ENDIF
+
+                 ENDDO
+
+                 IF(dist < 2000.)THEN
+                    WRITE(indic,'(i7.7)')syn_num(nval)
+                 ELSE
+                    found = .FALSE.
+                 ENDIF
+              
+              ENDIF
+            
+           CASE ('SE ')
+
+              IF (rval_station < 1000.) THEN
+              
+                 WRITE(indic(1:4),'(i4.4)')INT(rval_block)
+                 WRITE(indic(5:7),'(i3.3)')INT(rval_station)
+
+              ELSE
+
+                 indic(1:1) = '2'
+                 WRITE(indic(2:7),'(i6.6)')INT(rval_station)           
+              
+              ENDIF
+            
+           CASE ('FI ')
+
+              IF(rval_station < 1000.)THEN
+
+                 WRITE(indic(1:4),'(i4.4)')INT(rval_block)
+                 WRITE(indic(5:7),'(i3.3)')INT(rval_station)
+
+              ELSE
+
+                 indic(1:1) = '2'
+                 WRITE(indic(2:7),'(i6.6)')INT(rval_station)           
+
+              ENDIF
+
+           CASE ('IC ')
+
+              DO i = 1,nice
+
+                 IF(INT(rval_station) == ice1_num(i))THEN
+                    WRITE(indic,'(i7.7)')ice2_num(i)
+                 ENDIF
+
+              ENDDO
+              
+              IF(INT(rval_station) == 4180 .and. idate <= 2007123123)indic = '4005552'
+
+           END SELECT
+
+           READ(cval_date(1:4),*)  rval_year
+           READ(cval_date(5:6),*)  rval_month
+           READ(cval_date(7:8),*)  rval_day
+           READ(cval_date(9:10),*) rval_hour
+
+                        
+           IF (rval_hei .LT. 0.) THEN
+              rval_hei = rmissval
+           ENDIF
+            
+           IF (rval_baro .LT. 0.) THEN
+              rval_baro = rmissval
+           ENDIF
+            
+           IF (rval_t2m .LT. 0.) THEN
+              rval_t2m = rmissval
+           ENDIF
+
+           IF (rval_td2m .LT. 0.) THEN
+              rval_td2m = rmissval
+           ENDIF
+
+           IF (rval_rh .LT. 0.) THEN
+              rval_rh = rmissval
+           ENDIF
+
+           IF (rval_dd .LT. 0.) THEN
+              rval_dd = rmissval
+           ENDIF
+
+           IF (rval_ff .LT. 0.) THEN
+              rval_ff = rmissval
+           ENDIF
+
+           IF (rval_snow .LT. 0.) THEN
+              rval_snow = rmissval
+           ENDIF
+
+           IF (rval_pmsl .LT. 0.) THEN
+              rval_pmsl = rmissval
+           ENDIF
+
+           IF (rval_ps .LT. 0.) THEN
+              rval_ps = rmissval
+           ENDIF
+
+           IF (rval_prec24 .LT. 0.) THEN
+              rval_prec24 = rmissval
+           ENDIF
+
+           IF (rval_prec12 .LT. 0.) THEN
+              rval_prec12 = rmissval
+           ENDIF
+
+           IF (rval_prec1 .LT. 0.) THEN
+              rval_prec1 = rmissval
+           ENDIF
+
+           IF (rval_n .LT. 0.) THEN
+              rval_n = rmissval
+           ENDIF
+
+           IF (rval_vi .LT. 0.) THEN
+              rval_vi = rmissval
+           ENDIF
+
+           IF(rval_rh == rmissval)THEN
+
+              IF(rval_td2m .NE. rmissval .AND. rval_t2m .NE. rmissval)THEN
+                 IF(rval_ps .NE. rmissval)THEN
+                    pp = rval_ps
+                 ELSE
+                    pp = 101326.
+                 ENDIF
+                 rval_rh = 100.*qsat(pp,rval_td2m)/qsat(pp,rval_t2m)
+                 rval_rh = MIN(rval_rh,100.)
+              ENDIF
+           ENDIF
+
+           IF(found)THEN
+
+              icount = icount + 1
+
+              cstn(icount) = indic
+
+              lat(icount) = rval_lat
+              lon(icount) = rval_lon
+              hei(icount) = rval_hei
+              pmsl(icount) = rval_pmsl
+              ps(icount) = rval_ps
+              t2m(icount) = rval_t2m
+              td2m(icount) = rval_td2m
+              dd(icount) = rval_dd
+              ff(icount) = rval_ff
+              rh(icount) = rval_rh
+              prec24(icount) = rval_prec24
+              prec12(icount) = rval_prec12
+              prec1(icount) = rval_prec1
+              n8(icount) = rval_n 
+              vis(icount) = rval_vi
+
+           ENDIF
+            
+        ENDIF
+
+     ENDDO
+
+  ENDDO
+
+  WRITE(*,*)'Number of obs: ',icount
+
+  OPEN(38,FILE=COUT,ACTION='write', STATUS='unknown',IOSTAT=IEOF)
+
+  WRITE(38,'(1X,3i6)')icount,ntemp,iversion
+  WRITE(38,*)nvar
+  WRITE(38,*)' NN               0'
+  WRITE(38,*)' DD               0'
+  WRITE(38,*)' FF               0'
+  WRITE(38,*)' TT               0'
+  WRITE(38,*)' TD               0'
+  WRITE(38,*)' RH               0'
+  WRITE(38,*)' PSS              0'
+  WRITE(38,*)' PS               0'
+  WRITE(38,*)' VI               0'
+  WRITE(38,*)' PE24            24'
+  WRITE(38,*)' PE              12'
+  WRITE(38,*)' PE1              1'
+ 
+  DO i = 1,icount
+
+     WRITE(38,'(1x,a,2f9.3,f8.1,12en13.4e2)')cstn(i),lat(i),lon(i),hei(i), &
+          n8(i),dd(i),ff(i),t2m(i),td2m(i),rh(i),ps(i),pmsl(i),vis(i),prec24(i), &
+          prec12(i),prec1(i)
+
+  ENDDO
+
+  CLOSE(38)
+     
+  STOP
+
+END PROGRAM carra2vobs
+
+SUBROUTINE split(iymdh,iy,im,id,ih)
+
+  IMPLICIT NONE
+
+  INTEGER iymdh,iy,im,id,ih
+
+  iy = iymdh/1000000
+  im = (iymdh-iy*1000000)/10000
+  id = (iymdh-iy*1000000-im*10000)/100
+  ih = iymdh-iy*1000000-im*10000-id*100
+
+  RETURN
+
+END SUBROUTINE split
+
+FUNCTION QSAT(P,T)
+
+!   COMPUTES SATURATED SPECIFIC HUMIDITY (QS=EPSILON*EMAX(T)/P)
+!   EVAPORATION PRESSURE CALCULATED OVER WATER
+
+  PARAMETER (c1es=610.78,c2es=17.269,c2is=21.875,c3es=273.16,c4es=35.86,c4is=7.66)
+  zc2(t)=c2is + MIN(1.,MAX(0.,(t-c3es+15)/15.))*( c2es-c2is )
+  zc4(t)=c4is + MIN(1.,MAX(0.,(t-c3es+15)/15.))*( c4es-c4is )
+  zes(t)=zc2(t)*(t-c3es)/(t-zc4(t))
+  E=c1es*EXP(zes(t))
+  QSAT=0.622*E/P
+  RETURN
+END FUNCTION QSAT