C ================================================================================================== PROGRAM REGRID_to_v5d C ------------------------------------------------------------------------------------------------------- C This code is used to transform pressure-coordinate REGRID output to V5D format C COMPILATION: make -f REGRID_to_v5d.f.m C EXECUTION: REGRID_to_v5d pars.REGRID output.v5d C See v5df.h for parameters MAXVARS, MAXTIMES, MAXROWS, MAXCOLUMNS, MAXLEVELS, MISSING, IMISSING C See below ROMU: For personal choices on variables etc... C ------------------------------------------------------------------------------------------------------- Implicit none Include "/usr/local/vis5d/src/v5df.h" Integer v5dsetlowlev,iargc C LOCAL VARIABLES (ROMU: nvarV5D is number of fields listed in pars.REGRID) Integer nvarV5D Parameter(nvarV5D=18) Integer imonth,julian,days_month(12) Integer NY(100),FILT(100) Integer varV5D(100),filtV5D(100) Real sumV5D(100),facV5D(100) Integer ivarV5D,iREGRIDunit,nvarREGRID, $ ivarREGRID,idimREGRID(100),ivar,ilev, $ intSS,YYYYo,YYYYf,MMo,MMf,DDo,DDf, $ HHo,HHf,MIo,MIf,SSo,SSf,num60SS,inum60SS, $ n,i,maxnl,time_means, $ mean_Tm,mean_Um,mean_Vm,mean_Hm,mean_RHm, $ numtimes_means Real a,conv,xn,xlatc,phi1,phic,cnst,yyc,ymcentri,centri,x,y Real*8 YYYYMMDDHHMISSo,YYYYMMDDHHMISSf Character*9 NAMEV5D(100),NAMEREGRID(100) Character*24 date(100),date_means(100) Character*100 parsREGRIDname,REGRIDname,v5dname C VARIABLES DESCRIBING REGRID OUTPUT (flag, BIG/SUB HEADERS AND FIELDS) Integer flag Integer BHI(50,20) Real BHR(20,20) Character*80 BHIC(50,20), BHRC(20,20) Integer ndim, start_index(4), end_index(4) Real xtime Character*4 staggering, ordering Character*24 current_date Character*9 name Character*25 unit Character*46 description Real f3d(MAXROWS,MAXCOLUMNS,MAXLEVELS),f2d(MAXROWS,MAXCOLUMNS), $ f1d(MAXLEVELS) C THE FOLLOWING VARIABLES DESCRIBE THE DATASET TO BE CONVERTED. YOU MUST INITIALIZE ALL C THESE VARIABLES WITH VALUES FROM YOUR FILE OR ASSIGN SUITABLE CONSTANTS. C SEE THE README FILE FOR DESCRIPTIONS OF THESE VARIABLES. Real G(MAXROWS,MAXCOLUMNS,MAXLEVELS), $ Tm(MAXROWS,MAXCOLUMNS,MAXLEVELS), $ Um(MAXROWS,MAXCOLUMNS,MAXLEVELS), $ Vm(MAXROWS,MAXCOLUMNS,MAXLEVELS), $ Hm(MAXROWS,MAXCOLUMNS,MAXLEVELS), $ RHm(MAXROWS,MAXCOLUMNS,MAXLEVELS) Integer ir,ic,ilev,nr,nc,nl(MAXVARS) Integer numtimes Integer numvars Character*10 varname(MAXVARS) Character*20 units(MAXVARS) Integer lowlev(MAXVARS) Integer dates(MAXTIMES) Integer times(MAXTIMES) Integer compressmode Integer projection Real proj_args(100) Integer vertical Real vert_args(MAXLEVELS) C INITIALIZE THE VARIABLES TO MISSING VALUES Data nr,nc / IMISSING, IMISSING / Data (nl(i),i=1,MAXVARS) / MAXVARS*IMISSING / Data (lowlev(i),i=1,MAXVARS) / MAXVARS*IMISSING / Data numtimes,numvars / IMISSING, IMISSING / Data (varname(i),i=1,MAXVARS) / MAXVARS*" " / Data (dates(i),i=1,MAXTIMES) / MAXTIMES*IMISSING / Data (times(i),i=1,MAXTIMES) / MAXTIMES*IMISSING / Data compressmode / 1 / Data projection / IMISSING / Data (proj_args(i),i=1,100) / 100*MISSING / Data vertical / IMISSING / Data (vert_args(i),i=1,MAXLEVELS) / MAXLEVELS*MISSING / C DAYS WIHIN THE 12 MONTHS (FOR NORMAL YEARS) Data days_month / 31,28,31,30,31,30,31,31,30,31,30,31 / C GET COMMAND LINE ARGUMENTS If (iargc() .ne. 2) then Print *,'Error: two filename arguments are needed.' Call Exit(1) Else Call Getarg(1, parsREGRIDname) Call Getarg(2, v5dname) Endif Print * Print * Print 10,'Input file (pars REGRID): ', parsREGRIDname Print 10,'Output file (V5D): ', v5dname 10 Format (A,A) C READ REGRID PARAMETERS FROM parsREGRIDname (SEE FOR USER CHOICES) Open(11,file=parsREGRIDname,status='old') Read(11,*) Read(11,*) Read(11,*) Read(11,*) REGRIDname Read(11,*) Read(11,*) Read(11,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') $ YYYYo,MMo,DDo,HHo,MIo,SSo julian=0 Do imonth=1,MMo-1 julian=julian+days_month(imonth) Enddo If ((MMo.Gt.2).And.(Mod(YYYYo-1964,4).Eq.0)) $ julian=julian+1 julian=julian+DDo YYYYMMDDHHMISSo= $ YYYYo*10000000000.d0+MMo*100000000.d0+ $ DDo*1000000.d0+HHo*10000.d0+MIo*100.d0+SSo Read(11,*) Read(11,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') $ YYYYf,MMf,DDf,HHf,MIf,SSf YYYYMMDDHHMISSf= $ YYYYf*10000000000.d0+MMf*100000000.d0+ $ DDf*1000000.d0+HHf*10000.d0+MIf*100.d0+SSf Read(11,*) Read(11,*) intSS num60SS=intSS/60 numtimes=0 Do While (YYYYMMDDHHMISSo.Le.YYYYMMDDHHMISSf) numtimes=numtimes+1 date(numtimes)=Char(YYYYo/1000 +48)// $ Char(YYYYo/100-(YYYYo/1000)*10 +48)// $ Char(YYYYo/10-(YYYYo/100)*10 +48 )// $ Char(YYYYo-(YYYYo/10)*10 +48)// $ '-'// $ Char(MMo/10 +48)// $ Char(MMo-(MMo/10)*10 +48)// $ '-'// $ Char(DDo/10 +48)// $ Char(DDo-(DDo/10)*10 +48)// $ '_'// $ Char(HHo/10 +48)// $ Char(HHo-(HHo/10)*10 +48)// $ ':'// $ Char(MIo/10 +48)// $ Char(MIo-(MIo/10)*10 +48)// $ ':'// $ Char(SSo/10 +48)// $ Char(SSo-(SSo/10)*10 +48)// $ '.0000' dates(numtimes)=(YYYYo-(YYYYo/100)*100)*1000+julian times(numtimes)=HHo*10000+MIo*100+SSo Do inum60SS=1,num60SS+1 If (inum60SS.Le.num60SS) Then SSo=SSo+60 Else SSo=SSo+intSS-num60SS*60 Endif If (SSo.Ge.60) Then SSo=SSo-60 MIo=Mio+1 If (MIo.Ge.60) Then MIo=MIo-60 HHo=HHo+1 If (HHo.Ge.24) Then HHo=HHo-24 DDo=DDo+1 julian=julian+1 If ( ((DDo.Gt.28).And.(MMo.Eq.2).And. $ (Mod(YYYYo-1964,4).Ne.0)) $ .Or. $ ((DDo.Gt.29).And.(MMo.Eq.2).And. $ (Mod(YYYYo-1964,4).Eq.0)) $ .Or. $ ((DDo.Gt.30).And. $ (MMo.Eq.4.Or.MMo.Eq.6.Or. $ MMo.Eq.9.Or.MMo.Eq.11)) $ .Or. $ (DDo.Gt.31) ) Then DDo=1 MMo=MMo+1 Endif If (MMo.Gt.12) Then MMo=1 YYYYo=YYYYo+1 Endif Endif Endif Endif Enddo YYYYMMDDHHMISSo= $ YYYYo*10000000000.d0+MMo*100000000.d0+ $ DDo*1000000.d0+HHo*10000.d0+MIo*100.d0+SSo Enddo Read(11,*) Read(11,*) Read(11,*) compressmode Read(11,*) Read(11,*) Read(11,*) time_means Read(11,*) Read(11,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') $ YYYYo,MMo,DDo,HHo,MIo,SSo YYYYMMDDHHMISSo= $ YYYYo*10000000000.d0+MMo*100000000.d0+ $ DDo*1000000.d0+HHo*10000.d0+MIo*100.d0+SSo Read(11,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') $ YYYYf,MMf,DDf,HHf,MIf,SSf YYYYMMDDHHMISSf= $ YYYYf*10000000000.d0+MMf*100000000.d0+ $ DDf*1000000.d0+HHf*10000.d0+MIf*100.d0+SSf Read(11,*) intSS If (time_means.Eq.1) Then num60SS=intSS/60 numtimes_means=0 Do While (YYYYMMDDHHMISSo.Le.YYYYMMDDHHMISSf) numtimes_means=numtimes_means+1 date_means(numtimes_means)= $ Char(YYYYo/1000 +48)// $ Char(YYYYo/100-(YYYYo/1000)*10 +48)// $ Char(YYYYo/10-(YYYYo/100)*10 +48 )// $ Char(YYYYo-(YYYYo/10)*10 +48)// $ '-'// $ Char(MMo/10 +48)// $ Char(MMo-(MMo/10)*10 +48)// $ '-'// $ Char(DDo/10 +48)// $ Char(DDo-(DDo/10)*10 +48)// $ '_'// $ Char(HHo/10 +48)// $ Char(HHo-(HHo/10)*10 +48)// $ ':'// $ Char(MIo/10 +48)// $ Char(MIo-(MIo/10)*10 +48)// $ ':'// $ Char(SSo/10 +48)// $ Char(SSo-(SSo/10)*10 +48)// $ '.0000' Do inum60SS=1,num60SS+1 If (inum60SS.Le.num60SS) Then SSo=SSo+60 Else SSo=SSo+intSS-num60SS*60 Endif If (SSo.Ge.60) Then SSo=SSo-60 MIo=Mio+1 If (MIo.Ge.60) Then MIo=MIo-60 HHo=HHo+1 If (HHo.Ge.24) Then HHo=HHo-24 DDo=DDo+1 If ( ((DDo.Gt.28).And.(MMo.Eq.2).And. $ (Mod(YYYYo-1964,4).Ne.0)) $ .Or. $ ((DDo.Gt.29).And.(MMo.Eq.2).And. $ (Mod(YYYYo-1964,4).Eq.0)) $ .Or. $ ((DDo.Gt.30).And. $ (MMo.Eq.4.Or.MMo.Eq.6.Or. $ MMo.Eq.9.Or.MMo.Eq.11)) $ .Or. $ (DDo.Gt.31) ) Then DDo=1 MMo=MMo+1 Endif If (MMo.Gt.12) Then MMo=1 YYYYo=YYYYo+1 Endif Endif Endif Endif Enddo YYYYMMDDHHMISSo= $ YYYYo*10000000000.d0+MMo*100000000.d0+ $ DDo*1000000.d0+HHo*10000.d0+MIo*100.d0+SSo Enddo Endif Read(11,*) Read(11,*) Read(11,*) Read(11,*) Do ivarV5D=1,nvarV5D Read(11,'(a9,1x,i10,1x,i9)') NAMEV5D(ivarV5D),NY(ivarV5D), $ FILT(ivarV5D) Enddo Close(11) C OPEN YOUR REGRID DATAFILE (REGRIDname) HERE AND READ ITS HEADER INFORMATION TO C INITIALIZE SOME OF THE ABOVE V5D PARAMETERS iREGRIDunit=20 Open(iREGRIDunit,file=REGRIDname,form='unformatted') Rewind(iREGRIDunit) Read(iREGRIDunit,err=9001,end=9002) flag Read(iREGRIDunit,err=9001,end=9002) BHI,BHR,BHIC,BHRC If (BHI(1,1).ne.2) Go To 9003 If (BHI(15,1).Ne.0) Then nr=BHI(16,1)-1 nc=BHI(17,1)-1 Else nr=BHI(5,1)-1 nc=BHI(6,1)-1 Endif maxnl=BHI(12,2)-1 nvarREGRID=0 Read(iREGRIDunit) flag Do While (flag.Eq.1) Read(iREGRIDunit) $ ndim,start_index,end_index,xtime,staggering, $ ordering,current_date,name,unit,description nvarREGRID=nvarREGRID+1 idimREGRID(nvarREGRID)=ndim NAMEREGRID(nvarREGRID)=name If (ndim.Eq.3) Then Read(iREGRIDunit) $ (((f3d(ir,ic,ilev),ir=1,nr+1),ic=1,nc+1),ilev=1,maxnl+1) Else If (ndim.Eq.2) Then Read(iREGRIDunit) ((f2d(ir,ic),ir=1,nr+1),ic=1,nc+1) Else If (ndim.Eq.1) Then Read(iREGRIDunit) (f1d(ilev),ilev=1,maxnl+1) If (ordering.Eq.'P ') Then vertical=3 Do ilev=1,maxnl vert_args(ilev)=0.01*f1d(ilev+1) Enddo Endif Endif Read(iREGRIDunit) flag Enddo Close(iREGRIDunit) numvars=0 Do ivarREGRID=1,nvarREGRID varV5D(ivarREGRID)=0 Enddo mean_Tm=0 mean_Um=0 mean_Vm=0 mean_Hm=0 mean_RHm=0 Print *, 'ID FIELD SMOOTHING UNITS' Print *, '---------------------------------------' Do ivarV5D=1,nvarV5D If (NY(ivarV5D).Eq.0) Go to 9004 Do ivarREGRID=1,nvarREGRID If (NAMEV5D(ivarV5D).Eq.NAMEREGRID(ivarREGRID)) Then numvars=numvars+1 varV5D(ivarREGRID)=numvars filtV5D(ivarREGRID)=FILT(ivarV5D) If (idimREGRID(ivarREGRID).Eq.3) Then nl(numvars)=maxnl Else nl(numvars)=1 Endif varname(numvars)=NAMEREGRID(ivarREGRID) lowlev(numvars)=0 C ROMU: Estan son las unidades que se tomaran (precisandose a veces para ello C transformaciones, a base de sumas o productos) If (NAMEREGRID(ivarREGRID).Eq.'U '.Or. $ NAMEREGRID(ivarREGRID).Eq.'V ') Then units(numvars)='m/s' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'RH ') Then units(numvars)='%' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'T '.Or. $ NAMEREGRID(ivarREGRID).EQ.'TSEASFC ') Then units(numvars)='C' sumV5D(ivarREGRID)=-273.16 facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'H '.Or. $ NAMEREGRID(ivarREGRID).Eq.'TERRAIN ') Then units(numvars)='m' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'PSEALVLD '.Or. $ NAMEREGRID(ivarREGRID).EQ.'PSEALVLC ') Then units(numvars)='mb' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=0.01 Else If (NAMEREGRID(ivarREGRID).Eq.'MAPFACCR '.Or. $ NAMEREGRID(ivarREGRID).Eq.'MAPFACDT '.Or. $ NAMEREGRID(ivarREGRID).Eq.'SNOWCOVR ') Then units(numvars)='no dim.' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'CORIOLIS ') Then units(numvars)='1/s' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'LATITCRS '.Or. $ NAMEREGRID(ivarREGRID).Eq.'LONGICRS '.Or. $ NAMEREGRID(ivarREGRID).Eq.'LATITDOT '.Or. $ NAMEREGRID(ivarREGRID).EQ.'LONGIDOT ') Then units(numvars)='degrees' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Else If (NAMEREGRID(ivarREGRID).Eq.'LAND USE ') Then units(numvars)='categ.' sumV5D(ivarREGRID)=0. facV5D(ivarREGRID)=1. Endif Print 578, numvars,NAMEREGRID(ivarREGRID), $ filtV5D(ivarREGRID),units(numvars) If (time_means.Eq.1) Then If (NAMEREGRID(ivarREGRID).Eq. 'T ') Then numvars=numvars+1 mean_Tm=numvars varname(numvars)= 'Tm ' nl(numvars)=maxnl lowlev(numvars)=0 units(numvars)=units(numvars-1) Print 578, numvars,'Tm ', $ filtV5D(ivarREGRID),units(numvars) Else If (NAMEREGRID(ivarREGRID).Eq.'U ') Then numvars=numvars+1 mean_Um=numvars varname(numvars)= 'Um ' nl(numvars)=maxnl lowlev(numvars)=0 units(numvars)=units(numvars-1) Print 578, numvars,'Um ', $ filtV5D(ivarREGRID),units(numvars) Else If (NAMEREGRID(ivarREGRID).Eq.'V ') Then numvars=numvars+1 mean_Vm=numvars varname(numvars)= 'Vm ' nl(numvars)=maxnl lowlev(numvars)=0 units(numvars)=units(numvars-1) Print 578, numvars,'Vm ', $ filtV5D(ivarREGRID),units(numvars) Else If (NAMEREGRID(ivarREGRID).Eq.'H ') Then numvars=numvars+1 mean_Hm=numvars varname(numvars)= 'Hm ' nl(numvars)=maxnl lowlev(numvars)=0 units(numvars)=units(numvars-1) Print 578, numvars,'Hm ', $ filtV5D(ivarREGRID),units(numvars) Else If (NAMEREGRID(ivarREGRID).Eq.'RH ') Then numvars=numvars+1 mean_RHm=numvars varname(numvars)= 'RHm ' nl(numvars)=maxnl lowlev(numvars)=0 units(numvars)=units(numvars-1) Print 578, numvars,'RHm ', $ filtV5D(ivarREGRID),units(numvars) Endif Endif 578 Format (1x,i2,5x,a9,5x,i5,9x,a20) Endif Enddo 9004 Continue Enddo C ROMU: Se asume que la proyeccion es siempre LAMBERT CONFORMAL C (de momento aun no estan incluidas las posibilidades de C Polar Stereo. y Mercator !!!!!!!!) If (BHI(7,1).Eq.1) Then projection=2 proj_args(1)=BHR(5,1) proj_args(2)=BHR(6,1) a=6370. conv=57.29578 xn=BHR(4,1) xlatc=BHR(2,1) phi1=BHR(6,1)/conv phic=(90.-xlatc)/conv cnst=a*sin(phi1)/(xn*tan(phi1/2.)**xn) yyc=-1.*cnst*(tan(phic/2.)**xn) ymcentri=-yyc/(0.001*BHR(1,1)) centri=(float(BHI(5,1))+1.)/2. y=ymcentri+centri If (BHI(15,1).Ne.0) Then y=y-BHR(10,1)+1. y=(y-1.)*float(BHI(20,1))+1. proj_args(3)=float(nr)-(y-0.5)+1. x=(float(BHI(6,1))+1.)/2. x=x-BHR(11,1)+1. x=(x-1.)*float(BHI(20,1))+1. proj_args(4)=x-0.5 proj_args(5)=-BHR(3,1) proj_args(6)=0.001*BHR(9,1) Else proj_args(3)=float(nr)-(y-0.5)+1. proj_args(4)=(float(BHI(6,1))+1.)/2.-0.5 proj_args(5)=-BHR(3,1) proj_args(6)=0.001*BHR(1,1) Endif C ROMU: Por alguna extranya razon se les ha de restar 1. al row y column del polo C (es como si el vis5d al calcular la proyeccion pensase que el punto superior C (izquierdo no es (1,1) sino (0,0) !!!!!!!!!) proj_args(3)=proj_args(3)-1. proj_args(4)=proj_args(4)-1. Endif C CREATE THE V5D FILE n = v5dcreate (v5dname, numtimes, numvars, nr, nc, nl, $ varname, times, dates, compressmode, $ projection, proj_args, vertical, vert_args) If (n.Eq.0) Then Print*,'PROBLEM: CAN NOT CREATE V5D FILE' Call Exit(1) Endif Do ivar=1,numvars Call v5dsetunits(ivar,units(ivar)) Enddo n=v5dsetlowlev(lowlev) C PRINT SOME V5D PARAMETERS Print * Print *, 'Number of rows = ',nr Print *, 'Number of columns = ',nc Print *, 'Number of levels for 3D fields = ',maxnl Print *, 'Vertical coordinate = ',vertical,' (pressure)' Print *, 'Map projection = ', $ projection,' (Lambert conformal)' Print *, 'Compressmode = ',compressmode Print * C CALL SUBROUTINE CONVERT TO WRITE REGRID OUTPUT AS V5D Call CONVERT (nr,nc,maxnl,G,f3d,f2d,f1d, $ REGRIDname,nvarREGRID, $ varV5D,filtV5D,sumV5D,facV5D,time_means, $ mean_Tm,mean_Um,mean_Vm,mean_Hm,mean_RHm, $ Tm,Um,Vm,Hm,RHm, $ numtimes,date, $ numtimes_means,date_means) C ERROR IN READING HEADER ......... 9001 Continue Print *,' ' Print *,'ERROR IN READING HEADER' Stop 9001 9002 Continue Print *,' ' Print *,'EOF DETECTED INITIALLY IN READING HEADER' Print *,'THIS DATA SET IS EMPTY, OR INCORRECTLY ASSIGNED' Stop 9002 C ERROR: THIS IS NOT A REGRID OUTPUT ......... 9003 Continue Print *,' ' Print *,'ERROR: NOT A REGRID OUTPUT' Stop 9003 END C ================================================================================================== C ================================================================================================== SUBROUTINE CONVERT (nr,nc,maxnl,G,f3d,f2d,f1d, $ REGRIDname,nvarREGRID, $ varV5D,filtV5D,sumV5D,facV5D,time_means, $ mean_Tm,mean_Um,mean_Vm,mean_Hm,mean_RHm, $ Tm,Um,Vm,Hm,RHm, $ numtimes,date, $ numtimes_means,date_means) Implicit none Include "/usr/local/vis5d/src/v5df.h" C LOCAL VARIABLES Integer nr,nc,maxnl,n,ir,ic,ilev,ifilt, $ ivarREGRID,nvarREGRID, $ iREGRIDunit,time_means, $ mean_Tm,mean_Um,mean_Vm,mean_Hm,mean_RHm, $ itime,numtimes, $ itime_means,numtimes_means, $ iprocess_time,iprocess_time_means Integer varV5D(100),filtV5D(100) Character*9 NAMEREGRID Real sumV5D(100),facV5D(100) Character*24 date(100),date_means(100) Character*100 REGRIDname C VARIABLES DESCRIBING REGRID OUTPUT (flag, BIG/SUB HEADERS AND FIELDS) Integer flag Integer BHI(50,20) Real BHR(20,20) Character*80 BHIC(50,20), BHRC(20,20) Integer ndim, start_index(4), end_index(4) Real xtime Character*4 staggering, ordering Character*24 current_date Character*9 name Character*25 unit Character*46 description Real G(nr,nc,maxnl), $ Tm(nr,nc,maxnl),Um(nr,nc,maxnl),Vm(nr,nc,maxnl), $ Hm(nr,nc,maxnl),RHm(nr,nc,maxnl), $ f3d(nr+1,nc+1,maxnl+1),f2d(nr+1,nc+1),f1d(maxnl+1) C OPEN REGRID OUTPUT FILE iREGRIDunit=20 Open(unit=iREGRIDunit,file=REGRIDname,form='unformatted') Rewind(iREGRIDunit) Read(iREGRIDunit) flag Read(iREGRIDunit) BHI,BHR,BHIC,BHRC C RESET TIME-MEAN FIELDS If (time_means.Eq.1) Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc Tm(ir,ic,ilev)=0. Um(ir,ic,ilev)=0. Vm(ir,ic,ilev)=0. Hm(ir,ic,ilev)=0. RHm(ir,ic,ilev)=0. Enddo Enddo Enddo Endif C RESET COUNTER OF TIMES REQUESTED FOR FIELDS AND MEAN FIELDS itime=0 itime_means=0 C LOOP OVER time series (READING FIELDS FROM REGRIDname, and PROCESSING C AND WRITTING IN v5dname if time period was requested) DO WHILE (itime.Lt.numtimes.Or.itime_means.Lt.numtimes_means) Read(iREGRIDunit,end=10000) flag ivarREGRID=0 DO WHILE (flag.Eq.1) Read(iREGRIDunit) $ ndim,start_index,end_index,xtime,staggering, $ ordering,current_date,name,unit,description ivarREGRID=ivarREGRID+1 NAMEREGRID=name If (ivarREGRID.Eq.1) Then Print *, 'Reading date = ',current_date iprocess_time=0 iprocess_time_means=0 If (current_date.Eq.date(itime+1)) Then itime=itime+1 iprocess_time=1 Print *, ' > Processing FIELDS ' Endif If (current_date.Eq.date_means(itime_means+1)) Then itime_means=itime_means+1 iprocess_time_means=1 Print *, ' > Processing MEAN FIELDS ' Endif Endif C 3D FIELDS If (ndim.Eq.3) Then C Dots (3D) If (staggering.Eq.'D ') Then Read(iREGRIDunit) f3d If (iprocess_time.Eq.0.And.iprocess_time_means.Eq.0) $ Goto 11 If (varV5D(ivarREGRID).Gt.0) Then If (sumV5D(ivarREGRID).Ne.0.) Then Do ilev=2,maxnl+1 Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,ilev-1)=.25* $ (f3d(ir,ic,ilev)+ $ f3d(ir+1,ic,ilev)+ $ f3d(ir,ic+1,ilev)+ $ f3d(ir+1,ic+1,ilev)) $ +sumV5D(ivarREGRID) Enddo Enddo Enddo Else If (facV5D(ivarREGRID).Ne.1.) Then Do ilev=2,maxnl+1 Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,ilev-1)=.25* $ (f3d(ir,ic,ilev)+ $ f3d(ir+1,ic,ilev)+ $ f3d(ir,ic+1,ilev)+ $ f3d(ir+1,ic+1,ilev)) $ *facV5D(ivarREGRID) Enddo Enddo Enddo Else Do ilev=2,maxnl+1 Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,ilev-1)=.25* $ (f3d(ir,ic,ilev)+ $ f3d(ir+1,ic,ilev)+ $ f3d(ir,ic+1,ilev)+ $ f3d(ir+1,ic+1,ilev)) Enddo Enddo Enddo Endif Do ifilt=1,filtV5D(ivarREGRID) Call Smooth(G,nr,nc,maxnl,maxnl,0.0625,0.0625) Enddo If (iprocess_time.Eq.0) Goto 111 n = V5dwrite( itime, varV5D(ivarREGRID), G ) If (n.Eq.0) Then Call Exit(1) Endif 111 Continue If (iprocess_time_means.Eq.0) Goto 1111 If (NAMEREGRID.Eq.'U ' $ .And.mean_Um.Gt.0) Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc Um(ir,ic,ilev)=Um(ir,ic,ilev)+G(ir,ic,ilev)/ $ Real(numtimes_means) Enddo Enddo Enddo Else If (NAMEREGRID.Eq.'V ' $ .And.mean_Vm.Gt.0) Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc Vm(ir,ic,ilev)=Vm(ir,ic,ilev)+G(ir,ic,ilev)/ $ Real(numtimes_means) Enddo Enddo Enddo Endif 1111 Continue Endif 11 Continue C Crosses (3D) Else If (staggering.Eq.'C ') Then Read(iREGRIDunit) f3d If (iprocess_time.Eq.0.And.iprocess_time_means.Eq.0) $ Goto 22 If (varV5D(ivarREGRID).Gt.0) Then If (sumV5D(ivarREGRID).Ne.0.) Then Do ilev=2,maxnl+1 Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,ilev-1)= $ f3d(ir,ic,ilev) $ +sumV5D(ivarREGRID) Enddo Enddo Enddo Else If (facV5D(ivarREGRID).Ne.1.) Then Do ilev=2,maxnl+1 Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,ilev-1)= $ f3d(ir,ic,ilev) $ *facV5D(ivarREGRID) Enddo Enddo Enddo Else Do ilev=2,maxnl+1 Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,ilev-1)= $ f3d(ir,ic,ilev) Enddo Enddo Enddo Endif Do ifilt=1,filtV5D(ivarREGRID) Call Smooth(G,nr,nc,maxnl,maxnl,0.0625,0.0625) Enddo If (NAMEREGRID.Eq.'RH ') Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc G(ir,ic,ilev)=Max(G(ir,ic,ilev),0.) Enddo Enddo Enddo Endif If (NAMEREGRID.Eq.'RH ') Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc G(ir,ic,ilev)=Min(G(ir,ic,ilev),100.) Enddo Enddo Enddo Endif If (iprocess_time.Eq.0) Goto 222 n = V5dwrite( itime, varV5D(ivarREGRID), G ) If (n.Eq.0) Then Call Exit(1) Endif 222 Continue If (iprocess_time_means.Eq.0) Goto 2222 If (NAMEREGRID.Eq.'T ' $ .And.mean_Tm.Gt.0) Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc Tm(ir,ic,ilev)=Tm(ir,ic,ilev)+G(ir,ic,ilev)/ $ Real(numtimes_means) Enddo Enddo Enddo Else If (NAMEREGRID.Eq.'H ' $ .And.mean_Hm.Gt.0) Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc Hm(ir,ic,ilev)=Hm(ir,ic,ilev)+G(ir,ic,ilev)/ $ Real(numtimes_means) Enddo Enddo Enddo Else If (NAMEREGRID.Eq.'RH ' $ .And.mean_RHm.Gt.0) Then Do ilev=1,maxnl Do ir=1,nr Do ic=1,nc RHm(ir,ic,ilev)=RHm(ir,ic,ilev)+ $ G(ir,ic,ilev)/ $ Real(numtimes_means) Enddo Enddo Enddo Endif 2222 Continue Endif 22 Continue Endif C 2D FIELDS Else If (ndim.Eq.2) Then C Dots (2D) If (staggering.Eq.'D ') Then Read(iREGRIDunit) f2d If (iprocess_time.Eq.0) Goto 33 If (varV5D(ivarREGRID).Gt.0) Then If (sumV5D(ivarREGRID).Ne.0.) Then Do ir=1,nr Do ic=1,nc G(nr-1+ir,ic,1)=.25* $ (f2d(ir,ic)+ $ f2d(ir+1,ic)+ $ f2d(ir,ic+1)+ $ f2d(ir+1,ic+1)) $ +sumV5D(ivarREGRID) Enddo Enddo Else If (facV5D(ivarREGRID).Ne.1.) Then Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,1)=.25* $ (f2d(ir,ic)+ $ f2d(ir+1,ic)+ $ f2d(ir,ic+1)+ $ f2d(ir+1,ic+1)) $ *facV5D(ivarREGRID) Enddo Enddo Else Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,1)=.25* $ (f2d(ir,ic)+ $ f2d(ir+1,ic)+ $ f2d(ir,ic+1)+ $ f2d(ir+1,ic+1)) Enddo Enddo Endif Do ifilt=1,filtV5D(ivarREGRID) Call Smooth(G,nr,nc,maxnl,1,0.0625,0.0625) Enddo n = V5dwrite( itime, varV5D(ivarREGRID), G ) If (n.Eq.0) Then Call Exit(1) Endif Endif 33 Continue C Crosses (2D) Else If (staggering.Eq.'C ') Then Read(iREGRIDunit) f2d If (iprocess_time.Eq.0) Goto 44 If (varV5D(ivarREGRID).Gt.0) Then If (sumV5D(ivarREGRID).Ne.0.) Then Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,1)= $ f2d(ir,ic) $ +sumV5D(ivarREGRID) Enddo Enddo Else If (facV5D(ivarREGRID).Ne.1.) Then Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,1)= $ f2d(ir,ic) $ *facV5D(ivarREGRID) Enddo Enddo Else Do ir=1,nr Do ic=1,nc G(nr+1-ir,ic,1)= $ f2d(ir,ic) Enddo Enddo Endif Do ifilt=1,filtV5D(ivarREGRID) Call Smooth(G,nr,nc,maxnl,1,0.0625,0.0625) Enddo n = V5dwrite( itime, varV5D(ivarREGRID), G ) If (n.Eq.0) Then Call Exit(1) Endif Endif 44 Continue Endif C 1D FIELDS Else If (ndim.Eq.1) Then Read(iREGRIDunit) f1d Endif Read(iREGRIDunit) flag ENDDO ENDDO C AT THIS POINT ALL REQUESTED DATES WERE PROCESSED Print * Print *, $ '... SUCCESSFUL PROCESS OF ALL REQUESTED DATES !!!' C ERROR IN FINDING REQUESTED DATES ............ Goto 10001 10000 Continue Print *,' ' Print *,'ERROR !!! EOF BEFORE ALL REQUESTED DATES WERE FOUND' Stop 10000 10001 Continue C TIME-MEANS AND WRITE If (mean_Tm.Gt.0) Then Do itime=1,numtimes n = V5dwrite( itime, mean_Tm, Tm ) If (n.Eq.0) Then Call Exit(1) Endif Enddo Endif If (mean_Um.Gt.0) Then Do itime=1,numtimes n = V5dwrite( itime, mean_Um, Um ) If (n.Eq.0) Then Call Exit(1) Endif Enddo Endif If (mean_Vm.Gt.0) Then Do itime=1,numtimes n = V5dwrite( itime, mean_Vm, Vm ) If (n.Eq.0) Then Call Exit(1) Endif Enddo Endif If (mean_Hm.Gt.0) Then Do itime=1,numtimes n = V5dwrite( itime, mean_Hm, Hm ) If (n.Eq.0) Then Call Exit(1) Endif Enddo Endif If (mean_RHm.Gt.0) Then Do itime=1,numtimes n = V5dwrite( itime, mean_RHm, RHm ) If (n.Eq.0) Then Call Exit(1) Endif Enddo Endif C CLOSE THE V5D FILE AND EXIT Close(iREGRIDunit) n = V5dclose() If (n.Eq.0) Then Print*,'V5D CLOSE FAILED' Call Exit(1) Else Print*,'V5D CLOSE OK' Call Exit(0) Endif RETURN END C ================================================================================================== C ================================================================================================== SUBROUTINE SMOOTH(HT,NXS,NYS,NZS,MAXZ,HSX4,HSY4) C ============================================================== C C SUBRUTINA SUAVIZACION DE CUARTO ORDEN PARA VALORES DIGITALIZADOS C EN UNA MALLA C REMITIDA AMABLEMENTE POR DALE R. DURRAN (DEPART. OF METEOROLOGY C UNIVERSITY OF UTAH. SALT LAKE CITY. UT 84112) C C MODIFICADA POR C. RAMIS (1991) C =============================================================== Implicit none Include "/usr/local/vis5d/src/v5df.h" integer k,j,i,nxs,nys,nzs,maxz REAL HSX4,HSY4,HSX2,HSY2 REAL HT(NXS,NYS,NZS), HTS(MAXROWS,MAXCOLUMNS) C C C COEFICIENTES PARA FILTRO DE 4 ORDEN C C HSX2=-4.*HSX4 HSY2=-4.*HSY4 C C ========================================================= C HSX4=HSY4=.06 (.0625 MEJOR, ELIMINA 2DX EN UNA PASADA) C ========================================================= C C C SUAVIZACION DE LOS CAMPOS C DO 1000 K=1,MAXZ DO 204 J=1,NYS DO 203 I=3,NXS-2 HTS(I,J)=HT(I,J,K)-HSX4*(HT(I+2,J,K)+HT(I-2,J,K) * -4.*(HT(I+1,J,K)+HT(I-1,J,K))+6.*HT(I,J,K)) 203 CONTINUE HTS(2,J)=HT(2,J,K)-HSX2*(HT(3,J,K)+HT(1,J,K)-2.*HT(2,J,K)) HTS(NXS-1,J)=HT(NXS-1,J,K)-HSX2*(HT(NXS,J,K)+HT(NXS-2,J,K) * -2.*HT(NXS-1,J,K)) 204 CONTINUE DO 206 I=1,NXS DO 205 J=3,NYS-2 HTS(I,J)=HTS(I,J)-HSY4*(HT(I,J+2,K)+HT(I,J-2,K) * -4.*(HT(I,J+1,K)+HT(I,J-1,K))+6.*HT(I,J,K)) 205 CONTINUE HTS(I,2)=HTS(I,2)-HSY2*(HT(I,3,K)+HT(I,1,K)-2.*HT(I,2,K)) HTS(I,NYS-1)=HTS(I,NYS-1)-HSY2*(HT(I,NYS,K)+HT(I,NYS-2,K) * -2.*HT(I,NYS-1,K)) 206 CONTINUE C C REASIGNACION DE VALORES A LA MATRIZ ORIGINAL C DO 208 I=2,NXS-1 DO 208 J=2,NYS-1 HT(I,J,K)=HTS(I,J) 208 CONTINUE 1000 CONTINUE RETURN END C ==================================================================================================