C I_Vp1.f
C
C VIS-5D function to calculate COMPONENT 1 OF PERTURBATION balanced v-wind
C (I_Vp1) 

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 VARS:
      integer var, ir, ic, il
      integer imap, ii_sp1, itbp1
      real grav, po, cp, rd, kap, ds, ds2

      parameter (grav=9.81,po=100000.,cp=1004.5,rd=2.*cp/7.,kap=rd/cp)

      real Fndim,PI2ndim,PIndim,Lndim,Sndim  
      real p(22),ph(0:22),exnf(22),exnfh(0:22),dexnfh(0:22),
     $     v(200,200,22),
     $     map(200,200),
     $     psi(200,200,0:22),
     $     thetaB(200,200)
     
C ==============================================================================
      
C Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C Find the I_Sp1, TBp1 and MAPFACCR variables
      ii_sp1 = -1
      itbp1 = -1
      imap = -1      
      do var=1,nvars
        if (names(var) .eq. 'I_Sp1') ii_sp1 = var                  
        if (names(var) .eq. 'TBp1') itbp1 = var 
        if (names(var) .eq. 'MAPFACCR') imap = var                  
      end do

C If I_Sp1, TBp1 or MAPFACCR not found, return error 1
      if ( ii_sp1.eq.-1 .or. itbp1.eq.-1 .or. imap.eq.-1 ) then
        print *, "Couldn't find I_Sp1, TBp1 or MAPFACCR vars!"
        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 Definition ds2
      ds2=2.*ds
      
C Constants for nondimensionalization of the system
      Fndim=1.E-4
      PI2ndim=500./Float(maxnl)
      PIndim=PI2ndim**0.5
      Lndim=ds
      Sndim=Fndim*Fndim*Lndim*Lndim/PI2ndim  
      
C ==============================================================================

C NOTE: In VIS5D (ir:N-->S, ic:W-->E, il:B-->T)
      
C Fill p/ph(Pa), exnf/exnfh(J kg-1 K-1) and map(adim.)
       
      do il=1,maxnl
        p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
        exnf(il)=cp*((p(il)/po)**kap)  
      enddo
      do il=1,maxnl-1
        exnfh(il)=0.5*(exnf(il)+exnf(il+1))
        ph(il)=(exnfh(il)*(po**kap)/cp)**(1./kap)
      enddo 
      exnfh(0)=2.*exnf(1)-exnfh(1)
      ph(0)=(exnfh(0)*(po**kap)/cp)**(1./kap)      
      exnfh(maxnl)=2.*exnf(maxnl)-exnfh(maxnl-1)      
      ph(maxnl)=(exnfh(maxnl)*(po**kap)/cp)**(1./kap)

      do 10 ir=1,nr
      do 10 ic=1,nc
        map(ir,ic)=ingrid(ir,ic,1,imap)   
10    continue   

C Definition of Neumann boundary condition thetaB (K)
      do 100 ir=1,nr
      do 100 ic=1,nc
        thetaB(ir,ic)=ingrid(ir,ic,1,itbp1)*
     $                cp/exnf(1)     
100   continue    
               
C Definition dexnfh
      dexnfh(0)=exnfh(1)-exnfh(0)
      do 120 il=1,maxnl-1
       dexnfh(il)=exnfh(il+1)-exnfh(il) 
120   continue      

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

C NON-DIMENSIONALIZATION OF SYSTEM
      ds=ds/Lndim
      ds2=2.*ds
      
      do 910 ir=1,nr
      do 910 ic=1,nc      
        thetaB(ir,ic)=thetaB(ir,ic)/(Sndim*PIndim)          
910   continue     

      exnfh(0)=exnfh(0)/PIndim
      dexnfh(0)=dexnfh(0)/PIndim
      do 930 il=1,maxnl-1
        exnf(il)=exnf(il)/PIndim 
        exnfh(il)=exnfh(il)/PIndim
        dexnfh(il)=dexnfh(il)/PIndim                          
930   continue 
      exnf(maxnl)=exnf(maxnl)/PIndim     
      exnfh(maxnl)=exnfh(maxnl)/PIndim
        
C ==============================================================================

C Retrieve I_Sp1 (psi), used to solve v-wind

      Do 44 ir=1,nr
      Do 44 ic=1,nc
        psi(ir,ic,0)=(grav/Fndim)*ingrid(ir,ic,1,ii_sp1)*
     $               (Fndim/(Sndim*PI2ndim))+
     $               thetaB(ir,ic)*(exnf(1)-exnfh(0))
      Do 44 il=1,maxnl
        psi(ir,ic,il)=psi(ir,ic,il-1)+dexnfh(il-1)*
     $                ((grav/Fndim)*ingrid(ir,ic,il,ii_sp1)*
     $                (Fndim/(Sndim*PI2ndim))-
     $                psi(ir,ic,il-1))/(exnf(il)-exnfh(il-1))
44    Continue
      print*,'Retrieved PSI'

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

C SOLVE V-WIND u (m s-1)

      Do 3500 ir=1,nr
      Do 3500 ic=1,nc
      Do 3500 il=1,maxnl
        v(ir,ic,il)=(Sndim*PI2ndim/Fndim)*map(ir,ic)* 
     $              (
     $              (psi(ir,ic+1,il)-(exnfh(il)-exnf(il))*
     $              (psi(ir,ic+1,il)-psi(ir,ic+1,il-1))/
     $               dexnfh(il-1))
     $              - 
     $              (psi(ir,ic-1,il)-(exnfh(il)-exnf(il))*
     $              (psi(ir,ic-1,il)-psi(ir,ic-1,il-1))/
     $               dexnfh(il-1)) 
     $              ) / (Lndim*ds2)
3500  Continue      

C Extrapolate values at the lateral boundaries
      Do 3550 il=1,maxnl
        Do ir=1,nr
          v(ir,1,il)=2.*v(ir,2,il)-v(ir,3,il)
          v(ir,nc,il)=2.*v(ir,nc-1,il)-v(ir,nc-2,il)          
        Enddo
        Do ic=1,nc
          v(1,ic,il)=2.*v(2,ic,il)-v(3,ic,il)
          v(nr,ic,il)=2.*v(nr-1,ic,il)-v(nr-2,ic,il)          
        Enddo
        v(1,1,il)=-v(2,2,il)+v(2,1,il)+v(1,2,il)
        v(nr,1,il)=-v(nr-1,2,il)+v(nr-1,1,il)+v(nr,2,il)
        v(1,nc,il)=-v(2,nc-1,il)+v(1,nc-1,il)+v(2,nc,il)
        v(nr,nc,il)=-v(nr-1,nc-1,il)+v(nr-1,nc,il)+v(nr,nc-1,il)                              
3550  Continue

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

C OUTPUT OF COMPONENT 1 OF PERTURBATION BALANCED V-WIND I_Vp1 (=v, in m s-1)

      print*,
     $'Output of COMPONENT 1 OF PERTURBATION balanced v-wind ',
     $'I_Vp1 (m s-1)'
      DO 9999 ir=1,nr
      DO 9999 ic=1,nc
      DO 9999 il=1,maxnl
        outgrid(ir,ic,il)=v(ir,ic,il)
9999  CONTINUE     

C ==============================================================================
               
      USERFUNC = 0
      RETURN
      END

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