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 ==================================================================================================

