#include "BULK_FORCE_OPTIONS.h"

CBOP
C     !ROUTINE: BULKF_FORMULA_LAY
C     !INTERFACE:
      SUBROUTINE BULKF_FORMULA_LAY(
     I                           uw, vw, ws, Ta, Qa, tsfCel,
     O                           flwupa, flha, fsha, df0dT,
     O                           ust, vst, evp, ssq, dEvdT,
     I                           iceornot, i,j,bi,bj,myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE BULKF_FORMULA_LAY
C     | o Calculate bulk formula fluxes over open ocean or seaice
C     |   Large and Yeager, 2004, NCAR/TN-460+STR.
C     *==========================================================*
C     \ev
C
C === Turbulent Fluxes :
C  * use the approach "B": shift coeff to height & stability of the
C    atmosphere state (instead of "C": shift temp & humid to the height
C    of wind, then shift the coeff to this height & stability of the atmos).
C  * similar to EXF (except over sea-ice) ; default parameter values
C    taken from Large & Yeager.
C  * assume that Qair & Tair inputs are from the same height (zq=zt)
C  * formulae in short:
C     wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
C     Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
C     Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
C                      = -Evap * Lvap
C   with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
C        del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ;
C        Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number
C              respectively [no-units], function of height & stability

C     !USES:
       IMPLICIT NONE
C     === Global variables ===
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "BULKF_PARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     input:
      _RL uw                 ! zonal wind speed (at grid center) [m/s]
      _RL vw                 ! meridional wind speed (at grid center) [m/s]
      _RL ws                 ! wind speed        [m/s]   at height zwd
      _RL Ta                 ! air temperature   [K]     at height zth
      _RL Qa                 ! specific humidity [kg/kg] at heigth zth
      _RL tsfCel             ! sea-ice or sea surface temperature [oC]
      INTEGER iceornot       ! 0=open water, 1=sea-ice, 2=sea-ice with snow
      INTEGER i,j, bi,bj     !current grid point indices
      INTEGER myThid         ! my Thread Id number
C     output:
      _RL flwupa             ! upward long wave radiation (>0 upward) [W/m2]
      _RL flha               ! latent heat flux         (>0 downward) [W/m2]
      _RL fsha               ! sensible heat flux       (>0 downward) [W/m2]
      _RL df0dT              ! derivative of heat flux with respect to Tsf [W/m2/K]
      _RL ust                ! zonal wind stress (at grid center)     [N/m2]
      _RL vst                ! meridional wind stress (at grid center)[N/m2]
      _RL evp                ! evaporation rate (over open water) [kg/m2/s]
      _RL ssq                ! surface specific humidity          [kg/kg]
      _RL dEvdT              ! derivative of evap. with respect to Tsf [kg/m2/s/K]
CEOP

#ifdef ALLOW_BULK_FORCE

C     == Local variables ==
      _RL dflhdT             ! derivative of latent heat with respect to T
      _RL dfshdT             ! derivative of sensible heat with respect to T
      _RL dflwupdT           ! derivative of long wave with respect to T

      _RL Tsf                ! surface temperature [K]
      _RL Ts2                ! surface temperature square [K^2]
c     _RL ht                 ! height for air temperature [m]
c     _RL hq                 ! height for humidity [m]
c     _RL hu                 ! height for wind speed [m]
c     _RL zref               ! reference height [m]
      _RL wsm                ! limited wind speed [m/s] (> umin)
      _RL usn                ! neutral, zref (=10m) wind speed [m/s]
      _RL usm                ! usn but limited [m/s] (> umin)
c     _RL umin               ! minimum wind speed used for drag-coeff [m/s]
      _RL lath               ! latent heat of vaporization or sublimation [J/kg]
      _RL t0                 ! virtual temperature [K]
      _RL delth              ! potential temperature diff [K]
      _RL delq               ! specific humidity difference [kg/kg]
      _RL ustar              ! friction velocity [m/s]
      _RL tstar              ! temperature scale [K]
      _RL qstar              ! humidity scale  [kg/kg]
      _RL rd                 ! = sqrt(Cd)          [-]
      _RL re                 ! = Ce / sqrt(Cd)     [-]
      _RL rh                 ! = Ch / sqrt(Cd)     [-]
      _RL rdn, ren, rhn      ! neutral, zref (=10m) values of rd, re, rh
      _RL stable             ! = 1 if stable ; = 0 if unstable
      _RL huol               ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
      _RL htol               ! stability parameter at zth [-]
      _RL x                  ! stability function  [-]
      _RL xsq                ! = x^2               [-]
      _RL psimh              ! momentum stability function
      _RL psixh              ! latent & sensib. stability function
      _RL czol               ! = zref*Karman_cst*gravity
      _RL zwln               ! = log(zwd/zref)
      _RL ztln               ! = log(zth/zref)
c     _RL cdalton            ! coeff to evaluate Dalton Number
c     _RL mixratio
c     _RL ea
c     _RL psim_fac
      _RL tau                ! surface stress  coef = rhoA * Ws * Cd
      _RL csha               ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir
      _RL clha               ! latent heat flx coef = rhoA * Ws * Ce * Lvap
c     _RL zice
c     _RL ssq0, ssq1, ssq2   ! constant used in saturated specific humidity
c     _RL p0                 ! reference sea-level atmospheric pressure [mb]
      _RL qs1w, qs2w         !   above freezing saturated specific humidity
      _RL qs1i, qs2i         !   below freezing saturated specific humidity
      _RL tmpBlk
      _RL half, one, two
      INTEGER iter

C     == external Functions

C--   Constant
      DATA   half,      one,      two
     &     / 0.5 _d 0 , 1. _d 0 , 2. _d 0 /
c     DATA   ssq0,           ssq1,           ssq2
c    &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
c     DATA   p0 / 1013. _d 0 /
      DATA   qs1w,           qs2w
     &     /   640.38 _d 3 , 5107.0 _d -0 /
      DATA   qs1i,           qs2i
     &     / 11637.80 _d 3 , 5897.8 _d -0 /

C-- Set surface parameters :
c             zice = 0.0005 _d 0
              zwln = LOG(zwd/zref)
              ztln = LOG(zth/zref)
              czol = zref*xkar*gravity

C-   Surface Temp.
              Tsf = tsfCel+Tf0kel
              Ts2 = Tsf*Tsf
C-   Wind speed
              IF (ws.EQ.0. _d 0) THEN
                ws = SQRT(uw*uw + vw*vw)
              ENDIF
              wsm = MAX(ws,umin)

C--- Compute turbulent surface fluxes
C-   Pot. Temp and saturated specific humidity
              t0     = Ta*(one + humid_fac*Qa)
              IF ( iceornot.EQ.0 ) THEN
                lath=Lvap
                ssq = saltQsFac*qs1w*EXP( -qs2w/Tsf ) / rhoA
                dEvdT = qs2w
              ELSE
                lath = Lvap+Lfresh
                ssq =           qs1i*EXP( -qs2i/Tsf ) / rhoA
                dEvdT = qs2i
              ENDIF
c             ssq = ssq0*EXP( lath*(ssq1-ssq2/Tsf) ) / p0
c             dEvdT = lath*ssq2

              delth  = Ta + gamma_blk*zth - Tsf
              delq   = Qa - ssq

C--  initial guess for exchange coefficients:
C    take U_N = del.U ; stability from del.Theta ;
              stable = half + SIGN(half, delth)
              tmpBlk = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
              rdn = SQRT(tmpBlk)
              rhn = stable*cStantonS + (one-stable)*cStantonU
              ren = cDalton
c             rdn=xkar/(LOG(zref/zice))
c             rhn=rdn
c             ren=rdn
C--  calculate turbulent scales
              ustar=rdn*wsm
              tstar=rhn*delth
              qstar=ren*delq

C--- iterate with psi-functions to find transfer coefficients
              DO iter=1,blk_nIter

                 huol   = ( tstar/t0
     &                     +qstar/(Qa + one/humid_fac)
     &                    )*czol/(ustar*ustar)
                 huol   = SIGN( MIN(abs(huol),10. _d 0), huol)
                 stable = half + SIGN(half, huol)
                 xsq    = SQRT( ABS(one - huol*16. _d 0) )
                 x      = SQRT(xsq)
                 psimh = -5. _d 0*huol*stable
     &             + (one-stable)*
     &                    ( LOG( (one + two*x + xsq)*(one+xsq)*.125 )
     &                     -two*ATAN(x) + half*pi )
                 htol   = huol*zth/zwd
                 xsq    = SQRT( ABS(one - htol*16. _d 0) )
                 psixh  = -5. _d 0*htol*stable
     &             + (one-stable)*( two*LOG(half*(one+xsq)) )

C-   Shift wind speed using old coefficient
                 usn = ws/(one + rdn/xkar*(zwln-psimh) )
                 usm = MAX(usn, umin)

C-   Update the 10m, neutral stability transfer coefficients
                 tmpBlk = cdrag_1/usm + cdrag_2 + cdrag_3*usm
                 rdn = SQRT(tmpBlk)
                 rhn = stable*cStantonS + (one-stable)*cStantonU
                 ren = cDalton

C-   Shift all coefficients to the measurement height and stability.
                 rd = rdn/(1. _d 0 + rdn*(zwln-psimh)/xkar)
                 rh = rhn/(1. _d 0 + rhn*(ztln-psixh)/xkar)
                 re = ren/(1. _d 0 + ren*(ztln-psixh)/xkar)

C--  Update ustar, tstar, qstar using updated, shifted coefficients.
                 ustar = rd*wsm
                 qstar = re*delq
                 tstar = rh*delth

              ENDDO

C-   Coeff:
              tau   = rhoA*rd*ws
              csha  = cpAir*tau*rh
              clha  =  lath*tau*re

C-   Turbulent Fluxes
              fsha  = csha*delth
              flha  = clha*delq
              evp   = -flha/lath
              ust   = tau*rd*uw
              vst   = tau*rd*vw

C-   surf.Temp derivative of turbulent Fluxes
              dEvdT  =  tau*re*ssq*dEvdT/Ts2
              dflhdT = -lath*dEvdT
              dfshdT = -csha

C--- Upward long wave radiation
              IF ( iceornot.EQ.0 ) THEN
                flwupa  = ocean_emissivity*stefan*Ts2*Ts2
                dflwupdT= ocean_emissivity*stefan*Ts2*Tsf*4. _d 0
              ELSEIF (iceornot.EQ.2) THEN
                flwupa   = snow_emissivity*stefan*Ts2*Ts2
                dflwupdT = snow_emissivity*stefan*Ts2*Tsf*4. _d 0
              ELSE
                flwupa   =  ice_emissivity*stefan*Ts2*Ts2
                dflwupdT =  ice_emissivity*stefan*Ts2*Tsf*4. _d 0
              ENDIF

C-   Total derivative with respect to surface temperature
              df0dT = -dflwupdT+dfshdT+dflhdT

#endif /*ALLOW_BULK_FORCE*/

      RETURN
      END
