C ==================================================================================================
      PROGRAM MM5p_to_v5d

C -------------------------------------------------------------------------------------------------------
C This code is used to transform pressure-coordinate MM5p output to V5D format

C COMPILATION: make -f MM5p_to_v5d.f.m

C EXECUTION: MM5p_to_v5d pars.MM5p 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.MM5p)
       Integer nvarV5D
       Parameter(nvarV5D=54)
       Integer imonth,julian,days_month(12)
       Integer NY(100),FILT(100)
       Integer varV5D(100),filtV5D(100)
       Real sumV5D(100),facV5D(100)                
       Integer ivarV5D,iMM5punit,nvarMM5p,
     $         ivarMM5p,idimMM5p(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),NAMEMM5p(100) 
       Character*24 date(100),date_means(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),        
     $      Tm(MAXROWS,MAXCOLUMNS,MAXLEVELS),
     $      Um(MAXROWS,MAXCOLUMNS,MAXLEVELS),
     $      Vm(MAXROWS,MAXCOLUMNS,MAXLEVELS),
     $      Hm(MAXROWS,MAXCOLUMNS,MAXLEVELS),
     $      RHm(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)


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, 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 parsMM5pname (SEE FOR USER CHOICES)
       Open(11,file=parsMM5pname,status='old')
          Read(11,*)
          Read(11,*)
          Read(11,*)
          Read(11,*) MM5pname
          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 MM5p 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
   
          If (BHI(1,1).ne.8) 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,8)-1

          nvarMM5p=0       
          Read(iMM5punit) flag
          Do While (flag.Eq.1)           
            Read(iMM5punit) 
     $           ndim,start_index,end_index,xtime,staggering,
     $           ordering,current_date,name,unit,description
            nvarMM5p=nvarMM5p+1
            idimMM5p(nvarMM5p)=ndim
            NAMEMM5p(nvarMM5p)=name
            If (ndim.Eq.3) Then
              Read(iMM5punit) 
     $        (((f3d(ir,ic,ilev),ir=1,nr+1),ic=1,nc+1),ilev=1,maxnl+1)
            Else If (ndim.Eq.2) Then
              Read(iMM5punit) ((f2d(ir,ic),ir=1,nr+1),ic=1,nc+1)           
            Else If (ndim.Eq.1) Then
              Read(iMM5punit) (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(iMM5punit) flag
          Enddo         
       
       Close(iMM5punit)

       numvars=0 
       Do ivarMM5p=1,nvarMM5p
          varV5D(ivarMM5p)=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 ivarMM5p=1,nvarMM5p

             If (NAMEV5D(ivarV5D).Eq.NAMEMM5p(ivarMM5p)) Then
                numvars=numvars+1
                varV5D(ivarMM5p)=numvars
                filtV5D(ivarMM5p)=FILT(ivarV5D)
                                       
                If (idimMM5p(ivarMM5p).Eq.3) Then
                   nl(numvars)=maxnl  
                Else
                   nl(numvars)=1
                Endif 
                varname(numvars)=NAMEMM5p(ivarMM5p)       
		If (varname(numvars).Eq.'RAD TEND  ') varname(numvars)='RADTEND	  '
		If (varname(numvars).Eq.'RAIN CON  ') varname(numvars)='RAINCON   '	
		If (varname(numvars).Eq.'RAIN NON  ') varname(numvars)='RAINNON   '
		If (varname(numvars).Eq.'RES TEMP  ') varname(numvars)='RESTEMP   '
		If (varname(numvars).Eq.'LAND USE  ') varname(numvars)='LANDUSE   '	
		If (varname(numvars).Eq.'GROUND T  ') varname(numvars)='GROUNDT   '	
		If (varname(numvars).Eq.'PBL HGT   ') varname(numvars)='PBLHGT    '
		If (varname(numvars).Eq.'SOIL T 1  ') varname(numvars)='SOILT1    '
		If (varname(numvars).Eq.'SOIL T 2  ') varname(numvars)='SOILT2    '
		If (varname(numvars).Eq.'SOIL T 3  ') varname(numvars)='SOILT3    '	
		If (varname(numvars).Eq.'SOIL T 4  ') varname(numvars)='SOILT4    '
		If (varname(numvars).Eq.'SOIL T 5  ') varname(numvars)='SOILT5    '
		If (varname(numvars).Eq.'SOIL T 6  ') varname(numvars)='SOILT6    '	
		If (varname(numvars).Eq.'RH SFC    ') varname(numvars)='RHSFC     '	
                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      (NAMEMM5p(ivarMM5p).Eq.'U        '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'V        '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'UST      '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'U10      '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'V10      ') Then
                                         units(numvars)='m/s'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.
                Else If (NAMEMM5p(ivarMM5p).Eq.'W        ') Then
                                         units(numvars)='cm/s'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=100.
                Else If (NAMEMM5p(ivarMM5p).Eq.'RAD TEND ') Then
                                         units(numvars)='C/day'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'Q        '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'CLW      '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'RNW      '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'ICE      '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SNOW     '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'GRAUPEL  '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'Q2       ') Then                                      
                                         units(numvars)='g/kg'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1000.
                Else If (NAMEMM5p(ivarMM5p).Eq.'RH       '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'RH SFC   ') Then
                                         units(numvars)='%'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.
                Else If (NAMEMM5p(ivarMM5p).Eq.'NCI      ') Then
                                         units(numvars)='numb/m3'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.
                Else If (NAMEMM5p(ivarMM5p).Eq.'T        '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'GROUND T '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'RES TEMP '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SOIL T 1 '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SOIL T 2 '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SOIL T 3 '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SOIL T 4 '.Or. 
     $                   NAMEMM5p(ivarMM5p).Eq.'SOIL T 5 '.Or.               
     $                   NAMEMM5p(ivarMM5p).Eq.'SOIL T 6 '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'TSEASFC  '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'T2       ') Then                                      
                                         units(numvars)='C'
                                         sumV5D(ivarMM5p)=-273.16
                                         facV5D(ivarMM5p)=1.
                Else If (NAMEMM5p(ivarMM5p).Eq.'H        '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'TERRAIN  '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'PBL HGT  ') Then
                                         units(numvars)='m'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.
                Else If (NAMEMM5p(ivarMM5p).Eq.'TKE      ') Then
                                         units(numvars)='J/kg'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'PSFC     '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'PSEALVLC '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'PSEALVLD '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'PSTARCRS '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'PP       ') Then
                                         units(numvars)='mb'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=0.01
                Else If (NAMEMM5p(ivarMM5p).Eq.'MAPFACCR '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'MAPFACDT '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SNOWCOVR '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'REGIME   ') Then                                      
                                         units(numvars)='no dim.'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'CORIOLIS ') Then
                                         units(numvars)='1/s'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'LATITCRS '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'LATITDOT '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'LONGICRS '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'LONGIDOT ') Then                                      
                                         units(numvars)='degrees'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                               facV5D(ivarMM5p)=1.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'LAND USE ') Then
                                         units(numvars)='categ.'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'RAIN NON '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'RAIN CON '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'RAIN TOT ') Then
                                         units(numvars)='mm'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=10.                                                       
                Else If (NAMEMM5p(ivarMM5p).Eq.'SHFLUX   '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'LHFLUX   '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SWDOWN   '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'LWDOWN   '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'SWOUT    '.Or.
     $                   NAMEMM5p(ivarMM5p).Eq.'LWOUT    ') Then                                      
                                         units(numvars)='W/m2'
                                         sumV5D(ivarMM5p)=0.
                                         facV5D(ivarMM5p)=1.                                                              
                Endif                  

                Print 578, numvars,NAMEMM5p(ivarMM5p),
     $                     filtV5D(ivarMM5p),units(numvars)
              
                If (time_means.Eq.1) Then
                   If (NAMEMM5p(ivarMM5p).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(ivarMM5p),units(numvars)     
                   Else If (NAMEMM5p(ivarMM5p).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(ivarMM5p),units(numvars)     
                   Else If (NAMEMM5p(ivarMM5p).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(ivarMM5p),units(numvars)     
                   Else If (NAMEMM5p(ivarMM5p).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(ivarMM5p),units(numvars)     
                   Else If (NAMEMM5p(ivarMM5p).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(ivarMM5p),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-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 *, '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,G,f3d,f2d,f1d,
     $               MM5pname,nvarMM5p,
     $               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 MM5p OUTPUT .........
9003   Continue   
       Print *,'  '                                                             
       Print *,'ERROR: NOT A MM5p OUTPUT'                      
       Stop 9003                          
       
       
      END
C ==================================================================================================      


C ==================================================================================================
      SUBROUTINE CONVERT (nr,nc,maxnl,G,f3d,f2d,f1d,
     $                    MM5pname,nvarMM5p,
     $                    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,
     $           ivarMM5p,nvarMM5p,
     $           iMM5punit,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 NAMEMM5p
         Real sumV5D(100),facV5D(100)                                                             
         Character*24 date(100),date_means(100)
         Character*100 MM5pname             

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 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 MM5p OUTPUT FILE
         iMM5punit=20
         Open(unit=iMM5punit,file=MM5pname,form='unformatted')
         Rewind(iMM5punit)    
         Read(iMM5punit) flag
         Read(iMM5punit) 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 MM5pname, and PROCESSING 
C AND WRITTING IN v5dname if time period was requested)
          DO WHILE (itime.Lt.numtimes.Or.itime_means.Lt.numtimes_means)

          
          Read(iMM5punit,err=10000) flag 
          ivarMM5p=0                      
          DO WHILE (flag.Eq.1)           
             Read(iMM5punit) 
     $            ndim,start_index,end_index,xtime,staggering,
     $            ordering,current_date,name,unit,description
             ivarMM5p=ivarMM5p+1
             NAMEMM5p=name
             

             If (ivarMM5p.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(iMM5punit) f3d

                   If (iprocess_time.Eq.0.And.iprocess_time_means.Eq.0) 
     $             Goto 11
                                       
                   If (varV5D(ivarMM5p).Gt.0) Then

                      If (sumV5D(ivarMM5p).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(ivarMM5p) 
                            Enddo   
                         Enddo 
                         Enddo

                      Else If (facV5D(ivarMM5p).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(ivarMM5p)                        
                            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(ivarMM5p)
                         Call Smooth(G,nr,nc,maxnl,maxnl,0.0625,0.0625)
                      Enddo

                      If (iprocess_time.Eq.0) Goto 111
     
                      n = V5dwrite( itime, varV5D(ivarMM5p), G )
                      If (n.Eq.0) Then
                         Call Exit(1)
                      Endif
                    
111                   Continue                      

                      If (iprocess_time_means.Eq.0) Goto 1111
                                            
                      If (NAMEMM5p.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 (NAMEMM5p.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(iMM5punit) f3d

                   If (iprocess_time.Eq.0.And.iprocess_time_means.Eq.0) 
     $             Goto 22
                                      
                   If (varV5D(ivarMM5p).Gt.0) Then

                      If (sumV5D(ivarMM5p).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(ivarMM5p) 
                            Enddo   
                         Enddo 
                         Enddo

                      Else If (facV5D(ivarMM5p).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(ivarMM5p)                        
                            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(ivarMM5p)
                         Call Smooth(G,nr,nc,maxnl,maxnl,0.0625,0.0625)
                      Enddo
 
                      If (NAMEMM5p.Eq.'Q        '.Or.
     $                    NAMEMM5p.Eq.'Q2       '.Or.
     $                    NAMEMM5p.Eq.'CLW      '.Or.  
     $                    NAMEMM5p.Eq.'RNW      '.Or.  
     $                    NAMEMM5p.Eq.'ICE      '.Or.  
     $                    NAMEMM5p.Eq.'SNOW     '.Or.
     $                    NAMEMM5p.Eq.'GRAUPEL  '.Or.  
     $                    NAMEMM5p.Eq.'NCI      '.Or.
     $                    NAMEMM5p.Eq.'RH       '.Or.   
     $                    NAMEMM5p.Eq.'RH SFC   ')    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 (NAMEMM5p.Eq.'RH       '.Or.
     $                    NAMEMM5p.Eq.'RH SFC   ')    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(ivarMM5p), G )
                      If (n.Eq.0) Then
                         Call Exit(1)
                      Endif
               
222                   Continue                      

                      If (iprocess_time_means.Eq.0) Goto 2222 
                                         
                      If (NAMEMM5p.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 (NAMEMM5p.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 (NAMEMM5p.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(iMM5punit) f2d                   

                   If (iprocess_time.Eq.0) Goto 33
                                                                    
                   If (varV5D(ivarMM5p).Gt.0) Then

                      If (sumV5D(ivarMM5p).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(ivarMM5p) 
                            Enddo   
                         Enddo 

                      Else If (facV5D(ivarMM5p).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(ivarMM5p)                        
                            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(ivarMM5p)
                         Call Smooth(G,nr,nc,maxnl,1,0.0625,0.0625)
                      Enddo

                      n = V5dwrite( itime, varV5D(ivarMM5p), G )
                      If (n.Eq.0) Then
                         Call Exit(1)
                      Endif

                   Endif      

33                 Continue
                                                                                                                                                   
C Crosses (2D)                  
                Else If (staggering.Eq.'C   ') Then
 
                   Read(iMM5punit) f2d
 
                   If (iprocess_time.Eq.0) Goto 44
                       
                   If (varV5D(ivarMM5p).Gt.0) Then

                      If (sumV5D(ivarMM5p).Ne.0.) Then   
                         Do ir=1,nr
                            Do ic=1,nc
                               G(nr+1-ir,ic,1)=
     $                                        f2d(ir,ic)   
     $                                       +sumV5D(ivarMM5p) 
                            Enddo   
                         Enddo 

                      Else If (facV5D(ivarMM5p).Ne.1.) Then
                         Do ir=1,nr
                            Do ic=1,nc                                                                                   
                               G(nr+1-ir,ic,1)=
     $                                        f2d(ir,ic)
     $                                       *facV5D(ivarMM5p)                        
                            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(ivarMM5p)
                         Call Smooth(G,nr,nc,maxnl,1,0.0625,0.0625)
                      Enddo

                      n = V5dwrite( itime, varV5D(ivarMM5p), G )
                      If (n.Eq.0) Then
                         Call Exit(1)
                      Endif

                   Endif 

44                 Continue
                                                     
                Endif 

C 1D FIELDS
             Else If (ndim.Eq.1) Then 
                
                Read(iMM5punit) f1d
                
             Endif
             
                       
             Read(iMM5punit) 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(iMM5punit)

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

