c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SFCFLX                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE SFCFLX(nx,ny,nz, u,v,ptprt,qv,rhobar,ptbar,qvbar,
     :                  usflx,vsflx,ptsflx,qvsflx)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the surface momentum , heat and surface moisture fluxes
c     using a simple bulk aerodynamic formulation. These fluxes will
c     be passed into turbulent mixing subroutines.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     4/07/1992.  
c
c     MODIFICATION HISTORY:
c
c     5/06/92 (M. Xue)
c     Added full documentation.
c
c     6/4/92 (M. Xue and H. Jin)
c     Further facelift.
c
c     9/14/1992 (M. Xue)
c     Different drag coefficients defined for momentum, temrperature and
c     moisture. Definition of ptsfc0 and qvsfc0 changed. qvbar added to
c     the argument list.
c
c     9/28/1992 (M. Xue)
c     Sign correction in all the flux formulations.
C
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)
c       ptprt    Perturbation potential temperature at a given time level (K)
c       qv       Water vapor specific humidity at a given time level (kg/kg)
c       rhobar   Base state air density (kg/m**3)
c       ptbar    Base state potential temperature (K)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c     OUTPUT:
c
c       usflx    Surface flux of u-momentum
c       vsflx    Surface flux of v-momentum
c       ptsflx   Surface flux of heat (K*kg/(m*s**2))
c       qvsflx   Surface flux of moisture (K*kg/(m*s**2))
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! The number grid points in 3 directions
c
      real u     (nx,ny,nz)     ! Total u-velocity (m/s)
      real v     (nx,ny,nz)     ! Total v-velocity (m/s)
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)

c
      real usflx (nx,ny)        ! surface flux of u-momentum (kg/(m*s**2))
      real vsflx (nx,ny)        ! surface flux of v-momentum (kg/(m*s**2))
      real ptsflx(nx,ny)        ! surface flux of heat (K*kg/(m*s**2))
      real qvsflx(nx,ny)        ! surface flux of moisture (kg/(m**2*s))
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j
      real    usuf, vsuf, cdrho, rhoavg, ptsfc0, qvsfc0
c
c#######################################################################
c
c     Include files: 
c
c     Include file globcst.inc passes cdm, cdh, cdq
c     and UMOVE and VMOVE into this subroutine.
c
c#######################################################################
c
      include 'globcst.inc'
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      IF ( sfcphy .eq. 0 ) THEN
c
c#######################################################################
c
c     If sfcphy = 0, usflx, vsflx, ptsflx and qvsflx are set to zero and 
c     will not be used in the mixing subroutines.  
c     For safety, we set them to zero.
c
c#######################################################################
c
        DO 150 j=1,ny
        DO 150 i=1,nx
          usflx(i,j) = 0.0
          vsflx(i,j) = 0.0
          ptsflx(i,j) = 0.0
          qvsflx(i,j) = 0.0
150     CONTINUE

      ELSE
c
c#######################################################################
c
c     Compute the surface flux of u-momentum, i.e. the surface
c     drag in x-direction, according to a simple drag law:
c
c     usflx = rhobar*bar( (w'u') ) = rhobar*Cdm*V*u
c
c     where Cdm is the drag coefficient, V the wind speed at the first 
c     scalar level above ground, and u the x-component of wind at the
c     same level. 
c
c     usflx is defined half dz below the u-point, i.e. at the ground level. 
c
c#######################################################################
c

        DO 100 j=2,ny-2
        DO 100 i=2,nx-1

          usuf = u(i,j,2) + umove 

          vsuf = (v(i,j,2)+v(i-1,j,2)+v(i,j+1,2)+v(i-1,j+1,2))*0.25  
     :            +vmove

          usflx(i,j) = cdm*(rhobar(i-1,j,1)+rhobar(i,j,1)
     :                      +rhobar(i-1,j,2)+rhobar(i,j,2))*0.25
     :               * sqrt( usuf*usuf + vsuf*vsuf ) 
     :               * usuf 

100     CONTINUE
c
c#######################################################################
c
c     Compute the surface flux of v-momentum, i.e. the surface
c     drag in y-direction, according to a simple drag law:
c
c     vsflx = - rhobar*bar( (w'v') ) = rhobar*Cdm*V*v
c
c     where Cdm is the drag coefficient, V the wind speed at the first
c     scalar level above ground, and v the y-component of wind at the
c     same level.
c
c     vsflx is defined half dz below the v-point, i.e. at the ground level.
c
c#######################################################################
c
        DO 200 j=2,ny-1
        DO 200 i=2,nx-2

          usuf = (u(i,j,2)+u(i,j-1,2)+u(i+1,j,2)+u(i+1,j-1,2))*0.25 
     :           +umove

          vsuf = v(i,j,2) + vmove 

          vsflx(i,j) = cdm* (rhobar(i,j,1)+rhobar(i,j-1,1)
     :                      +rhobar(i,j,2)+rhobar(i,j-1,2))*0.25 
     :               * sqrt( usuf*usuf + vsuf*vsuf ) 
     :               * vsuf 

200     CONTINUE

c
c#######################################################################
c
c     Compute the surface heat flux according to a simple drag law:
c
c          ptsflx = - rhobar*bar( (w' pt') ) = rhobar*Cdh*V*(pt-ptsfc0)
c
c     where Cdh is the drag coefficient, V the wind speed at the first
c     scalar level above ground and pt the potential temperature at that
c     level.  
c
c     ptsfc0 should be the potential temperature at the ground level,
c     which is taken as the base state temperature at the ground for now.
c
c     ptsflx is defined half dz below the scalar point, therefore is at
c     the ground level. 
c
c     Similarly, the vertical moisture flux at the surface is defined as 
c
c          qvsflx = - rhobar*bar( (w' qv') ) = rhobar*Cdq*V*(qv-qvsfc0)
c
c     where qv is the specific humidity at the first model scalar level.
c     qvsfc0 should be the ground level moisture, and is taken as the 
c     base state moisture at the moment.  
c
c#######################################################################
c
      DO 300 j=2,ny-2
      DO 300 i=2,nx-2

        usuf = (u(i,j,2)+u(i+1,j,2))*0.5 + umove

        vsuf = (v(i,j,2)+v(i,j+1,2))*0.5 + vmove

        rhoavg = (rhobar(i,j,1)+rhobar(i,j,2)) *0.5

        ptsfc0 = (ptbar(i,j,1)+ptbar(i,j,2))*0.5
        qvsfc0 = (qvbar(i,j,1)+qvbar(i,j,2))*0.5

        ptsflx(i,j) = cdh * rhoavg * sqrt( usuf*usuf + vsuf*vsuf )
     :              *( ptprt(i,j,2)+ptbar(i,j,2) - ptsfc0 )

        qvsflx(i,j) = cdq * rhoavg * sqrt( usuf*usuf + vsuf*vsuf )
     :              *( qv(i,j,2) - qvsfc0 )

300   CONTINUE

      ENDIF

      RETURN
      END

