C TFPWPT.f
C
C VIS-5D analysis function to compute Thermal Frontal Parameter, based
C on Wet bulb 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
      integer it,ihum,iwpt,imap,npas
      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)
      real tpi,ews,td,VAPFW,DPTOF,tcp,tconof,wobf,escal
      PARAMETER(rd=287.04,cp=1004.64,npas=0,escal=1.e4)   
C Aunque escal valga 1.e4, las derivadas estan hechas en Km, por tanto
C nos reescalan 1.e6
	print*,'Reescalado en 1.e10,unidades de  1.e-4 °C/m^2'
c	print*,'Durran filter passed ',npas,' times'
C     Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
      
C   Find the T, RH, WPT and  MAPFACCR variables
      it = -1
      ihum=-1
      iwpt=-1
      imap=-1
      do var=1,nvars
         if (NAMES(var).eq."T") it = var
         if (NAMES(var).eq."RH") ihum = var
         if (NAMES(var).eq."WPT") iwpt = var
         if (names(var) .eq. 'MAPFACCR') imap = var           
      end do
C     If Tem not found, return error 1
      if (it .eq. -1.or.ihum.eq.-1 .or. imap.eq.-1) then
         print *, "Couldn't find T, RH 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
      
            
	if (iwpt.eq.-1) then
	print*,'Calculating WPT internally'
C     Calculo de la temperatura potencial del termometro humedo
      do il=1,maxnl
         do ic=1,NC
            do ir=1,NR
        p=1012.5*exp(vert_args(il)/-7.2)
C Para aprovechar las variables del programa TFPPT se ha renombrado 
C tp->tpi i tp ha pasado a ser la temperatura del termometro humedo
        tpi=(ingrid(ir,ic,il,it)+273.16)*(1000./p)
     &  **(rd/cp)-273.16
	EWS=VAPFW(ingrid(ir,ic,il,it))
	TD=DPTOF(EWS*ingrid(ir,ic,il,ihum)/100.)
	TCP=TCONOF(ingrid(ir,ic,il,it),TD)
	tp(ir,ic,il)=tpi-WOBF(tpi)+WOBF(TCP)
            end do
         end do
      end do
	else
	print*,'Reading calculated WPT'
      do il=1,maxnl
         do ic=1,NC
            do ir=1,NR
	tp(ir,ic,il)=ingrid(ir,ic,il,iwpt)
            end do
         end do
      end do
	endif
C Fin del calculo de la temperatura potencial de termometro humedo

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))*escal
	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 --------------------------------------------------------------------
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
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 ===========================================================================
C	----------------------------------------------------------
C	ALGORITMO PARA CALCULAR LA FUNCION DE WOBUS
C	ENTRADA TEMPERATURA EN C

	FUNCTION WOBF(T)
	
	X=T-20.
	IF(X) 10,10,20
10	CONTINUE

	POL=1.+X*(-8.8416605E-3+X*(1.4714143E-4+X*(-9.671989E-7
     #  +X*(-3.2607217E-8+X*(-3.8598073E-10)))))
	POL=POL*POL
	WOBF=15.130/(POL*POL)
	RETURN

20	POL=1.+X*(3.6182989E-3+X*(-1.3603273E-5+X*(4.9618922E-7
     #  +X*(-6.1059365E-9+X*(3.9401551E-11+X*(-1.2588129E-13
     #  +X*(1.668828E-16)))))))
	POL=POL*POL
	WOBF=29.930/(POL*POL)+.96*X-14.8
	RETURN
	END
c ===========================================================================



c ===========================================================================
C	ALGORITMO PARA CALCULAR LA PRESION DE VAPOR SATURANTE
C	RESPECTO AL AGUA EN hPa.
c	ENTRADA TEMPERATURA EN C.

	FUNCTION VAPFW(T)

	X=T
	IF(100.-X)3,4,4

3	X=X-273.16

C	AJUSTE DE LA CURVA ENTRE -50 < T < 100
4	POL=.99999683E-0+X*(-.90826951E-2+X*(.78736169E-4+
     1	X*(-.61117958E-6+X*(.43884187E-8+X*(-.29883885E-10+
     2	X*(.21874425E-12+X*(-.17892321E-14+X*(.11112018E-16+
     3	X*(-.30994571E-19)))))))))
	POL=POL*POL
	POL=POL*POL
	VAPFW=6.1078/(POL*POL)
	RETURN
	END
c ===========================================================================


c ===========================================================================
C	ALGORITMO PARA INVERTIR LA FORMULA DE TETENS
C	CALCULA LA TEMPERATURA EN GRADOS C DADA LA PRESION DE VAPOR
C	SATURANTE EN hPa

	FUNCTION DPTOF(EW)

	TOL=.0001
	IF (EW-.21382876E-9)20,20,30

20	DPTOF=-10000.
	RETURN

30	IF (1013.-EW)20,100,100

100	X=ALOG(EW/6.1078)
	BOT=17.269388-X
	DPTOF=(273.3*X)/BOT
	BOT=BOT*EW
200	EDP=VAPFW(DPTOF)
	DTDE=(DPTOF+237.7)/BOT
	DELT=DTDE*(EW-EDP)
	DPTOF=DPTOF+DELT

	IF(ABS(DELT)-TOL)300,300,200

300	DPTOF=DPTOF-TOL
	RETURN
	END
c ===========================================================================


c ===========================================================================
C	-------------------------------------------------------------
C	ALGORITMO PARA CALCULAR LA TEMPERATURA DE CONDENSACION
C	PO ELEVACION. ENTRADA TEMPERATURA Y PUNTO DE ROCIO AMBOS
C	EN C o K. LA TEMPERATURA DE CONDENSACION SE DEVUELVE EN C.


	FUNCTION TCONOF(TEMP,DEWPT)
	S=TEMP-DEWPT
	T=TEMP
	IF(100.0-TEMP)4,5,5
4	T=TEMP-273.16
C	AJUSTE DE UNA CURVA
5	DLT=S*(1.2185+.001278*T+S*(-.00219+11.73E-6*S-5.2E-6*T))
	TCONOF=T-DLT
	RETURN
	END
c ===========================================================================



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