C QWSLOPE.f
C
C VIS-5D analysis function to compute Upslope vertical moisture flux
C (in 1.E-2 g m-2 s-1),  using (U,V) and QV at 1000 hPa, and TERRAIN
 
	
      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 ir, ic, il, imap
      integer iu,iv,var,i1000,npas,ioro,it,ihum,iqv
      real a1,a2,a4,q(300,300),wslope(300,300)
      real pi,rt,p,xincy,xincx,escal,npasf,ds
      parameter (Pi=3.141596,Rt=6371227.,npas=2,npasf=3,escal=1.e5)
      
C     Specify number of levels in OUTGRID
      OUTNL = 1
      OUTLOWLEV = 0	
	
	print*,'Unidades de  1.e5 Kg/(m2 s)'
	
C Find the U, V, T, RH, TERRAIN and MAPFACCR variables
      iu = -1
      iv = -1
      it = -1
      ihum=-1
      ioro=-1
      iqv= -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."TERRAIN") ioro = var
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'RH') ihum = var
         if (names(var) .eq. 'QV') iqv = var
         if (names(var) .eq. 'MAPFACCR') imap = var          
      end do
      
C  If  not found, return error 1
      if (iu .eq.-1.or.iv.eq.-1.or.ioro.eq.-1.or.it.eq.-1
     $    .or.ihum.eq.-1 .or. imap.eq.-1) then
         print *, "Couldn't find U,V,TERRAIN,T,RH,MAPFACCR s!"
         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   
           
	do il=1,maxnl
	p=1012.5*exp(vert_args(il)/-7.2)
	if (p.gt.999.and.p.lt.1001) i1000=il
	enddo

C Calculo de WSLOPE
        xincy=2.*ds
	do 100 ir=2,nr-1
	xincx=2.*ds
	do 100 ic=2,nc-1
	wslope(ir,ic)=ingrid(ir,ic,i1000,iu)*(ingrid(ir,ic+1,1,ioro)
     &  -ingrid(ir,ic-1,1,ioro))/xincx+ingrid(ir,ic,i1000,iv)*
     &  (ingrid(ir-1,ic,1,ioro)-ingrid(ir+1,ic,1,ioro))/xincy
        wslope(ir,ic)=wslope(ir,ic)*ingrid(ir,ic,1,imap)     
100	enddo

c	do il=1,npas
c	call smooth(wslope,nr,nc,0.0625,0.0625)
c	enddo

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

c Fin calculo de WSLOPE


C Calculo de QV

c Comprovam si tenim QV calculat
	if (iqv.eq.-1) then
	print*,'Calculating QV1000 internally'
c Calcula los niveles de presion (hPA)
	p=1012.5*exp(vert_args(i1000)/-7.2)	
	do 142 ir=1,nr
	do 142 ic=1,nc
C ---------------------------------------------------------------
C Calcula la humedad especifica en kg/Kg
	a1=ingrid(ir,ic,i1000,it)
	a2=ingrid(ir,ic,i1000,ihum)
	call procio(a1,a2,a4)
	call humesp(p,a4,q(ir,ic))
142	enddo
	else
	print*,'Reading calculated QV'
	do 157 ir=1,nr
	do 157 ic=1,nc
	q(ir,ic)=ingrid(ir,ic,i1000,iqv)/1000.
157	enddo
	print*,'QV llegit',q(5,5)
c Endif del if de iqv
	endif
C Calculo de QVSLOPE
        p=1012.5*exp(vert_args(i1000)/-7.2)*100.
	do 163 ir=1,nr
	do 163 ic=1,nc
	do 163 il=1,outnl       
	  outgrid(ir,ic,il)=(wslope(ir,ic)*q(ir,ic)*p/
     & (287.04*(ingrid(ir,ic,i1000,it)+273.16)))*escal
163	enddo

c	do il=1,npasf
c	call smooth3d(outgrid,nc,nr,maxnl,.0625,.0626)
c	enddo

	do 103 ir=1,nr
	do 103 ic=1,nc
	do 103 il=1,outnl
	if (ingrid(ir,ic,il,ioro).Le.0.) outgrid(ir,ic,il)=0.
103	enddo

      USERFUNC = 0
      RETURN
      END



	SUBROUTINE SMOOTH(HT,NXS,NYS,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
	REAL HSX4,HSY4,HSX2,HSY2
	REAL*4 HT(NXS,NYS), HTS(300,300)

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
c		do 103 il=1,nzs
		DO 204 J=1,NYS

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

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

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


		DO 206 I=1,NXS

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

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

      HTS(I,NYS-1)=HTS(I,NYS-1)-HSY2*(HT(I,NYS)+HT(I,NYS-2)
     *        -2.*HT(I,NYS-1))
  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)=HTS(I,J)
  208		CONTINUE
c103 	enddo
	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
	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 ===========================================================================






c -----------------------------------------------------------------
	subroutine procio(t,u,e)

	es=6.112*EXP((17.67*t)/(t+243.5))
	e=es*u/100.
	return
	end

c ------------------------------------------------------------------
	subroutine humesp (p,e,qq)

	r=.622*e/(p-e)
	qq=r/(1+r)
	return
	end
