C ErPV_ad.f
C
C VIS-5D function to calculate horizontal ADVECTION of 3-D Ertel PV field 
C (ErPV_ad). 

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 iu, iv, ierpv, imap 
      real ds, ds2
      real u(200,200,22),v(200,200,22),
     $     map(200,200), 
     $     q(200,200,22),
     $     q_ad(200,200,22)

C ==============================================================================
      
C Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C Find the U, V, ErPV variables (ROMU: Set appropiate WIND and PV fields here !)
      iu = -1
      iv = -1 
      ierpv = -1
      do var=1,nvars
        if (names(var) .eq. 'I_U') iu = var  
        if (names(var) .eq. 'I_V') iv = var                 
        if (names(var) .eq. 'ErPV') ierpv = var            
      end do

C Find the MAPFACCR variable 
      imap = -1
      do var=1,nvars
        if (names(var) .eq. 'MAPFACCR') imap = var           
      end do

C If U, V, ErPV or MAPFACCR not found, return error 1
      if ( iu.eq.-1 .or. iv.eq.-1 .or. ierpv.eq.-1 .or.
     $     imap.eq.-1 ) then
        print *, "Couldn't find U, V, ErPV 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 ==============================================================================

C NOTE: In VIS5D (ir:N-->S, ic:W-->E, il:B-->T)

C Fill u(m/s), v(m/s), q(PVU) and map(adim.)
           
      do 10 ir=1,nr
      do 10 ic=1,nc
        map(ir,ic)=ingrid(ir,ic,1,imap) 
      do 10 il=1,maxnl
        u(ir,ic,il)=ingrid(ir,ic,il,iu)
        v(ir,ic,il)=ingrid(ir,ic,il,iv)
        q(ir,ic,il)=ingrid(ir,ic,il,ierpv)           
10    continue        

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

C CALCULATION of horizontal ADVECTION of 3-D Ertel PV field q_ad (PVU s-1) 
               
      DO 20 ir=2,nr-1
      DO 20 ic=2,nc-1
      DO 20 il=1,maxnl
        q_ad(ir,ic,il)=-map(ir,ic)*(
     $                  u(ir,ic,il)*(q(ir,ic+1,il)-q(ir,ic-1,il))+
     $                  v(ir,ic,il)*(q(ir-1,ic,il)-q(ir+1,ic,il))
     $                  )/ds2     
20    CONTINUE

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

      print*,'Calculated q_ad'
 
C ==============================================================================

C OUTPUT of horizontal ADVECTION of 3-D Ertel PV field q_ad (in PVU/12 h)

      print*,
     $'Output of horizontal ADVECTION of 3-D Ertel PV field ',
     $'q_ad (PVU/12 h)'
      DO 9999 ir=1,nr
      DO 9999 ic=1,nc
      DO 9999 il=1,maxnl
        outgrid(ir,ic,il)=12.*3600.*q_ad(ir,ic,il)
9999  CONTINUE    

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

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