C WSLOPE.f
C
C VIS-5D analysis function to compute Upslope vertical wind
C (in cm/s), using (U,V) at 1000 hPa, and TERRAIN

C See example.f for more information about writing functions.

	
      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
      integer iu,iv,var,i1000,npas,ioro,imap
      real pi,rt,p,xincy,xincx,escal,ds
      parameter (Pi=3.141596,Rt=6371227.,npas=2,escal=1.e2)
      
C     Specify number of levels in OUTGRID
      OUTNL = 1
      OUTLOWLEV = 0
      
	print*,'Reescalado en 1.e2'
C  Find the U, V, TERRAIN and MAPFACCR variables
      iu = -1
      iv = -1
      ioro =-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. 'MAPFACCR') imap = var 
      end do
C  If U, V, TERRAIN or MAPFACCR not found, return error 1
      if (iu .eq. -1.or.iv .eq. -1
     $    .or. ioro.eq.-1 .or. imap.eq.-1) then
         print *, "Couldn't find U, V, TERRAIN, 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   
                  
       do il=1,maxnl
	p=1012.5*exp(vert_args(il)/-7.2)
	if (p.gt.999.and.p.lt.1001) i1000=il
       enddo
	
        xincy=2.*ds
	do 100 ir=2,nr-1
	xincx=2.*ds
	do 100 ic=2,nc-1
	do 100 il=1,OUTNL
        outgrid(ir,ic,il)=ingrid(ir,ic,i1000,iu)*
     &  (ingrid(ir,ic+1,il,ioro)-ingrid(ir,ic-1,il,ioro))/xincx+
     &  ingrid(ir,ic,i1000,iv)*
     &  (ingrid(ir-1,ic,il,ioro)-ingrid(ir+1,ic,il,ioro))/xincy
        outgrid(ir,ic,il)=outgrid(ir,ic,il)*ingrid(ir,ic,1,imap)*
     &  escal
100	enddo

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

	do il=1,OUTNL
        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)
	enddo
	
	do 107 ir=1,nr
	do 107 ic=1,nc
	do 107 il=1,OUTNL
	if (ingrid(ir,ic,il,ioro).Le.0.) outgrid(ir,ic,il)=0.
107	enddo

      USERFUNC = 0
      RETURN

      END




	SUBROUTINE SMOOTH(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 ===========================================================================
