C ErPV_Is.f
C
C VIS-5D function to calculate INTERPOLATED 3-D Ertel PV field to ISENTROPIC
C surface (ErPV_Is)

C See below ROMU: For personal 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 VARS:
      integer var, ir, ic, il
      integer it, ierpv 
      real po, cp, rd, kap

      parameter (po=100000.,cp=1004.5,rd=2.*cp/7.,kap=rd/cp)
      
      real p(22),exnf(22),theta(200,200,22),q(200,200,22),
     $    thetaval,thetadif1,thetadif2
      
C ==============================================================================
      
C Specify number of levels in OUTGRID
      OUTNL = 1
      OUTLOWLEV = 0

C ROMU: Set which Isentropic surface (in K) here !
      thetaval=330.
      
C Find the T, ErPV variables (ROMU: Set appropiate ErPV field here !)
      it = -1
      ierpv = -1
      do var=1,nvars
        if (names(var) .eq. 'T') it = var
        if (names(var) .eq. 'ErPV') ierpv = var        
      end do

C If T or ErPV not found, return error 1
      if ( it.eq.-1 .or. ierpv.eq.-1 ) then
        print *, 
     $  "Couldn't find T or ErPV vars!"
        userfunc = 1
        return
      endif

C ==============================================================================
C NOTE: In VIS5D (ir:N-->S, ic:W-->E, il:B-->T)
      
C Fill p(Pa), exnf(J kg-1 K-1), theta(K) and q(PVU)
       
      do il=1,maxnl
        p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
        exnf(il)=cp*((p(il)/po)**kap)  
      enddo
               
      do 10 ir=1,nr
      do 10 ic=1,nc
      do 10 il=1,maxnl           
        theta(ir,ic,il)=(ingrid(ir,ic,il,it)+273.16)*cp/exnf(il)
        q(ir,ic,il)=ingrid(ir,ic,il,ierpv)        
10    continue   

C ==============================================================================
C OUTPUT OF INTERPOLATED 3-D ERTEL PV FIELD TO ISENTROPIC SURFACE ErPV_Is (PVU)

      print*,
     $'Output of INTERPOLATED 3-D Ertel PV field to ISENTROPIC surface',
     $' ErPV_Is (PVU)'
      DO 9999 ir=1,nr
      DO 9999 ic=1,nc
        thetadif1=theta(ir,ic,1)-thetaval
        Do il=2,maxnl
          thetadif2=theta(ir,ic,il)-thetaval
          If (thetadif1*thetadif2.Le.0.) Then
            outgrid(ir,ic,OUTNL)=(q(ir,ic,il-1)*thetadif2-
     $                            q(ir,ic,il)*thetadif1)/ 
     $                           (thetadif2-thetadif1)
            Goto 9999
          Endif
          thetadif1=thetadif2
        Enddo
9999  CONTINUE     

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

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