C I_RAWp2.f
C
C VIS-5D function used to SUBSTRACT COMPONENT 2 OF PERTURBATION balanced fields
C I_Hp2, I_Tp2, I_Up2 and I_Vp2 that resulted from inversion of COMPONENT 2 
C OF PERTURBATION 3-D Ertel PV distribution into RAWINS file. 

C NOTE: Vis5d output is the new sea level pressure (hPa)

C See below ROMU: For choices on variables etc...
C ==============================================================================

      INTEGER FUNCTION USERFUNC( OUTGRID, OUTNL, OUTLOWLEV,
     *                           INGRID, NR, NC, NL, LOWLEV, MAXNL,
     *                           NVARS, NAMES,
     *                           DATE, TIME,
     *                           PROJECTION, PROJ_ARGS,
     *                           VERTICAL, VERT_ARGS )
C ARGUMENTS:
      implicit none
      INTEGER NVARS
      INTEGER NR, NC, NL(NVARS), LOWLEV(NVARS), MAXNL
      REAL OUTGRID(NR,NC,MAXNL)
      INTEGER OUTNL, OUTLOWLEV
      REAL INGRID(NR,NC,MAXNL,NVARS)
      CHARACTER*8 NAMES(NVARS)
      INTEGER DATE, TIME
      INTEGER PROJECTION
      REAL PROJ_ARGS(*)
      INTEGER VERTICAL
      REAL VERT_ARGS(*)

C LOCAL VARIABLES:
      integer var, ir, ic, il, if, nf
      integer ii_up2, ii_vp2, ii_hp2, ii_tp2
      integer imonth,julian,days_month(12)
      integer icycle,iRAWINS_inputunit,iRAWINS_outputunit,
     $        iRAWINS_outputexists,IYYDDD,
     $        date_found,YYYY,MM,DD,HH,MI,SS         
      character*24 date_v5d
      real grav, cp, rd, pi, conv, ds

      parameter (grav=9.81,cp=1004.5,rd=2.*cp/7.,
     $           pi=3.141592,conv=180./pi)

      real p(22),ud(200,200,22),uc(200,200,22),
     $     vd(200,200,22),vc(200,200,22),
     $     hc(200,200,22),
     $     tc(200,200,22)
      character*100 RAWINS_input,RAWINS_output

C VARIABLES DESCRIBING RAWINS 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(50), start_index(4,50), end_index(4,50)
       Real         xtime(50)
       Character*4  staggering(50), ordering(50)
       Character*24 current_date(50)
       Character*9  name(50)
       Character*25 unit(50)
       Character*46 description(50)                                                                  
       Real t_RAW(200,200,22),u_RAW(200,200,22),v_RAW(200,200,22),
     $      h_RAW(200,200,22),rh_RAW(200,200,22)     
       Real ter_RAW(200,200),land_RAW(200,200),mapc_RAW(200,200),
     $      mapd_RAW(200,200),cor_RAW(200,200),latc_RAW(200,200),
     $      lonc_RAW(200,200),latd_RAW(200,200),lond_RAW(200,200),
     $      snow_RAW(200,200),slpd_RAW(200,200),slpc_RAW(200,200),
     $      sst_RAW(200,200),slpc0_RAW(200,200),
     $      weasd_RAW(200,200),st010_RAW(200,200),st200_RAW(200,200),   
     $      st400_RAW(200,200),sm010_RAW(200,200),sm200_RAW(200,200),
     $      seaice_RAW(200,200),shgt_RAW(200,200),v01_RAW(200,200),
     $      v02_RAW(200,200),v03_RAW(200,200),v04_RAW(200,200),
     $      v05_RAW(200,200),v06_RAW(200,200),v07_RAW(200,200),
     $      v08_RAW(200,200),v09_RAW(200,200),v10_RAW(200,200),
     $      v11_RAW(200,200),v12_RAW(200,200),tgrd_RAW(200,200),
     $      lmask_RAW(200,200),sindx_RAW(200,200)
       Real pressure_RAW(22)
       Real difH,difT,tsl,tavg,
     $      tsfc,ter,rho,rot,fric,spd,dir,dirrot            

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

C ROMU: SPECIFY HERE RAWINS_input AND RAWINS_output FILES

      iRAWINS_inputunit=10
      RAWINS_input=
     $'/disk109/people/grupo2/MM5V3-linux/RAWINS-linux/'//
     $'RAWINS_DOMAIN1'

      iRAWINS_outputunit=20 
      RAWINS_output=
     $'/disk109/people/grupo2/MM5V3-linux/RAWINS-linux/'//
     $'RAWINS_DOMAIN1_noNE'

      print*,'RAWINS_input =',RAWINS_input
      print*,'RAWINS_output =',RAWINS_output
     
C ==============================================================================

C Specify number of levels in OUTGRID
      OUTNL = 1
      OUTLOWLEV = 0

C Find the I_Up2, I_Vp2, I_Hp2, I_Tp2 variables
      ii_up2 = -1
      ii_vp2 = -1
      ii_hp2 = -1
      ii_tp2 = -1
      do var=1,nvars
        if (names(var) .eq. 'I_Up2') ii_up2 = var
        if (names(var) .eq. 'I_Vp2') ii_vp2 = var
        if (names(var) .eq. 'I_Hp2') ii_hp2 = var         
        if (names(var) .eq. 'I_Tp2') ii_tp2 = var                    
      end do

C If I_Up2, I_Vp2, I_Hp2 or I_Tp2 not found, return error 1
      if ( ii_up2.eq.-1 .or. ii_vp2.eq.-1 .or. 
     $     ii_hp2.eq.-1 .or. ii_tp2.eq.-1 ) then
        print *, "Couldn't find I_Up2, I_Vp2, I_Hp2 or I_Tp2 !"
        userfunc = 1
        return
      endif

C Grid-length ds (m)
      if (projection.Eq.2) then
        ds=proj_args(6)*1000.      
      else if (projection.Eq.3) then  
        ds=proj_args(5)*1000.
      else
        print *, "Uncomprensible projection data!"
        userfunc = 1
        return
      endif      
      
C ==============================================================================

C NOTE: In VIS5D    (ir:N-->S, ic:W-->E, il:B-->T)
C NOTE: In RAWINS (ir:S-->N, ic:W-->E, il:B-->T, il=1 is surface)     

C Fill p(Pa); as in RAWINS       
      do il=1,maxnl
        p(il+1)=100.*1012.5*exp(vert_args(il)/-7.2)
      enddo    
       
C Fill uc(m/s), vc(m/s), hc(m), tc(K); as in RAWINS (cross)      
               
      do 10 ir=1,nr
      do 10 ic=1,nc
      do 10 il=1,maxnl
        uc(nr+1-ir,ic,il+1)=ingrid(ir,ic,il,ii_up2)
        vc(nr+1-ir,ic,il+1)=ingrid(ir,ic,il,ii_vp2) 
        hc(nr+1-ir,ic,il+1)=ingrid(ir,ic,il,ii_hp2)                   
        tc(nr+1-ir,ic,il+1)=ingrid(ir,ic,il,ii_tp2) 
10    continue        

C Calculate ud(m/s), vd(m/s); as in RAWINS (dots) 

      DO 20 il=2,maxnl+1      
        do ir=2,nr
        do ic=2,nc
          ud(ir,ic,il)=0.25*(uc(ir,ic,il)+uc(ir-1,ic,il)+
     $                       uc(ir,ic-1,il)+uc(ir-1,ic-1,il))         
          vd(ir,ic,il)=0.25*(vc(ir,ic,il)+vc(ir-1,ic,il)+
     $                       vc(ir,ic-1,il)+vc(ir-1,ic-1,il))      
        enddo
        enddo
        do ir=2,nr
          ud(ir,1,il)=0.5*(uc(ir,1,il)+uc(ir-1,1,il))
          ud(ir,nc+1,il)=0.5*(uc(ir,nc,il)+uc(ir-1,nc,il))          
          vd(ir,1,il)=0.5*(vc(ir,1,il)+vc(ir-1,1,il))
          vd(ir,nc+1,il)=0.5*(vc(ir,nc,il)+vc(ir-1,nc,il))   
        enddo
        do ic=2,nc
          ud(1,ic,il)=0.5*(uc(1,ic,il)+uc(1,ic-1,il))
          ud(nr+1,ic,il)=0.5*(uc(nr,ic,il)+uc(nr,ic-1,il))          
          vd(1,ic,il)=0.5*(vc(1,ic,il)+vc(1,ic-1,il))
          vd(nr+1,ic,il)=0.5*(vc(nr,ic,il)+vc(nr,ic-1,il))    
        enddo
        ud(1,1,il)=uc(1,1,il)
        ud(1,nc+1,il)=uc(1,nc,il) 
        ud(nr+1,nc+1,il)=uc(nr,nc,il)
        ud(nr+1,1,il)=uc(nr,1,il)          
        vd(1,1,il)=vc(1,1,il)
        vd(1,nc+1,il)=vc(1,nc,il) 
        vd(nr+1,nc+1,il)=vc(nr,nc,il)
        vd(nr+1,1,il)=vc(nr,1,il)                                           
20    CONTINUE

C ==============================================================================

C SET/CHECK SOME RAWINS PARAMETERS

      Open(iRAWINS_inputunit,file=RAWINS_input,form='unformatted')
         Rewind(iRAWINS_inputunit)     
         Read(iRAWINS_inputunit,err=9001,end=9002) flag
         Read(iRAWINS_inputunit,err=9001,end=9002) BHI,BHR,BHIC,BHRC
      Close(iRAWINS_inputunit)
      
      If (BHI(1,1).ne.3) Goto 9003
     
      If (BHI(15,1).Ne.0) Then
         If (nr.Ne.BHI(16,1)-1.Or.nc.Ne.BHI(17,1)-1) Goto 9004
      Else
         If (nr.Ne.BHI(5,1)-1.Or.nc.Ne.BHI(6,1)-1) Goto 9004
      Endif 
      If (maxnl.Ne.BHI(12,3)-1) Goto 9004       

      iRAWINS_outputexists=0
      Open(iRAWINS_outputunit,file=RAWINS_output,form='unformatted')
        Rewind(iRAWINS_outputunit)    
        Read(iRAWINS_outputunit,end=9005) flag
        Read(iRAWINS_outputunit,end=9005) BHI,BHR,BHIC,BHRC    
      Close(iRAWINS_outputunit)
      iRAWINS_outputexists=1
9005  Continue
     
C ==============================================================================

C PERFORM CYCLES OVER time series
     
      YYYY=IYYDDD(date)/1000+1900
      If (YYYY.Le.1950) YYYY=YYYY+100
      If (Mod(YYYY-1964,4).Eq.0) days_month(2)=29
      DD=IYYDDD(date)-(IYYDDD(date)/1000)*1000
      MM=1
      Do While (DD.Gt.days_month(MM))
       DD=DD-days_month(MM)
       MM=MM+1
      Enddo
      HH=time/3600
      MI=(time-HH*3600)/60
      SS=(time-HH*3600-MI*60)/60
      date_v5d=Char(YYYY/1000 +48)//
     $         Char(YYYY/100-(YYYY/1000)*10 +48)//
     $         Char(YYYY/10-(YYYY/100)*10 +48 )//
     $         Char(YYYY-(YYYY/10)*10 +48)//
     $         '-'//
     $         Char(MM/10 +48)//
     $         Char(MM-(MM/10)*10 +48)//     
     $         '-'//
     $         Char(DD/10 +48)//
     $         Char(DD-(DD/10)*10 +48)//     
     $         '_'//
     $         Char(HH/10 +48)//
     $         Char(HH-(HH/10)*10 +48)//     
     $         ':'//
     $         Char(MI/10 +48)//
     $         Char(MI-(MI/10)*10 +48)//     
     $         ':'//
     $         Char(SS/10 +48)//
     $         Char(SS-(SS/10)*10 +48)//     
     $         '.0000'                    
      date_found=0 

C Read/Write? HEADER section 

      Open(iRAWINS_inputunit,file=RAWINS_input,form='unformatted')
      Rewind(iRAWINS_inputunit) 
      Read(iRAWINS_inputunit) flag
      Read(iRAWINS_inputunit) BHI,BHR,BHIC,BHRC

      Open(iRAWINS_outputunit,file=RAWINS_output,form='unformatted')
      Rewind(iRAWINS_outputunit) 
      If (iRAWINS_outputexists.Eq.0) Then 
        Write(iRAWINS_outputunit) flag
        Write(iRAWINS_outputunit) BHI,BHR,BHIC,BHRC
      Else
        Read(iRAWINS_outputunit) flag
        Read(iRAWINS_outputunit) BHI,BHR,BHIC,BHRC
      Endif       

5000  CONTINUE
      nf=0

      Read(iRAWINS_inputunit,end=6000) flag        
      DO 1000 While (flag.Eq.1)
         nf=nf+1

C ==============================================================================

C SUBHEADER SECTION for current nf field

         Read(iRAWINS_inputunit) 
     $        ndim(nf),
     $        start_index(1,nf),start_index(2,nf),
     $        start_index(3,nf),start_index(4,nf),
     $        end_index(1,nf),end_index(2,nf),
     $        end_index(3,nf),end_index(4,nf),
     $        xtime(nf),staggering(nf),ordering(nf),
     $        current_date(nf),name(nf),unit(nf),description(nf)       

C ==============================================================================

C FIELD SECTION for current nf field

C ------------------------------------------------------------------------------
        
C READ 3D FIELDS, from DATAGRID_input

C T (K), cross
         If  (name(nf).Eq.'T        ') Read(iRAWINS_inputunit) 
     $   (((t_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)
C U (m/s), dot
         If  (name(nf).Eq.'U        ') Read(iRAWINS_inputunit) 
     $   (((u_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)
C V (m/s), dot
         If  (name(nf).Eq.'V        ') Read(iRAWINS_inputunit) 
     $   (((v_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)
C H (m), cross
         If  (name(nf).Eq.'H        ') Read(iRAWINS_inputunit) 
     $   (((h_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)
C RH (%), cross
         If  (name(nf).Eq.'RH       ') Read(iRAWINS_inputunit) 
     $   (((rh_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)
                    
C ------------------------------------------------------------------------------       

C READ 2D FIELDS, from RAWINS_input 

C TERRAIN (m), cross          
         If  (name(nf).Eq.'TERRAIN  ') Read(iRAWINS_inputunit) 
     $   ((ter_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C LAND USE (category), cross          
         If  (name(nf).Eq.'LAND USE ') Read(iRAWINS_inputunit) 
     $   ((land_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C MAPFACCR (no dim.), cross          
         If  (name(nf).Eq.'MAPFACCR ') Read(iRAWINS_inputunit) 
     $   ((mapc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C MAPFACDT (no dim.), dot          
         If  (name(nf).Eq.'MAPFACDT ') Read(iRAWINS_inputunit) 
     $   ((mapd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C CORIOLIS (1/s), dot          
         If  (name(nf).Eq.'CORIOLIS ') Read(iRAWINS_inputunit) 
     $   ((cor_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C LATITCRS (degree), cross          
         If  (name(nf).Eq.'LATITCRS ') Read(iRAWINS_inputunit) 
     $   ((latc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C LONGICRS (degree), cross          
         If  (name(nf).Eq.'LONGICRS ') Read(iRAWINS_inputunit) 
     $   ((lonc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C LATITDOT (degree), dot       
         If  (name(nf).Eq.'LATITDOT ') Read(iRAWINS_inputunit) 
     $   ((latd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C LONGIDOT (degree), dot          
         If  (name(nf).Eq.'LONGIDOT ') Read(iRAWINS_inputunit) 
     $   ((lond_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SNOWCOVR (0/1 flag), cross          
         If  (name(nf).Eq.'SNOWCOVR ') Read(iRAWINS_inputunit) 
     $   ((snow_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C PSEALVLD (Pa), dot       
         If  (name(nf).Eq.'PSEALVLD ') Read(iRAWINS_inputunit) 
     $   ((slpd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C PSEALVLC (Pa), cross       
         If  (name(nf).Eq.'PSEALVLC ') Read(iRAWINS_inputunit) 
     $   ((slpc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C TSEASFC (K), cross          
         If  (name(nf).Eq.'TSEASFC  ') Read(iRAWINS_inputunit) 
     $   ((sst_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C WEASD (kg m-2), dot       
         If  (name(nf).Eq.'WEASD    ') Read(iRAWINS_inputunit) 
     $   ((weasd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SOILT010 (K), cross          
         If  (name(nf).Eq.'SOILT010 ') Read(iRAWINS_inputunit) 
     $   ((st010_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SOILT200 (K), cross          
         If  (name(nf).Eq.'SOILT200 ') Read(iRAWINS_inputunit) 
     $   ((st200_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SOILT400 (K), cross          
         If  (name(nf).Eq.'SOILT400 ') Read(iRAWINS_inputunit) 
     $   ((st400_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SOILM010 (fraction), cross          
         If  (name(nf).Eq.'SOILM010 ') Read(iRAWINS_inputunit) 
     $   ((sm010_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SOILM200 (fraction), cross          
         If  (name(nf).Eq.'SOILM200 ') Read(iRAWINS_inputunit) 
     $   ((sm200_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SEAICE (0/1 flag), cross          
         If  (name(nf).Eq.'SEAICE   ') Read(iRAWINS_inputunit) 
     $   ((seaice_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C SOILHGT (m), cross          
         If  (name(nf).Eq.'SOILHGT  ') Read(iRAWINS_inputunit) 
     $   ((shgt_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC01 (%), cross          
         If  (name(nf).Eq.'VEGFRC01 ') Read(iRAWINS_inputunit) 
     $   ((v01_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC02 (%), cross          
         If  (name(nf).Eq.'VEGFRC02 ') Read(iRAWINS_inputunit) 
     $   ((v02_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC03 (%), cross          
         If  (name(nf).Eq.'VEGFRC03 ') Read(iRAWINS_inputunit) 
     $   ((v03_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC04 (%), cross          
         If  (name(nf).Eq.'VEGFRC04 ') Read(iRAWINS_inputunit) 
     $   ((v04_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC05 (%), cross          
         If  (name(nf).Eq.'VEGFRC05 ') Read(iRAWINS_inputunit) 
     $   ((v05_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC06 (%), cross          
         If  (name(nf).Eq.'VEGFRC06 ') Read(iRAWINS_inputunit) 
     $   ((v06_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC07 (%), cross          
         If  (name(nf).Eq.'VEGFRC07 ') Read(iRAWINS_inputunit) 
     $   ((v07_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC08 (%), cross          
         If  (name(nf).Eq.'VEGFRC08 ') Read(iRAWINS_inputunit) 
     $   ((v08_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC09 (%), cross          
         If  (name(nf).Eq.'VEGFRC09 ') Read(iRAWINS_inputunit) 
     $   ((v09_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC10 (%), cross          
         If  (name(nf).Eq.'VEGFRC10 ') Read(iRAWINS_inputunit) 
     $   ((v10_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC11 (%), cross          
         If  (name(nf).Eq.'VEGFRC11 ') Read(iRAWINS_inputunit) 
     $   ((v11_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C VEGFRC12 (%), cross          
         If  (name(nf).Eq.'VEGFRC12 ') Read(iRAWINS_inputunit) 
     $   ((v12_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C TEMPGRD (K), cross          
         If  (name(nf).Eq.'TEMPGRD  ') Read(iRAWINS_inputunit) 
     $   ((tgrd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
C LANDMASK (category), cross          
         If  (name(nf).Eq.'LANDMASK ') Read(iRAWINS_inputunit) 
     $   ((lmask_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)                                                       
C SOILINDX (category), cross          
         If  (name(nf).Eq.'SOILINDX ') Read(iRAWINS_inputunit) 
     $   ((sindx_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)   

C ------------------------------------------------------------------------------
      
C READ 1D FIELDS, from RAWINS_input 

C PRESSURE (Pa)          
         If  (name(nf).Eq.'PRESSURE ') Read(iRAWINS_inputunit) 
     $   (pressure_RAW(il),il=1,maxnl+1)   

C ==============================================================================

      Read(iRAWINS_inputunit) flag

1000  CONTINUE

C ==============================================================================

C COMPARISON OF current_date WITH date_v5d (YYYY-MM-DD_HH:MI:SS.0000 dates !!!)

      print*,'TIME PERIOD IN RAWINS_INPUT=',current_date(1)

      If (current_date(1).Eq.date_v5d) Then
         date_found=1
         print*,'...... Modifying fields'
      Endif       

C ------------------------------------------------------------------------------

C SUBSTRACTION OF COMPONENT 2 OF PERTURBATION BALANCED FIELDS, if 
C current_date=date_v5d (il=1 is surface)

      If (current_date(1).Ne.date_v5d) Goto 2500

C Previous slpc0_RAW (Pa), cross
        Do 2005 ir=1,nr
        Do 2005 ic=1,nc                     
          tsl=t_RAW(ir,ic,2)+6.5E-3*h_RAW(ir,ic,2)
          tavg=0.5*(tsl+t_RAW(ir,ic,2))
          slpc0_RAW(ir,ic)=p(2)*exp(grav*h_RAW(ir,ic,2)/(rd*tavg))
2005    Continue   
     
C t_RAW sfc (K), cross
        Do 2010 ir=1,nr
        Do 2010 ic=1,nc
          difH=-hc(ir,ic,2)   
          difT=-tc(ir,ic,2)                                     
          t_RAW(ir,ic,1)=t_RAW(ir,ic,1)+difH*6.5E-3+difT 
2010    Continue    

C t_RAW atm (K), cross
        Do 2020 ir=1,nr
        Do 2020 ic=1,nc
        Do 2020 il=2,maxnl+1              
          t_RAW(ir,ic,il)=t_RAW(ir,ic,il)-tc(ir,ic,il) 
2020    Continue     

C h_RAW atm (m), cross
        Do 2030 ir=1,nr
        Do 2030 ic=1,nc
        Do 2030 il=2,maxnl+1              
          h_RAW(ir,ic,il)=h_RAW(ir,ic,il)-hc(ir,ic,il) 
2030    Continue  

C slpc_RAW (Pa), cross
        Do 2040 ir=1,nr
        Do 2040 ic=1,nc                    
          tsl=t_RAW(ir,ic,2)+6.5E-3*h_RAW(ir,ic,2)
          tavg=0.5*(tsl+t_RAW(ir,ic,2))
          slpc_RAW(ir,ic)=slpc_RAW(ir,ic)+
     $                    p(2)*exp(grav*h_RAW(ir,ic,2)/(rd*tavg))-
     $                    slpc0_RAW(ir,ic)
2040    Continue   

C slpd_RAW (Pa), dots
        Do ir=2,nr
        Do ic=2,nc
          slpd_RAW(ir,ic)=0.25*(slpc_RAW(ir,ic)+slpc_RAW(ir-1,ic)+
     $                          slpc_RAW(ir,ic-1)+slpc_RAW(ir-1,ic-1))           
        Enddo
        Enddo
        Do ir=2,nr
          slpd_RAW(ir,1)=0.5*(slpc_RAW(ir,1)+slpc_RAW(ir-1,1))
          slpd_RAW(ir,nc+1)=0.5*(slpc_RAW(ir,nc)+slpc_RAW(ir-1,nc))          
        Enddo
        Do ic=2,nc
          slpd_RAW(1,ic)=0.5*(slpc_RAW(1,ic)+slpc_RAW(1,ic-1))
          slpd_RAW(nr+1,ic)=0.5*(slpc_RAW(nr,ic)+slpc_RAW(nr,ic-1))          
        Enddo
        slpd_RAW(1,1)=slpc_RAW(1,1)
        slpd_RAW(1,nc+1)=slpc_RAW(1,nc) 
        slpd_RAW(nr+1,nc+1)=slpc_RAW(nr,nc)
        slpd_RAW(nr+1,1)=slpc_RAW(nr,1)           
        
C u_RAW atm (m/s), dot
        Do 2050 ir=1,nr+1
        Do 2050 ic=1,nc+1
        Do 2050 il=2,maxnl+1              
          u_RAW(ir,ic,il)=u_RAW(ir,ic,il)-ud(ir,ic,il) 
2050    Continue  

C v_RAW atm (m/s), dot
        Do 2060 ir=1,nr+1
        Do 2060 ic=1,nc+1
        Do 2060 il=2,maxnl+1              
          v_RAW(ir,ic,il)=v_RAW(ir,ic,il)-vd(ir,ic,il) 
2060    Continue  

C u_RAW and v_RAW sfc (m/s), dot
        Do 2070 ir=2,nr
        Do 2070 ic=2,nc 
          tsfc=0.25*(t_RAW(ir,ic,1)+t_RAW(ir-1,ic,1)+
     $               t_RAW(ir,ic-1,1)+t_RAW(ir-1,ic-1,1))         
          ter=0.25*(ter_RAW(ir,ic)+ter_RAW(ir-1,ic)+
     $              ter_RAW(ir,ic-1)+ter_RAW(ir-1,ic-1))         
          rho=(slpd_RAW(ir,ic)/(rd*tsfc))*
     $        (((44330.-ter)/44330.)**5.256242) 
          u_RAW(ir,ic,1)=
     $                   -0.5*((slpc_RAW(ir,ic)+slpc_RAW(ir,ic-1))-
     $                   (slpc_RAW(ir-1,ic)+slpc_RAW(ir-1,ic-1)))/
     $                   (ds*rho*cor_RAW(ir,ic))   
          v_RAW(ir,ic,1)=
     $                   0.5*((slpc_RAW(ir,ic)+slpc_RAW(ir-1,ic))-
     $                   (slpc_RAW(ir,ic-1)+slpc_RAW(ir-1,ic-1)))/
     $                   (ds*rho*cor_RAW(ir,ic))
          If (ter.Lt.-1..or.ter.Gt.0.) Then
            rot=25./conv
            fric=0.5
          Else
            rot=10./conv
            fric=0.7           
          Endif   
          spd=Sqrt(u_RAW(ir,ic,1)*u_RAW(ir,ic,1)+
     $             v_RAW(ir,ic,1)*v_RAW(ir,ic,1)) 
          If (Abs(v_RAW(ir,ic,1)).Lt.1.E-5) v_RAW(ir,ic,1)=1.E-5 
          dir=Atan(u_RAW(ir,ic,1)/v_RAW(ir,ic,1))
          If (v_RAW(ir,ic,1).Gt.0.) dir=dir+pi
          If (dir.Lt.0.) dir=dir+2.*pi
          If (dir.Gt.2.*pi) dir=dir-2.*pi
          dirrot=dir-rot
          u_RAW(ir,ic,1)=-spd*Sin(dirrot)*fric
          v_RAW(ir,ic,1)=-spd*Cos(dirrot)*fric                                      
2070    Continue         
        Do 2075 ir=2,nr 
          u_RAW(ir,1,1)=u_RAW(ir,2,1) 
          v_RAW(ir,1,1)=v_RAW(ir,2,1) 
          u_RAW(ir,nc+1,1)=u_RAW(ir,nc,1)
          v_RAW(ir,nc+1,1)=v_RAW(ir,nc,1)
2075    Continue
        Do 2080 ic=1,nc+1 
          u_RAW(1,ic,1)=u_RAW(2,ic,1)
          v_RAW(1,ic,1)=v_RAW(2,ic,1)
          u_RAW(nr+1,ic,1)=u_RAW(nr,ic,1) 
          v_RAW(nr+1,ic,1)=v_RAW(nr,ic,1)         
2080    Continue

2500  CONTINUE

C ==============================================================================
  
C WRITE/PASS FIELDS in RAWINS_output 
 
      DO 8000 if=1,nf

        If (iRAWINS_outputexists.Eq.0.Or.date_found.Eq.1) Then
          Write(iRAWINS_outputunit) 1
          Write(iRAWINS_outputunit)
     $        ndim(if),
     $        start_index(1,if),start_index(2,if),
     $        start_index(3,if),start_index(4,if),
     $        end_index(1,if),end_index(2,if),
     $        end_index(3,if),end_index(4,if),
     $        xtime(if),staggering(if),ordering(if),
     $        current_date(if),name(if),unit(if),description(if)   
        Else
          Read(iRAWINS_outputunit) flag
          Read(iRAWINS_outputunit)
     $        ndim(if),
     $        start_index(1,if),start_index(2,if),
     $        start_index(3,if),start_index(4,if),
     $        end_index(1,if),end_index(2,if),
     $        end_index(3,if),end_index(4,if),
     $        xtime(if),staggering(if),ordering(if),
     $        current_date(if),name(if),unit(if),description(if) 
        Endif

C ------------------------------------------------------------------------------  
    
C 3D FIELDS
           
        If (iRAWINS_outputexists.Eq.0.Or.date_found.Eq.1) Then
          If  (name(if).Eq.'T        ') Write(iRAWINS_outputunit) 
     $    (((t_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)       
          If  (name(if).Eq.'U        ') Write(iRAWINS_outputunit) 
     $    (((u_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1) 
          If  (name(if).Eq.'V        ') Write(iRAWINS_outputunit) 
     $    (((v_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)      
          If  (name(if).Eq.'H        ') Write(iRAWINS_outputunit) 
     $    (((h_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1) 
          If  (name(if).Eq.'RH       ') Write(iRAWINS_outputunit) 
     $    (((rh_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)      
        Else
          If  (name(if).Eq.'T        ') Read(iRAWINS_outputunit) 
     $    (((t_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)       
          If  (name(if).Eq.'U        ') Read(iRAWINS_outputunit) 
     $    (((u_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1) 
          If  (name(if).Eq.'V        ') Read(iRAWINS_outputunit) 
     $    (((v_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)      
          If  (name(if).Eq.'H        ') Read(iRAWINS_outputunit) 
     $    (((h_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1) 
          If  (name(if).Eq.'RH       ') Read(iRAWINS_outputunit) 
     $    (((rh_RAW(ir,ic,il),ir=1,nr+1),ic=1,nc+1),il=1,maxnl+1)                    
        Endif
                           
C ------------------------------------------------------------------------------
     
C 2D FIELDS 
             
        If (iRAWINS_outputexists.Eq.0.Or.date_found.Eq.1) Then
          If  (name(if).Eq.'TERRAIN  ') Write(iRAWINS_outputunit) 
     $    ((ter_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LAND USE ') Write(iRAWINS_outputunit) 
     $    ((land_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'MAPFACCR ') Write(iRAWINS_outputunit) 
     $    ((mapc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'MAPFACDT ') Write(iRAWINS_outputunit) 
     $    ((mapd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'CORIOLIS ') Write(iRAWINS_outputunit) 
     $    ((cor_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LATITCRS ') Write(iRAWINS_outputunit) 
     $    ((latc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LONGICRS ') Write(iRAWINS_outputunit) 
     $    ((lonc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LATITDOT ') Write(iRAWINS_outputunit) 
     $    ((latd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LONGIDOT ') Write(iRAWINS_outputunit) 
     $    ((lond_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'SNOWCOVR ') Write(iRAWINS_outputunit) 
     $    ((snow_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'PSEALVLD ') Write(iRAWINS_outputunit) 
     $    ((slpd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'PSEALVLC ') Write(iRAWINS_outputunit) 
     $    ((slpc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'TSEASFC  ') Write(iRAWINS_outputunit) 
     $    ((sst_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'WEASD    ') Write(iRAWINS_outputunit) 
     $    ((weasd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILT010 ') Write(iRAWINS_outputunit) 
     $    ((st010_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILT200 ') Write(iRAWINS_outputunit) 
     $    ((st200_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILT400 ') Write(iRAWINS_outputunit) 
     $    ((st400_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILM010 ') Write(iRAWINS_outputunit) 
     $    ((sm010_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILM200 ') Write(iRAWINS_outputunit) 
     $    ((sm200_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SEAICE   ') Write(iRAWINS_outputunit) 
     $    ((seaice_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILHGT  ') Write(iRAWINS_outputunit) 
     $    ((shgt_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)         
          If  (name(if).Eq.'VEGFRC01 ') Write(iRAWINS_outputunit) 
     $    ((v01_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)         
          If  (name(if).Eq.'VEGFRC02 ') Write(iRAWINS_outputunit) 
     $    ((v02_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC03 ') Write(iRAWINS_outputunit) 
     $    ((v03_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC04 ') Write(iRAWINS_outputunit) 
     $    ((v04_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC05 ') Write(iRAWINS_outputunit) 
     $    ((v05_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC06 ') Write(iRAWINS_outputunit) 
     $    ((v06_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC07 ') Write(iRAWINS_outputunit) 
     $    ((v07_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC08 ') Write(iRAWINS_outputunit) 
     $    ((v08_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC09 ') Write(iRAWINS_outputunit) 
     $    ((v09_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC10 ') Write(iRAWINS_outputunit) 
     $    ((v10_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC11 ') Write(iRAWINS_outputunit) 
     $    ((v11_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC12 ') Write(iRAWINS_outputunit) 
     $    ((v12_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'TEMPGRD  ') Write(iRAWINS_outputunit) 
     $    ((tgrd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'LANDMASK ') Write(iRAWINS_outputunit) 
     $    ((lmask_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'SOILINDX ') Write(iRAWINS_outputunit) 
     $    ((sindx_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
        Else
          If  (name(if).Eq.'TERRAIN  ') Read(iRAWINS_outputunit) 
     $    ((ter_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LAND USE ') Read(iRAWINS_outputunit) 
     $    ((land_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'MAPFACCR ') Read(iRAWINS_outputunit) 
     $    ((mapc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'MAPFACDT ') Read(iRAWINS_outputunit) 
     $    ((mapd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'CORIOLIS ') Read(iRAWINS_outputunit) 
     $    ((cor_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LATITCRS ') Read(iRAWINS_outputunit) 
     $    ((latc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LONGICRS ') Read(iRAWINS_outputunit) 
     $    ((lonc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LATITDOT ') Read(iRAWINS_outputunit) 
     $    ((latd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'LONGIDOT ') Read(iRAWINS_outputunit) 
     $    ((lond_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'SNOWCOVR ') Read(iRAWINS_outputunit) 
     $    ((snow_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'PSEALVLD ') Read(iRAWINS_outputunit) 
     $    ((slpd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'PSEALVLC ') Read(iRAWINS_outputunit) 
     $    ((slpc_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'TSEASFC  ') Read(iRAWINS_outputunit) 
     $    ((sst_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)
          If  (name(if).Eq.'WEASD    ') Read(iRAWINS_outputunit) 
     $    ((weasd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILT010 ') Read(iRAWINS_outputunit) 
     $    ((st010_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILT200 ') Read(iRAWINS_outputunit) 
     $    ((st200_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILT400 ') Read(iRAWINS_outputunit) 
     $    ((st400_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILM010 ') Read(iRAWINS_outputunit) 
     $    ((sm010_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILM200 ') Read(iRAWINS_outputunit) 
     $    ((sm200_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SEAICE   ') Read(iRAWINS_outputunit) 
     $    ((seaice_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)            
          If  (name(if).Eq.'SOILHGT  ') Read(iRAWINS_outputunit) 
     $    ((shgt_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)         
          If  (name(if).Eq.'VEGFRC01 ') Read(iRAWINS_outputunit) 
     $    ((v01_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)         
          If  (name(if).Eq.'VEGFRC02 ') Read(iRAWINS_outputunit) 
     $    ((v02_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC03 ') Read(iRAWINS_outputunit) 
     $    ((v03_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC04 ') Read(iRAWINS_outputunit) 
     $    ((v04_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC05 ') Read(iRAWINS_outputunit) 
     $    ((v05_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC06 ') Read(iRAWINS_outputunit) 
     $    ((v06_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC07 ') Read(iRAWINS_outputunit) 
     $    ((v07_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC08 ') Read(iRAWINS_outputunit) 
     $    ((v08_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC09 ') Read(iRAWINS_outputunit) 
     $    ((v09_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC10 ') Read(iRAWINS_outputunit) 
     $    ((v10_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC11 ') Read(iRAWINS_outputunit) 
     $    ((v11_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'VEGFRC12 ') Read(iRAWINS_outputunit) 
     $    ((v12_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'TEMPGRD  ') Read(iRAWINS_outputunit) 
     $    ((tgrd_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'LANDMASK ') Read(iRAWINS_outputunit) 
     $    ((lmask_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)      
          If  (name(if).Eq.'SOILINDX ') Read(iRAWINS_outputunit) 
     $    ((sindx_RAW(ir,ic),ir=1,nr+1),ic=1,nc+1)                
        Endif          

C ------------------------------------------------------------------------------
 
C 1D FIELDS

        If (iRAWINS_outputexists.Eq.0.Or.date_found.Eq.1) Then
          If  (name(if).Eq.'PRESSURE ') Write(iRAWINS_outputunit) 
     $    (pressure_RAW(il),il=1,maxnl+1)
        Else
          If  (name(if).Eq.'PRESSURE ') Read(iRAWINS_outputunit) 
     $    (pressure_RAW(il),il=1,maxnl+1)            
        Endif          


C ------------------------------------------------------------------------------

8000  CONTINUE

      If (iRAWINS_outputexists.Eq.0.Or.date_found.Eq.1) Then
        Write(iRAWINS_outputunit) 2 
      Else
        Read(iRAWINS_outputunit) flag
      Endif

C ==============================================================================

C OUTPUT OF THE NEW SEA LEVEL PRESSURE slpc_RAW (hPa), if current_date=date_v5d

        If (current_date(1).Ne.date_v5d) Goto 3500

          print*,
     $    'Output of new sea level pressure slpc_RAW (hPa)'
          DO 9999 ir=1,nr
          DO 9999 ic=1,nc
            outgrid(nr+1-ir,ic,OUTNL)=slpc_RAW(ir,ic)/100.
9999      CONTINUE    

3500    Continue 
                              
C ==============================================================================

C END OF THIS CYCLE (Go to next time period) 

      GOTO 5000

6000  CONTINUE 

      If (date_found.Eq.0) print*,
     $'DATE_V5D NOT FOUND IN RAWINS_INPUT FILE !!!'

      Close(iRAWINS_inputunit)
      Close(iRAWINS_outputunit)       

C ==============================================================================

C RETURN AND ERROR MESSAGES IF ANY

      USERFUNC = 0
      RETURN     

C Error: Error in reading RAWINS_input header .........
9001  Continue   
      Print *,'  '                                                             
      Print *,'ERROR IN READING RAWINS_INPUT HEADER'                                          
      Return                                                              
                                                              
C Error: EOF detected in RAWINS_input .......... 
9002  Continue     
      Print *,'  '                                                              
      Print *,'EOF DETECTED INITIALLY IN RAWINS_INPUT'                             
      Return     

C Error: RAWINS_input is not a RAWINS file .........
9003  Continue   
      Print *,'  '                                                             
      Print *,'ERROR: RAWINS_INPUT NOT A RAWINS FILE'                      
      Return     

C Error: Grid dimensions in Vis5d not matching RAWINS_input ........
9004  Continue   
      Print *,'  '                                                             
      Print *,'ERROR: GRID DIMENSIONS IN VIS5D NOT MATCHING ',
     $        'THOSE OF RAWINS_INPUT'                              
      Return         

      END

C ==============================================================================
