C EFGDIV.f
C
C VIS-5D analysis function to compute Divergence Term of 
C Escalar Frontogenesis (in 1.E-9 K / s m)


      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,imap
      integer iu, iv, iw, it, ifilt
      real Pi,Rt,Per,thetax,thetay,thetap,thetagrad,ux,uy,vx,vy,
     $     omex,omey,exnf,tem,ds
      real p(22),rocp,rg,g,xincx,xincy
      parameter (Pi=3.141596,Rt=6371227.,Per=86400.,rocp=2./7.,
     $           rg=287.04,g=9.807)
      real ome(200,200,50), theta(200,200,50), u(200,200,50),
     $     v(200,200,50)  
      real EFGDIV,Div
      

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


C     Find the U, V, W, T, and MAPFACCR variables
      iu = -1
      iv = -1
      iw = -1
      it = -1
      imap = -1
      do var=1,nvars
         if (names(var) .eq. 'U') iu = var
         if (names(var) .eq. 'V') iv = var
         if (names(var) .eq. 'W') iw = var         
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'MAPFACCR') imap = var              
      end do

C     If U, V, W, T, or MAPFACCR not found, return error 1
      if (iu .eq. -1 .or. iv .eq. -1.or. iw .eq. -1 .or. it .eq. -1
     $    .or. imap .eq. -1) then
         print *, "Couldn't find U, V, W, T, or MAPFACCR vars!"
         userfunc = 1
         return
      endif

      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 Calculo de P (Pa) y Conversion de W a ome (barica, en Pa/s), y de T a theta (K)
C y relleno de u y v 
      do il=1,maxnl
       p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
       exnf=(100000./p(il))**rocp   
        do ic=1,nc
         do ir=1,nr
          tem=ingrid(ir,ic,il,it)+273.16
          ome(ir,ic,il)=-p(il)*g*ingrid(ir,ic,il,iw)/
     $                  (100.*rg*tem)
          theta(ir,ic,il)=tem*exnf
          u(ir,ic,il)=ingrid(ir,ic,il,iu)
          v(ir,ic,il)=ingrid(ir,ic,il,iv)          
         enddo
        enddo
      enddo  

C Smoothing de los campos
      do ifilt=1,300
         Call SMOOTH3D(ome,nc,nr,maxnl,.0625,.0625)
         Call SMOOTH3D(theta,nc,nr,maxnl,.0625,.0625)
         Call SMOOTH3D(u,nc,nr,maxnl,.0625,.0625)
         Call SMOOTH3D(v,nc,nr,maxnl,.0625,.0625)                        
      enddo 

      xincx=2.*ds
      xincy=2.*ds

c Calculo del EFGDIV en 1.E-9 K / m s 
        do 105 il=1,maxnl
        
         do ic=2,nc-1
           do ir=2,nr-1
	     ux=(u(ir,ic+1,il)-u(ir,ic-1,il))/xincx
             ux=ux*ingrid(ir,ic,1,imap)	               
	     uy=(u(ir-1,ic,il)-u(ir+1,ic,il))/xincy
             uy=uy*ingrid(ir,ic,1,imap)	               
	     vx=(v(ir,ic+1,il)-v(ir,ic-1,il))/xincx
             vx=vx*ingrid(ir,ic,1,imap)	               
	     vy=(v(ir-1,ic,il)-v(ir+1,ic,il))/xincy
             vy=vy*ingrid(ir,ic,1,imap)	            
	     omex=(ome(ir,ic+1,il)-ome(ir,ic-1,il))/xincx
             omex=omex*ingrid(ir,ic,1,imap)	               
	     omey=(ome(ir-1,ic,il)-ome(ir+1,ic,il))/xincy
             omey=omey*ingrid(ir,ic,1,imap)	                        
	     thetax=(theta(ir,ic+1,il)-theta(ir,ic-1,il))/xincx
             thetax=thetax*ingrid(ir,ic,1,imap)	               
	     thetay=(theta(ir-1,ic,il)-theta(ir+1,ic,il))/xincy
             thetay=thetay*ingrid(ir,ic,1,imap)	               
             thetagrad=Sqrt(thetax*thetax+thetay*thetay) 
	if (il.eq.1) then
	     thetap=(theta(ir,ic,1)-theta(ir,ic,2))/(p(1)-p(2))
	else if (il.eq.maxnl) then
	     thetap=(theta(ir,ic,maxnl-1)-theta(ir,ic,maxnl))/
     &              (p(maxnl-1)-p(maxnl))
	else
	     thetap=((p(il)-p(il+1))*(theta(ir,ic,il-1)-
     &           theta(ir,ic,il))/(p(il-1)-p(il))+(p(il-1)-p(il))*
     &           (theta(ir,ic,il)-theta(ir,ic,il+1))/(p(il)-p(il+1)))/
     &           (p(il-1)-p(il+1))
	endif
	  Div=ux+vy
	  EFGDIV=-0.5*thetagrad*Div	
	  outgrid(ir,ic,il) = EFGDIV*1.E9
          enddo
        enddo 
     

        do ir=2,nr-1
	 outgrid(ir,1,il)=2.*outgrid(ir,2,il)-outgrid(ir,3,il)
	 outgrid(ir,nc,il)=2.*outgrid(ir,nc-1,il)-outgrid(ir,nc-2,il)
	enddo
	do ic=2,nc-1
	outgrid(1,ic,il)=2.*outgrid(2,ic,il)-outgrid(3,ic,il)
	outgrid(nr,ic,il)=2.*outgrid(nr-1,ic,il)-outgrid(nr-2,ic,il)
	enddo
	outgrid(1,1,il)=2.*outgrid(2,2,il)-outgrid(3,3,il)
	outgrid(nr,nc,il)=2.*outgrid(nr-1,nc-1,il)-outgrid(nr-2,nc-2,il)
	outgrid(1,nc,il)=2.*outgrid(2,nc-1,il)-outgrid(3,nc-2,il)
	outgrid(nr,1,il)=2.*outgrid(nr-1,2,il)-outgrid(nr-2,3,il)

105      end do

      USERFUNC = 0
      RETURN
      END



c ===========================================================================


	SUBROUTINE SMOOTH3d(HT,NXS,NYS,NZS,HSX4,HSY4)

C	==============================================================
C
C	SUBRUTINA SUAVIZACION DE CUARTO ORDEN PARA VALORES DIGITALIZADOS
C	EN UNA MALLA
C	REMITIDA AMABLEMENTE POR DALE R. DURRAN (DEPART. OF METEOROLOGY
C	UNIVERSITY OF UTAH. SALT LAKE CITY. UT 84112)
C	
C	          MODIFICADA POR C. RAMIS (1991) 
C	===============================================================
	implicit none
	REAL HSX4,HSY4,HSX2,HSY2
	REAL*4 HT(200,200,50), HTS(200,200,50)
	integer j,i,nxs,nys,il,nzs
C	
C
C	COEFICIENTES PARA FILTRO DE 4 ORDEN
C
C
	HSX2=-4.*HSX4
	HSY2=-4.*HSY4
C
C	=========================================================
C	HSX4=HSY4=.06 (.0625 MEJOR, ELIMINA 2DX EN UNA PASADA)
C	=========================================================
C

C
C	SUAVIZACION DE LOS CAMPOS
C
		do 103 il=1,nzs
		DO 204 J=1,NYS

		DO 203 I=3,NXS-2
	HTS(I,J,il)=HT(I,J,il)-HSX4*(HT(I+2,J,il)+HT(I-2,J,il)
     *        -4.*(HT(I+1,J,il)+HT(I-1,J,il))+6.*HT(I,J,il))
  203		CONTINUE

      HTS(2,J,il)=HT(2,J,il)-HSX2*(HT(3,J,il)+HT(1,J,il)-2.*HT(2,J,il))

      HTS(NXS-1,J,il)=HT(NXS-1,J,il)-HSX2*(HT(NXS,J,il)+HT(NXS-2,J,il)
     *         -2.*HT(NXS-1,J,il))
  204		CONTINUE


		DO 206 I=1,NXS

		DO 205 J=3,NYS-2
	HTS(I,J,il)=HTS(I,J,il)-HSY4*(HT(I,J+2,il)+HT(I,J-2,il)
     *        -4.*(HT(I,J+1,il)+HT(I,J-1,il))+6.*HT(I,J,il))
  205		CONTINUE

      HTS(I,2,il)=HTS(I,2,il)-HSY2*(HT(I,3,il)+HT(I,1,il)-2.*HT(I,2,il))

      HTS(I,NYS-1,il)=HTS(I,NYS-1,il)-HSY2*(HT(I,NYS,il)+HT(I,NYS-2,il)
     *        -2.*HT(I,NYS-1,il))
  206		CONTINUE

C
C	REASIGNACION DE VALORES A LA MATRIZ ORIGINAL
C

		DO 208 I=2,NXS-1
		DO 208 J=2,NYS-1
	HT(I,J,il)=HTS(I,J,il)
  208		CONTINUE
103 	enddo
	RETURN
	END

c ===========================================================================
