C ================================================================================================== PROGRAM Adj_to_v5d C ------------------------------------------------------------------------------------------------------- C This code is used to transform sigma-coordinate MM5 Adj output to V5D format C ATTENTION: VIS5D files created using standard UIB utilities, are stored on cross grid. The v5d files generated with this program C uses DOTS to consistently build MM5 perturbations from sensitivity fields. C C ****** THIS PROGRAM GENERATES DOTS FIELDs ******* C C COMPILATION: make C EXECUTION: adj2v5d-MM5 pars.adj output.v5d c NEEDED FILES: The program uses the files GRAD.D, GLBCT.00[1,2,...] (Adjoint model output) and MMINPUT_DOMAIN1 C See v5df.h for parameters MAXVARS, MAXTIMES, MAXROWS, MAXCOLUMNS, MAXLEVELS, MISSING, IMISSING C See comments in the code. C ------------------------------------------------------------------------------------------------------- Implicit none Include "v5df.h" Integer v5dsetlowlev,iargc Integer imonth,julian,days_month(12) Integer iMM5punit, $ ivar,ilev, $ intSS,YYYYo,YYYYf,MMo,MMf,DDo,DDf, $ HHo,HHf,MIo,MIf,SSo,SSf,num60SS,inum60SS, $ n,i,maxnl Real a,conv,xn,xlatc,phi1,phic,cnst,yyc,ymcentri,centri,x,y Real*8 YYYYMMDDHHMISSo,YYYYMMDDHHMISSf Character*24 date(100) Character*100 parsMM5pname,MM5pname,v5dname C VARIABLES DESCRIBING MM5p 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) Integer ir,ic,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),sigma(MAXLEVELS),ps,pt parameter(MM5pname='MMINPUT_DOMAIN1') integer IEXMS,IICE 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 / 4 / 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 8395 Print *,'Error: two filename arguments are needed.' Print *,'Usage: adj2v5d parsfile file.v5d' Print *,'The program uses the files GRAD.D, GLBCT.00[1,2,...] and & MMINPUT_DOMAIN1' Call Exit(1) Else Call Getarg(1, parsMM5pname) Call Getarg(2, v5dname) Endif Print * Print * Print 10,'Input file (pars MM5p): ', parsMM5pname Print 10,'Output file (V5D): ', v5dname 10 Format (A,A) C READ MM5p PARAMETERS FROM pars file (SEE FOR USER CHOICES) Open(11,file=parsMM5pname,status='old',readonly,ERR=8395) Read(11,*) Read(11,*) Read(11,*) Read(11,*)IEXMS,IICE Read(11,*) C READ BEGIN DATE Read(11,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') $ YYYYo,MMo,DDo,HHo,MIo,SSo C C Convert to julian day 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 CC Convert to v5d date format YYYYMMDDHHMISSo= $ YYYYo*10000000000.d0+MMo*100000000.d0+ $ DDo*1000000.d0+HHo*10000.d0+MIo*100.d0+SSo Read(11,*) C READ END DATE Read(11,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') $ YYYYf,MMf,DDf,HHf,MIf,SSf C C Convert to v5d date format YYYYMMDDHHMISSf= $ YYYYf*10000000000.d0+MMf*100000000.d0+ $ DDf*1000000.d0+HHf*10000.d0+MIf*100.d0+SSf Read(11,*) C READ Interval of the timesteps (in seconds) Read(11,*) intSS if (intSS.eq.0) then numtimes=1 dates(numtimes)=(YYYYo-(YYYYo/100)*100)*1000+julian times(numtimes)=HHo*10000+MIo*100+SSo else num60SS=intSS/60 C CALCULATE numtimesteps 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 endif ! Of intSS=0 Close(11) C END pars FILE READ ** END pars FILE READ ** END pars FILE READ ** END pars FILE READ ** C END pars FILE READ ** END pars FILE READ ** END pars FILE READ ** END pars FILE READ ** C OPEN MMINPUT DATAFILE (MM5pname) HERE AND READ ITS HEADER INFORMATION TO C INITIALIZE SOME OF THE ABOVE V5D PARAMETERS iMM5punit=20 Open(iMM5punit,file=MM5pname,form='unformatted') Rewind(iMM5punit) Read(iMM5punit,err=9001,end=9002) flag Read(iMM5punit,err=9001,end=9002) BHI,BHR,BHIC,BHRC C Identify the MM5 file type If (BHI(1,1).ne.5) Go To 9003 C Refuse nested domains If (BHI(15,1).Ne.0) Then print*,'NESTED DOMAINS NOT SUPPORTED. CHANGE PROGRAM' C nr=BHI(16,1)-1 C nc=BHI(17,1)-1 Else nr=BHI(5,1) nc=BHI(6,1) Endif C Define Number of Sigma levels maxnl=BHI(12,BHI(1,1)) C Define the value of vert_args. The strategy is to define the levels as P levels, although C they aren't. The p value assigned to each sigma level will be the one obtained by the expression of C sigma definition, using Ptop and Ps from the 1D pressure field of MMINPUT (top and bottom levels C of RAWINS file). vertical=3 Read(iMM5punit) flag Do While (flag.Eq.1) Read(iMM5punit) $ ndim,start_index,end_index,xtime,staggering, $ ordering,current_date,name,unit,description If (ndim.Eq.3) Then Read(iMM5punit) $ (((f3d(ir,ic,ilev),ir=1,nr),ic=1,nc),ilev=1,maxnl) Else If (ndim.Eq.2) Then Read(iMM5punit) ((f2d(ir,ic),ir=1,nr),ic=1,nc) Else If (ndim.Eq.1) Then Read(iMM5punit) (f1d(ilev),ilev=1,end_index(1)) C Keep sigma levels values for later definition of P as vertical coord in V5d if (name.eq.'SIGMAH ') then do ilev=1,end_index(1) sigma(ilev)=f1d(ilev) enddo endif If (name.Eq.'PRESSURE') Then ps=f1d(1) pt=f1d(end_index(1)) Endif Endif Read(iMM5punit) flag Enddo Close(iMM5punit) Do ilev=1,maxnl vert_args(maxnl-ilev+1)=(sigma(ilev)*(ps-pt)+pt)/100. Enddo C END of VERTICAL levels DEFINITION C Setup of Variable names, units, and low_lev numvars=16+1 ! Grads+ Sigma level value varname(1)='UA' units(1)='R/(m/s)' varname(2)='UAsmo' units(2)=units(1) varname(3)='VA' units(3)=units(1) varname(4)='VAsmo' units(4)=units(1) varname(5)='WA' units(5)=units(2) varname(6)='WAsmo' units(6)=units(2) varname(7)='TA' units(7)='R/(C)' varname(8)='TAsmo' units(8)=units(7) varname(9)='QVA' units(9)='R/(Kg/Kg)' varname(10)='QVAsmo' units(10)=units(9) varname(11)='PPA' units(11)='R/(Pa)' varname(12)='PPAsmo' units(12)=units(11) varname(13)='QCA' units(13)=units(5) varname(14)='QCAsmo' units(14)=units(5) varname(15)='QRA' units(15)=units(5) varname(16)='QRAsmo' units(16)=units(5) varname(17)='Sigma' units(17)='no dim.' do i=1,numvars lowlev(i)=0 nl(i)=maxnl 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-5.1 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 *, 'Number of timesteps = ',numtimes Print *, 'Vertical coordinate = ',vertical,' (pressure)' Print *, 'Map projection = ', $ projection,' (Lambert conformal)' Print *, 'Compressmode = ',compressmode Print * C CALL SUBROUTINE CONVERT TO WRITE MM5p OUTPUT AS V5D Call CONVERT (nr,nc,maxnl,nl,G, $ numvars, $ numtimes,sigma,varname,IEXMS,IICE) 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 MM5p OUTPUT ......... 9003 Continue Print *,' ' Print *,'ERROR: NOT A MM5 INPUT file (MMINPUT)' Stop 9003 END C ================================================================================================== C ================================================================================================== SUBROUTINE CONVERT (nr,nc,maxnl,nl,G, $ numvars, $ numtimes,sigma,varname,IEXMS,IICE) Implicit none Include "v5df.h" C LOCAL VARIABLES Integer nr,nc,maxnl,n,ir,ic,il, $ numvars,numtimes,itime Integer nl(numvars),ivar,iu,nsmoo(numvars) Real G(nr,nc,maxnl) C The Adj output variables are written to the GRAD.D file as real*8. Do not modify! real*8 UA(nr,nc,maxnl),VA(nr,nc,maxnl), & WA(nr,nc,maxnl+1),TA(nr,nc,maxnl),QVA(nr,nc,maxnl), & PPA(nr,nc,maxnl),QCA(nr,nc,maxnl),QRA(nr,nc,maxnl), & TGA(nr,nc) real*8 & UWB9(nr,maxnl,5),UEB9(nr,maxnl,5), & USB9(nc,maxnl,5),UNB9(nc,maxnl,5), & VWB9(nr,maxnl,5),VEB9(nr,maxnl,5), & VSB9(nc,maxnl,5),VNB9(nc,maxnl,5), & TWB9(nr,maxnl,5),TEB9(nr,maxnl,5), & TSB9(nc,maxnl,5),TNB9(nc,maxnl,5), & WWB9(nr,maxnl+1,5),WEB9(nr,maxnl+1,5), & WSB9(nc,maxnl+1,5),WNB9(nc,maxnl+1,5), & QWB9(nr,maxnl,5),QEB9(nr,maxnl,5), & QSB9(nc,maxnl,5),QNB9(nc,maxnl,5), & PPWB9(nr,maxnl,5),PPEB9(nr,maxnl,5), & PPSB9(nc,maxnl,5),PPNB9(nc,maxnl,5), & QCWB9(nr,maxnl,5),QCEB9(nr,maxnl,5), & QCSB9(nc,maxnl,5),QCNB9(nc,maxnl,5), & QRWB9(nr,maxnl,5),QREB9(nr,maxnl,5), & QRSB9(nc,maxnl,5),QRNB9(nc,maxnl,5), & QIWB9(1,1,5),QIEB9(1,1,5), & QISB9(1,1,5),QINB9(1,1,5), & QNIWB9(1,1,5),QNIEB9(1,1,5), & QNISB9(1,1,5),QNINB9(1,1,5) character*9 FILENAME integer*8 IEXMS,IICE,nr8,nc8,maxnl8,u8,numtim8 real sigma(maxnl) Character*10 varname(numvars) character*6 adjname C Some adjoint parameters parameter(adjname='GRAD.D') logical inwards,domainij C COnvert from single to double precision variables nr8=nr nc8=nc maxnl8=maxnl u8=1 numtim8=numtimes C Temprary, we fix nsmoo for all vars do ivar=1,numvars nsmoo(ivar)=5 enddo C WRITE to V5D FILE Do itime = 1,numtimes if (itime.eq.1) then C OPEN Adj OUTPUT FILE write(*,'(A)')'Reading IC gradients from GRAD.D' Open(20,file=adjname,form='unformatted',status='OLD', &readonly) read(20)UA,VA,WA,TA,QVA,PPA,QCA,QRA !GRAD.D file-type close(20) Do ivar=1,numvars Do il=1,nl(ivar) Do ir=1,nr Do ic=1,nc if(ivar.eq.1) G(ir,ic,il)=real(UA(nr-ir+1,ic,nl(ivar)-il+1)) if(ivar.eq.3) G(ir,ic,il)=real(VA(nr-ir+1,ic,nl(ivar)-il+1)) if(ivar.eq.5) G(ir,ic,il)=real(WA(nr-ir+1,ic,nl(ivar)-il+2)) if(ivar.eq.7) G(ir,ic,il)=real(TA(nr-ir+1,ic,nl(ivar)-il+1)) if(ivar.eq.9) G(ir,ic,il)=real(QVA(nr-ir+1,ic,nl(ivar)-il+1)) if(ivar.eq.11) G(ir,ic,il)=real(PPA(nr-ir+1,ic,nl(ivar)-il+1)) if(ivar.eq.13) G(ir,ic,il)=real(QCA(nr-ir+1,ic,nl(ivar)-il+1)) !GRAD file if(ivar.eq.15) G(ir,ic,il)=real(QRA(nr-ir+1,ic,nl(ivar)-il+1)) if(ivar.eq.17) G(ir,ic,il)=sigma(nl(ivar)-il+1) Enddo Enddo Enddo C Smooth all variables but Sigma if ((ivar/2)*2.eq.ivar) then print*,'Smoothing ',varname(ivar),nsmoo(ivar),' times' do n=1,nsmoo(ivar) call smooth(G,nr,nc,maxnl,maxnl,0.0625,0.0625) enddo endif C Write fields to V5D FILE n=v5dwrite(itime,ivar,G) C Do 37 il=1,maxnl C Do 37 ir=1,nr C Do 37 ic=1,nc C G(ir,ic,il)=0. C37 Enddo Enddo !ivar do else C Open Gradients at the boundary files print* write(*,'(A,I3.3)')'Reading boundary gradients from GLBCT.', &itime-1 write(filename(1:9),'(A6,I3.3)')'GLBCT.',itime-1 CALL BDY_STORE(-1,FILENAME,itime-1, 1 UWB9, UEB9, USB9, UNB9, 2 VWB9, VEB9, VSB9, VNB9, 3 TWB9, TEB9, TSB9, TNB9, 4 WWB9, WEB9, WSB9, WNB9, 5 QWB9, QEB9, QSB9, QNB9, 6 PPWB9, PPEB9, PPSB9, PPNB9, 7 QCWB9, QCEB9, QCSB9, QCNB9, 8 QRWB9, QREB9, QRSB9, QRNB9, 9 QIWB9, QIEB9, QISB9, QINB9, & QNIWB9,QNIEB9,QNISB9,QNINB9, 1 nr, nc, maxnl,nr, nc, maxnl,maxnl+1, 5, 2 nr, nc,maxnl,1, 1, 1, 5, 3 IEXMS,IICE) CVHS How are boundary vertical layers ordered (inward, outwards or following domain I, J)?. inwards=.FALSE. domainij=.TRUE. C inwards=.TRUE. C domainij=.FALSE. C Compose all boundary fields into G matrix: First E and W and then overwrite values at the corners C with the N and S good ones. do ivar=1,numvars do il=1,nl(ivar) C West and East boundaries do ir=1,nr do ic=5,1,-1 !Sponge region around the boundary if (inwards) then if(ivar.eq. 1.or.ivar.eq. 2) &G(ir, ic,il)= UWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 1.or.ivar.eq. 2) &G(ir,nc+1-ic,il)= UEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 3.or.ivar.eq. 4) &G(ir, ic,il)= VWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 3.or.ivar.eq. 4) &G(ir,nc+1-ic,il)= VEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 5.or.ivar.eq. 6) &G(ir, ic,il)= WWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 5.or.ivar.eq. 6) &G(ir,nc+1-ic,il)= WEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 7.or.ivar.eq. 8) &G(ir, ic,il)= TWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 7.or.ivar.eq. 8) &G(ir,nc+1-ic,il)= TEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 9.or.ivar.eq.10) &G(ir, ic,il)= QWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 9.or.ivar.eq.10) &G(ir,nc+1-ic,il)= QEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.11.or.ivar.eq.12) &G(ir, ic,il)=PPWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.11.or.ivar.eq.12) &G(ir,nc+1-ic,il)=PPEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.13.or.ivar.eq.14) &G(ir, ic,il)=QCWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.13.or.ivar.eq.14) &G(ir,nc+1-ic,il)=QCEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.15.or.ivar.eq.16) &G(ir, ic,il)=QRWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.15.or.ivar.eq.16) &G(ir,nc+1-ic,il)=QREB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.17) G(ir,ic,il)=sigma(nl(ivar)-il+1) endif ! Inwards if (domainij) then if(ivar.eq. 1.or.ivar.eq. 2) &G(ir, ic,il)= UWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 1.or.ivar.eq. 2) &G(ir,nc-5+ic,il)= UEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 3.or.ivar.eq. 4) &G(ir, ic,il)= VWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 3.or.ivar.eq. 4) &G(ir,nc-5+ic,il)= VEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 5.or.ivar.eq. 6) &G(ir, ic,il)= WWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 5.or.ivar.eq. 6) &G(ir,nc-5+ic,il)= WEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 7.or.ivar.eq. 8) &G(ir, ic,il)= TWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 7.or.ivar.eq. 8) &G(ir,nc-5+ic,il)= TEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 9.or.ivar.eq.10) &G(ir, ic,il)= QWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq. 9.or.ivar.eq.10) &G(ir,nc-5+ic,il)= QEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.11.or.ivar.eq.12) &G(ir, ic,il)=PPWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.11.or.ivar.eq.12) &G(ir,nc-5+ic,il)=PPEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.13.or.ivar.eq.14) &G(ir, ic,il)=QCWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.13.or.ivar.eq.14) &G(ir,nc-5+ic,il)=QCEB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.15.or.ivar.eq.16) &G(ir, ic,il)=QRWB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.15.or.ivar.eq.16) &G(ir,nc-5+ic,il)=QREB9(nr-ir+1,nl(ivar)-il+1,ic) if(ivar.eq.17) G(ir,ic,il)=sigma(nl(ivar)-il+1) endif ! domainij enddo ! ic (sponge) enddo !ir C North and South Boundaries do ic=1,nc do ir=5,1,-1 !Sponge region if(inwards) then if(ivar.eq. 1.or.ivar.eq. 2) &G(ir, ic,il)= UNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 1.or.ivar.eq. 2) &G(nr+1-ir,ic,il)= USB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 3.or.ivar.eq. 4) &G(ir, ic,il)= VNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 3.or.ivar.eq. 4) &G(nr+1-ir,ic,il)= VSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 5.or.ivar.eq. 6) &G(ir, ic,il)= WNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 5.or.ivar.eq. 6) &G(nr+1-ir,ic,il)= WSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 7.or.ivar.eq. 8) &G(ir, ic,il)= TNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 7.or.ivar.eq. 8) &G(nr+1-ir,ic,il)= TSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 9.or.ivar.eq.10) &G(ir, ic,il)= QNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 9.or.ivar.eq.10) &G(nr+1-ir,ic,il)= QSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.11.or.ivar.eq.12) &G(ir, ic,il)=PPNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.11.or.ivar.eq.12) &G(nr+1-ir,ic,il)=PPSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.13.or.ivar.eq.14) &G(ir, ic,il)=QCNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.13.or.ivar.eq.14) &G(nr+1-ir,ic,il)=QCSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.15.or.ivar.eq.16) &G(ir, ic,il)=QRNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.15.or.ivar.eq.16) &G(nr+1-ir,ic,il)=QRSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.17) G(ir,ic,il)=sigma(nl(ivar)-il+1) endif ! Inwards N & S if(domainij) then if(ivar.eq. 1.or.ivar.eq. 2) &G(6-ir, ic,il)= UNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 1.or.ivar.eq. 2) &G(nr+1-ir,ic,il)= USB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 3.or.ivar.eq. 4) &G(6-ir, ic,il)= VNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 3.or.ivar.eq. 4) &G(nr+1-ir,ic,il)= VSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 5.or.ivar.eq. 6) &G(6-ir, ic,il)= WNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 5.or.ivar.eq. 6) &G(nr+1-ir,ic,il)= WSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 7.or.ivar.eq. 8) &G(6-ir, ic,il)= TNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 7.or.ivar.eq. 8) &G(nr+1-ir,ic,il)= TSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 9.or.ivar.eq.10) &G(6-ir, ic,il)= QNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq. 9.or.ivar.eq.10) &G(nr+1-ir,ic,il)= QSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.11.or.ivar.eq.12) &G(6-ir, ic,il)=PPNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.11.or.ivar.eq.12) &G(nr+1-ir,ic,il)=PPSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.13.or.ivar.eq.14) &G(6-ir, ic,il)=QCNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.13.or.ivar.eq.14) &G(nr+1-ir,ic,il)=QCSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.15.or.ivar.eq.16) &G(6-ir, ic,il)=QRNB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.15.or.ivar.eq.16) &G(nr+1-ir,ic,il)=QRSB9(ic,nl(ivar)-il+1,ir) if(ivar.eq.17) G(ir,ic,il)=sigma(nl(ivar)-il+1) endif ! Domainij N & S enddo !ir enddo! ic enddo ! il n=v5dwrite(itime,ivar,G) C Do 38 il=1,maxnl C Do 38 ir=1,nr C Do 38 ic=1,nc C G(ir,ic,il)=0. C38 Enddo enddo ! ivar endif !itime if enddo ! itime do C AT THIS POINT ALL REQUESTED DATES WERE PROCESSED Print * Print *, $ '... SUCCESSFUL PROCESS OF ALL REQUESTED DATES !!!' 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 "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 ==================================================================================================