C I_Tp.f
C
C VIS-5D function to calculate PERTURBATION balanced temperature (I_Tp)

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 ii_hp, it, itm  
      real grav, po, cp, rd, kap, ds

      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),
     $     t(200,200,22),
     $     phi(200,200,0:22),
     $     thetaB(200,200)

C ==============================================================================
      
C Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C Find the I_Hp, T and Tm variables
      ii_hp = -1
      it = -1
      itm = -1
      do var=1,nvars
        if (names(var) .eq. 'I_Hp') ii_hp = var      
        if (names(var) .eq. 'T') it = var
        if (names(var) .eq. 'Tm') itm = var        
      end do

C If I_Hp, T or Tm not found, return error 1
      if ( ii_hp.eq.-1 .or. it.eq.-1 .or. itm.eq.-1 ) then
        print *, 
     $  "Couldn't find I_Hp, T or Tm 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 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)
       
      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)

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,it)-
     $                ingrid(ir,ic,1,itm))*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

      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_Hp (phi), used to solve temperature

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

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

C SOLVE TEMPERATURE t (K)

      Do 3500 ir=1,nr
      Do 3500 ic=1,nc
      Do 3500 il=1,maxnl
        t(ir,ic,il)=-Sndim*PI2ndim*
     $              (exnf(il)/cp)*   
     $              (phi(ir,ic,il)-phi(ir,ic,il-1))/
     $              dexnfh(il-1)
      
3500  Continue      

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

C OUTPUT OF PERTURBATION BALANCED TEMPERATURE I_Tp (=t, in C) 

      print*,
     $'Output of PERTURBATION balanced temperature I_Tp (C)'
      DO 9999 ir=1,nr
      DO 9999 ic=1,nc
      DO 9999 il=1,maxnl
        outgrid(ir,ic,il)=t(ir,ic,il)
9999  CONTINUE     

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

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