C TFPPT.f
C
C VIS-5D analysis function to compute Thermal Frontal Parameter,
C based on Potential Temperature (in 1.E-4 C/m2 ?)

      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,npas
      integer it
      real p,rd,cp,ds
      real tp(100,100,22),gtpx(100,100,22),gtpy(100,100,22)
      real modul(100,100,22),gmx(100,100,22),gmy(100,100,22)
      PARAMETER(rd=287.04,cp=1004.64,npas=0)   

	print*,'Reescalado a unidades de  1.e-4 °C/m^2'
	print*,'Durran filter passed',npas,' times'

C     Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C     Find the T and MAPFACCR variables
      it = -1
      imap = -1
      do var=1,nvars
         if (NAMES(var).eq."T") it = var
         if (names(var) .eq. 'MAPFACCR') imap = var  
      end do
C     If T or MAPFACCR not found, return error 1
      if (it .eq. -1 .or. imap .eq. -1) then
         print *, "Couldn't find T or MAPFACCR variables!"
         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 la temperatura potencial
      do il=1,maxnl
         do ic=1,NC
            do ir=1,NR
        p=1012.5*exp(vert_args(il)/-7.2)
        tp(ir,ic,il)=(ingrid(ir,ic,il,it)+273.16)*(1000./p)
     &  **(rd/cp)-273.16
            end do
         end do
      end do
C Fin del calculo de la temperatura potencial

c Calculo del parametro frontal termico
	call grad(tp,gtpx,gtpy,nr,nc,maxnl,ds)
c	print*,gtpx,gtpy
	do 63 il=1,maxnl
	do 63 ir=1,nr
	do 63 ic=1,nc
	gtpx(ir,ic,il)=gtpx(ir,ic,il)*ingrid(ir,ic,1,imap)
	gtpy(ir,ic,il)=gtpy(ir,ic,il)*ingrid(ir,ic,1,imap)	
	modul(ir,ic,il)=SQRT(gtpx(ir,ic,il)**2+gtpy(ir,ic,il)**2)
63	enddo
	call grad(modul,gmx,gmy,nr,nc,maxnl,ds)
	do 638 il=1,maxnl
	do 638 ir=1,nr
	do 638 ic=1,nc
	if (modul(ir,ic,il).eq.0.) then
	outgrid(ir,ic,il)=0.
	else
	gmx(ir,ic,il)=gmx(ir,ic,il)*ingrid(ir,ic,1,imap)
	gmy(ir,ic,il)=gmy(ir,ic,il)*ingrid(ir,ic,1,imap)		
	outgrid(ir,ic,il)=(-(gtpx(ir,ic,il)*gmx(ir,ic,il)+
     &  gtpy(ir,ic,il)*gmy(ir,ic,il))/modul(ir,ic,il))*1.e4
	endif
638	enddo
	do il=1,npas
	call smooth3d(outgrid,nc,nr,maxnl,.0625,.0625)
	enddo
	do il=1,maxnl
        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(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)
       outgrid(nr,nc,il)=2.*outgrid(nr-1,nc-1,il)-outgrid(nr-2,nc-2,il)
      end do




      USERFUNC = 0
      RETURN
      END

C--- Calcula el gradient horitzontal (matox,matoy) d´un camp escalar
c--- tridimensional matin(nr,nc,nl). Omple els contorns amb un ajust 
c--- lineal
	subroutine grad(matin,matox,matoy,nr,nc,nl,ds)
	implicit none
	real matin(100,100,22),matox(100,100,22),matoy(100,100,22)
	real xincx,xincy,ds
	integer nr,nc,nl,il,ic,ir
c Treballem amb escales horitzontals de Km. es analeg a reescalar
c la variable en 1.e6
	do il=1,nl
	do 25 ic=2,nc-1
	xincy=(2.*ds)/1000.
	do 25 ir=2,nr-1
	xincx=(2.*ds)/1000.
	matox(ir,ic,il)=(matin(ir,ic+1,il)-matin(ir,ic-1,il))/xincx
	matoy(ir,ic,il)=(matin(ir-1,ic,il)-matin(ir+1,ic,il))/xincy
c	PRINT*,matox(ir,ic,il),matoy(ir,ic,il),ir,ic,il
25	enddo
        do ir=2,nr-1
	   matox(ir,1,il)=2.*matox(ir,2,il)-matox(ir,3,il)
	   matox(ir,nc,il)=2.*matox(ir,nc-1,il)-matox(ir,nc-2,il)
	   matoy(ir,1,il)=2.*matoy(ir,2,il)-matoy(ir,3,il)
	   matoy(ir,nc,il)=2.*matoy(ir,nc-1,il)-matoy(ir,nc-2,il)
	enddo
	do ic=2,nc-1
	  matox(1,ic,il)=2.*matox(2,ic,il)-matox(3,ic,il)
	  matox(nr,ic,il)=2.*matox(nr-1,ic,il)-matox(nr-2,ic,il)
	  matoy(1,ic,il)=2.*matoy(2,ic,il)-matoy(3,ic,il)
	  matoy(nr,ic,il)=2.*matoy(nr-1,ic,il)-matoy(nr-2,ic,il)
	enddo
	matox(1,1,il)=2.*matox(2,2,il)-matox(3,3,il)
	matox(1,nc,il)=2.*matox(2,nc-1,il)-matox(3,nc-2,il)
       matox(nr,1,il)=2.*matox(nr-1,2,il)-matox(nr-2,3,il)
       matox(nr,nc,il)=2.*matox(nr-1,nc-1,il)-matox(nr-2,nc-2,il)
	matoy(1,1,il)=2.*matoy(2,2,il)-matoy(3,3,il)
	matoy(1,nc,il)=2.*matoy(2,nc-1,il)-matoy(3,nc-2,il)
       matoy(nr,1,il)=2.*matoy(nr-1,2,il)-matoy(nr-2,3,il)
       matoy(nr,nc,il)=2.*matoy(nr-1,nc-1,il)-matoy(nr-2,nc-2,il)
	enddo
	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
	integer j,i,nxs,nys,il,nzs	
	REAL HSX4,HSY4,HSX2,HSY2
	REAL*4 HT(NXS,NYS,NZS), HTS(300,300,20)

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 ===========================================================================
