SUBROUTINE SCREENL(USTAR,ZO,TG,QG,VMOD,
     *           HSENS,HLAT,RHO,PSS,theta1,QS0,THETA0,
     & THETAS,QSTAR, PSC,RIB, DELU,DELTHE,DELQS0, UP,TSP,QSP,RL,
     *                 ILONGL,TSAIR,QSAIR,U10)

c $Id: screenl.f,v 2.4 1998/05/21 16:57:53 mwehner Exp $
C
C     subroutine to calculate the "screen height air" temperature and
C      "screen height" mixing ratio 
C      and "anemometer level" wind using the 
C      predictor corrector method of Hess and McAvaney
C
C     Input variables
C       USTAR - friction velocity (m/sec) 
C       ZO    - rougness length (m)
C       TG    - temperature (K) at first model level above surface
C       QG    - mixing ratio (kgm/kgm) at first model level above surface
C       VMOD  - magnitude of wind (m/sec) at first model level above surface
C       HSENS - sensible heat flux at first model level above surface (w/m^2)
C             - positive out of surface
C       HLAT  - latent heat flux at first model level above surface (w/m^2)
C             - positive out of surface
C       RHO   - air density at first model level above surface
C       PSS   - surface pressure (Pa)
C       QS0   - mixing ratio at surface (kgm/kgm)
C      THETA0 - potential temperature at surface (K)
C       ILONGL - number of points to be calculated
C     THETA1 - potential temperature at lowest model level
C
C
      implicit none
      real*8 cp, rgas, roncp, hl, grav, vkar, zsc, zan, eps
      PARAMETER (CP=1.00464d3,RGAS=287.04,RONCP=RGAS/CP,HL=2.5104d6)
      PARAMETER (GRAV=9.80616,VKAR=0.4,ZSC=2.0,ZAN=10.,EPS=1.0E-20)
C
C      CP
C      RGAS - gas constant
C      HL - latent heat of vaporisation
C      GRAV - gravity aceeleration
C      VKAR - vom karman constant
C      ZSC  - screen height (m)
C      ZAN - anemometer height (m)
C      EPS - small "test" number
C
C
C     input
      integer ilongl
      real*8 USTAR(ILONGL),ZO(ILONGL),TG(ILONGL),QG(ILONGL),
     *          HSENS(ILONGL),VMOD(ILONGL),
     *          HLAT(ILONGL),RHO(ILONGL),PSS(ILONGL)
      real*8 THETA0(ILONGL),QS0(ILONGL)
      real*8 THETA1(ILONGL)
C
C    Output variables
C       TSAIR - screen height air temperature at height ZSC 
C          (ZSC is set by parameter statement below)
C       QSAIR - screen height mixing ratio at height ZSC
C       U10   - magnitude of wind at anemometer level ZAN (set by parameter below)
C
      real*8 TSAIR(ILONGL),QSAIR(ILONGL),U10(ILONGL)
C
C     Local variables
C     THETAS - potential temperature scaling
C     QSTAR - mixing ratio scaling
C     PSC    -pressure at "screen" height (current)
C     RIB    - Richardson number at "screen " height (diagnostic)
C     UP     - current wind at "screen height"
C     TSP    -current temperature at screen height
C     QSP    -current mixing ratio at screen height
C     DELU   - correction to current guess for wind
C     DELTHE - correction to current guess of potential temperature
C     DELQS0 - correction to current mixing ratio
C     RL    - Monin Obukhov length - for diagnostic purposes
C
      real*8 THETAS(ILONGL),QSTAR(ILONGL)
      real*8 PSC(ILONGL),RIB(ILONGL)
      real*8 DELU(ILONGL),DELTHE(ILONGL),DELQS0(ILONGL)
      real*8 UP(ILONGL),TSP(ILONGL),QSP(ILONGL),RL(ILONGL)
      integer ik, iter
      real*8 THETASC


      integer niter
C
C     NITER sets the number of iterations of the corrector
      PARAMETER (NITER=2)
C
C     PREPARE INITIAL VALUES
C
      DO  IK=1,ILONGL
C
c      thetas is as in Eqn 11 of H&M
      THETAS(IK) = -SIGN(max(ABS(HSENS(IK)),EPS),HSENS(IK))
     & / (rho(ik) * cp * ustar(ik))
     
c      qstar is as in Eqn  12 of H&M
      QSTAR(IK) = -HLAT(IK)/(RHO(IK)*HL*USTAR(IK))
C      theta1 is potential temperature at lowest model level above surface
* I changed theta1 to be an input variable so no assumptions about the
* the vertical coordinate are made. (mfw 3/19/98)
c     THETA1(IK) = TG(IK) * (100000.d0/(S*PSS(IK)))**RONCP
C      rl is Monin-Obukhov length
      RL(IK) = (USTAR(IK)**2 * THETA1(IK)) /(VKAR * GRAV * THETAS(IK))
C     ensure RL is not too small
      if (abs(rl(ik)) .lt. eps) rl(ik) = eps
      END DO
C
C      At SCREEN HEIGHT - ZSC
C
C     COMPUTE FIRST PREDICTOR
C
      CALL SCREENP(RL,ZO,THETAS,USTAR,QSTAR,
     *             VMOD,THETA0,THETA1,QS0,QG, ILONGL,
     *             ZSC,
     *                       DELU,DELTHE,DELQS0)
C
      DO IK=1,ILONGL
      UP(IK)=DELU(IK)
      QSP(IK)= QS0(IK) + DELQS0(IK)
      THETASC = THETA0(IK) + DELTHE(IK)
      TSP(IK) = THETASC * (100000.d0/PSS(IK))**(-RONCP)
      END DO
C
C       THE ITERATION
C
      DO ITER=1,NITER
C
C     COMPUTE CORRECTOR
C
      CALL SCREENC(ZO,UP,ZSC,QSP,TSP,PSS,
     *           USTAR,THETAS,QSTAR,THETA0,QS0,ILONGL,
     *           PSC,RIB,DELU,DELTHE,DELQS0)
C
      DO IK=1,ILONGL
      UP(IK)=DELU(IK)
      QSP(IK)= QS0(IK) + DELQS0(IK)
      THETASC = THETA0(IK) + DELTHE(IK)
      TSP(IK) = THETASC * (100000.d0/PSC(IK))**(-RONCP)
      END DO
C
      END DO
C
C     END OF ITERATION
C
      DO IK=1,ILONGL
      TSAIR(IK)=TSP(IK)
      QSAIR(IK)=QSP(IK)
      END DO
C
C      At anemometer level - ZAN
C
      CALL SCREENP(RL,ZO,THETAS,USTAR,QSTAR,
     *             VMOD,THETA0,THETA1,QS0,QG, ILONGL,
     *             ZAN,
     *                       DELU,DELTHE,DELQS0)
C
      DO IK=1,ILONGL
      UP(IK)=DELU(IK)
      QSP(IK)= QS0(IK) + DELQS0(IK)
      THETASC = THETA0(IK) + DELTHE(IK)
      TSP(IK) = THETASC * (100000.d0/PSS(IK))**(-RONCP)
      END DO
C
C     ITERATION
C
      DO ITER=1,NITER
C
      CALL SCREENC(ZO,UP,ZAN,QSP,TSP,PSS,
     *           USTAR,THETAS,QSTAR,THETA0,QS0,ILONGL,
     *           PSC,RIB,DELU,DELTHE,DELQS0)
C
      DO IK=1,ILONGL
      UP(IK)=DELU(IK)
      QSP(IK)= QS0(IK) + DELQS0(IK)
      THETASC = THETA0(IK) + DELTHE(IK)
      TSP(IK) = THETASC * (100000.d0/PSC(IK))**(-RONCP)
      END DO
C
      END DO
C
C     END OF ITERATION
C
      DO IK=1,ILONGL
      U10(IK)=UP(IK)
      END DO
C
      RETURN
      END
      SUBROUTINE SCREENP(RL,ZO,THETAS,USTAR,QSTAR,
     *             VMOD,THETA0,THETA1,QS0,QG, ILONGL,HEIGHT,
     *                    DELU,DELTHE,DELQS0)
C
C     The Predictor equation for screen temps and wind
C      based on dyer-businger
C      As per Hess&Mcavaney, 1995.
C
      implicit none
      real*8 vkar
      PARAMETER (VKAR=0.4)
C
C     INPUT
C     RL   - Monin-Obukhov length
C     ZO   - roughness length
C     THETAS - potential temperature scaling factor
C     USTAR - friction velocity
C     QSTAR - mixing ratio scling factor
C     VMOD - wind magnitude at lowest model level
C     THETA0 - surface potential temperature
C     THETA1  - potential temperature of lowest model level
C     QSO    - mixing ratio at surface
C     QG - mixing ratio at lowest model level
C     ILONGL - number of points
C     HEIGHT - height of level for calculations 
C
      integer ilongl
      real*8 RL(ILONGL),ZO(ILONGL)
      real*8 THETAS(ILONGL),USTAR(ILONGL),QSTAR(ILONGL)
      real*8 VMOD(ILONGL),THETA0(ILONGL),THETA1(ILONGL),QS0(ILONGL),
     *        QG(ILONGL)
      real*8 HEIGHT
C
C     OUTPUT
C     DELU - correction to wind
C     DELTHE - correction to potential temperature
C     DELQS0 - correction to mixing ratio
C
      real*8 DELU(ILONGL),DELTHE(ILONGL),DELQS0(ILONGL)

c     working variables
      real*8 xm, xmh, xm0, xmh0
      integer ik
C
      DO 100 IK=1,ILONGL
C
      IF(RL(IK).GE.0.d00)THEN
C
C     STABLE
C        full formula only when wind exceeds a minimum value (1 m/sec) 
C             and stability not too large
C
        IF(VMOD(IK).GT.1.0.AND.RL(IK).LE.1.0)THEN
C
          DELU(IK) = (USTAR(IK)/VKAR)* 
     *                   ( log((HEIGHT+ZO(IK))/ZO(IK)) +
     *                  min(5.0d0,5.0d0 * (HEIGHT-ZO(IK))/RL(IK) ) )
          DELTHE(IK) = (THETAS(IK)/VKAR) * 
     *                  ( log((HEIGHT+ZO(IK))/ZO(IK)) +
     *                  min(5.0d0, 5.0d0 * (HEIGHT-ZO(IK))/RL(IK)) )
          DELQS0(IK) = (QSTAR(IK)/VKAR) * 
     *                  ( log((HEIGHT+ZO(IK))/ZO(IK)) +
     *                  min(5.0d0, 5.0d0 * (HEIGHT-ZO(IK))/RL(IK)) )
C
        ELSE
C
          DELTHE(IK)=0.1*(THETA1(IK)-THETA0(IK))
          DELQS0(IK)=0.1*(max(QG(IK),0.0)-QS0(IK))
          DELU(IK)=0.1*VMOD(IK)
        ENDIF
      ELSE
C
C     UNSTABLE
C
C      winds must not be too light and  not too unstable
        IF(VMOD(IK).GT.5.0d0.AND.ABS(RL(IK)).LE.50.0)THEN
C
          XM = SQRT(SQRT(1.0d0 - 16.0d0 * HEIGHT/RL(IK)) )    
          XMH = SQRT(1.0d0 - 16.0d0 * HEIGHT/RL(IK))
          XM0 = SQRT(SQRT(1.0d0 - 16.0d0 * ZO(IK)/RL(IK)) )
          XMH0 = SQRT(1.0d0 - 16.0d0 * ZO(IK)/RL(IK)) 
          DELU(IK) = (USTAR(IK)/VKAR) * 
     *       ( log((HEIGHT+ZO(IK))/ZO(IK)) -
     *       2.0d0 *log(0.5d0*(1.+XM)) + 2.0d0 *log(0.5d0*(1.d0+XM0)) -
     *       log(0.5d0*(1.d0+XM**2)) + log(0.5d0*(1.d0+XM0**2)) +
     *      2.0d0 * ATAN(XM) - 2.0d0 * ATAN(XM0)  )
          DELTHE(IK) = (THETAS(IK)/VKAR) * 
     *       ( log((HEIGHT+ZO(IK))/ZO(IK)) -
     *   2.0d0*log(0.5d0*(1.0d0+XMH)) + 2.0d0*log(0.5d0*(1.0d0+XMH0)) )
          DELQS0(IK) = (QSTAR(IK)/VKAR) * 
     *       ( log((HEIGHT+ZO(IK))/ZO(IK)) -
     *   2.0d0*log(0.5d0*(1.0d0+XMH)) + 2.0d0*log(0.5d0*(1.0d0+XMH0)) )
        ELSE
C
C      when winds are light or instability is extreme
          DELTHE(IK)=0.5*(THETA1(IK)-THETA0(IK))
          DELQS0(IK)=0.5*(max(QG(IK),0.0d0)-QS0(IK))
          DELU(IK)=0.5*VMOD(IK)
        ENDIF
C
      ENDIF
C
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE SCREENC(ZO,VMOD,HEIGHT,QG,TG,PSS,
     *           USTAR,THETAS,QSTAR,THETA0,QS0,ILONGL,
     *           PSC,RIBB,DELU,DELTHE,DELQS0)
C
C     THE corrector phase - of computation of screen
C        Temperatures
C      As per Hess&mcavaney, 1995.
C
      implicit none
      real*8 grav, vkar, cp, rgas, roncp
      PARAMETER (GRAV=9.80616,VKAR=0.4)
      PARAMETER (CP=1.00464d3,RGAS=287.04,RONCP=RGAS/CP)
C
C     INPUT
C     ZO - roughness length
C     VMOD - wind magnitude at lowest model level
C     HEIGHT - height at which output required
C     QG  - mixing ratio at lowest model level
C     TG  - temperature at lowest model level
C     PSS - surface pressure
C     USTAR - scaled friction velocity
C     THETAS - scaling for potential temperature
C     QSTAR - scaling for mixing ratio
C     THETA0 - potential temperature at surface
C     QS0    - mixing ratio ai surface
C     ILONGL - number of points
C
C
      integer ilongl
      real*8 ZO(ILONGL),QG(ILONGL),TG(ILONGL)
      real*8 PSS(ILONGL)
      real*8 VMOD(ILONGL)
      real*8 USTAR(ILONGL),THETAS(ILONGL),QSTAR(ILONGL)
      real*8 THETA0(ILONGL),QS0(ILONGL)
      real*8 HEIGHT

C     OUTPUT
C     PSC - pressure at level HEIGHT (diagnostic)
C     RIBB - Richardson muber at HEIGHT (iagnostic)
C     DELU
C     DELTHE
C     DELQS0
      real*8 DELU(ILONGL),DELTHE(ILONGL),DELQS0(ILONGL)
      real*8 PSC(ILONGL),RIBB(ILONGL)

c     working variables
      integer ik
      real*8 qz, qs0z, trm1, trm2, theta1, rib, af, fmu, fmtq
C
      DO 100 IK=1,ILONGL
C
      QZ = max(0.0d0,QG(IK))
      QS0Z = max(0.0d0,QS0(IK))
      TRM1 = 1.0d0 + 0.608 *QS0Z
      TRM2 = 1.0d0 + 0.608 * QZ
      PSC(IK) = EXP( log(PSS(IK)) -
     *      GRAV * HEIGHT / ( RGAS * TG(IK) *
     *     (1+.608*QZ) ) )
      THETA1 = TG(IK) *(100000./PSC(IK))**RONCP
C     Richardson number - Eqn 7 of H&M
      RIB =(GRAV*HEIGHT)/(VMOD(IK)**2) * (1.-THETA0(IK)/THETA1 *
     *               (1.0d0 - TRM1/TRM2) )
C
      AF=(VKAR/log((HEIGHT+ZO(IK))/ZO(IK)))**2
C
C     the following functions are the Louis transfer faunctions as used by BMRC
C     each model needs to use same form as in surface flux calculation
      IF(RIB.GE.0.)THEN
 	   RIB=min(20.,RIB)
 	   FMU = 1.0d0/(1.0d0+10.*RIB*(1.0d0+8.0d0*RIB))
	   FMTQ = FMU
      ELSE
           FMU = 1.0d0 - 10.0d0 * RIB / ( 1. + AF*75.*SQRT(-RIB *
     *                               (HEIGHT+ZO(IK))/ZO(IK)) )
           FMTQ = 1.0d0 - 15.0d0 * RIB / ( 1.d0 + AF*75.d0*SQRT(-RIB *
     *                               (HEIGHT+ZO(IK))/ZO(IK)) )
      ENDIF
C
      RIBB(IK)=RIB
C      Equations 11-13 of H&M
      DELU(IK) = USTAR(IK)/SQRT(AF*FMU)
      DELTHE(IK)=THETAS(IK)*(SQRT(FMU/AF)/FMTQ)
      DELQS0(IK)=QSTAR(IK)*(SQRT(FMU/AF)/FMTQ)
C
  100 CONTINUE
C
      RETURN
      END