#include "BULK_FORCE_OPTIONS.h"

CBOP
C     !ROUTINE: BULKF_FORMULA_LANL
C     !INTERFACE:
      SUBROUTINE BULKF_FORMULA_LANL(
     I                           uw, vw, us, Ta, Qa, nc, tsf_in,
     O                           flwupa, flha, fsha, df0dT,
     O                           ust, vst, evp, ssq, dEvdT,
     I                           iceornot, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE BULKF_FORMULA_LANL
c     | o Calculate bulk formula fluxes over open ocean or seaice
C     *==========================================================*
C     \ev
C swd -- bulkf formula used in bulkf and ice pkgs
C        taken from exf package
C
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 = transfer coefficient for momentum, sensible
C                    & latent heat flux [no units]
C     *==========================================================*

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 us                 ! wind speed        [m/s]   at height hu
      _RL Ta                 ! air temperature   [K]     at height ht
      _RL Qa                 ! specific humidity [kg/kg] at heigth ht
      _RL nc                 ! fraction cloud cover
      _RL tsf_in             ! sea-ice or sea surface temperature [oC]
      INTEGER iceornot       ! 0=open water, 1=sea-ice, 2=sea-ice with snow
      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 ht                 ! height for air temperature [m]
c     _RL hq                 ! height for humidity [m]
      _RL hu                 ! height for wind speed [m]
c     _RL zref               ! reference height [m]
      _RL usm                ! wind speed limited [m/s]
c     _RL umin               ! minimum wind speed used for drag-coeff [m/s]
      _RL lath               ! latent heat of vaporization or sublimation
      _RL t0                 ! virtual temperature [K]
      _RL deltap             ! 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      ! initial (neutral) values of rd, re, rh
      _RL stable             ! = 1 if stable ; = 0 if unstable
      _RL huol               ! stability parameter [-]
      _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 aln                ! = log(ht/zref)
c     _RL cdalton            ! coeff to evaluate Dalton Number
c     _RL mixratio
c     _RL ea
c     _RL psim_fac
      _RL tau                ! surface stress  coef = ?
      _RL csha               ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir
      _RL clha               ! latent heat flx coef = rhoA * Ws * Ce * Lvap
      _RL zice
      _RL ssq0, ssq1, ssq2   ! constant used in surface specific humidity
      _RL p0                 ! reference sea-level atmospheric pressure [mb]
      _RL bulkf_Cdn          ! drag coefficient
      INTEGER niter_bulk, iter

C     == external Functions
c     _RL       exf_BulkCdn
c     external  exf_BulkCdn
c     _RL       exf_BulkqSat
c     external  exf_BulkqSat
c     _RL       exf_BulkRhn
c     external  exf_BulkRhn

      DATA   ssq0,           ssq1,           ssq2
     &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
      DATA   p0 / 1013. _d 0 /

C--- Compute turbulent surface fluxes
              ht =  2. _d 0
c             hq =  2. _d 0
              hu = 10. _d 0
c             zref = 10. _d 0
              zice = 0.0005 _d 0
              aln = log(ht/zref)
              niter_bulk = 5
c             cdalton = 0.0346000 _d 0
              czol = zref*xkar*gravity
c             psim_fac=5. _d 0
c             umin=1. _d 0

              lath=Lvap
              if (iceornot.gt.0) lath=Lvap+Lfresh
              Tsf=Tsf_in+Tf0kel
C-   Wind speed
              if (us.eq.0. _d 0) then
                us = sqrt(uw*uw + vw*vw)
              endif
              usm = max(us,umin)
c
              t0     = Ta*(1. _d 0 + humid_fac*Qa)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
c             ssq = 3.797915 _d 0*exp(
c    &                lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
c    &                               )/p0
              ssq = ssq0*exp( lath*(ssq1-ssq2/Tsf) ) / p0

              deltap = ta  - tsf + gamma_blk*ht
              delq   = Qa - ssq

C--  initialize estimate exchange coefficients
              rdn=xkar/(log(zref/zice))
              rhn=rdn
              ren=rdn
C--  calculate turbulent scales
              ustar=rdn*usm
              tstar=rhn*deltap
              qstar=ren*delq

C--  iteration with psi-functions to find transfer coefficients
              do iter=1,niter_bulk
                 huol   = czol/ustar**2 *(tstar/t0 +
     &                    qstar/(1. _d 0/humid_fac+Qa))
                 huol   = sign( min(abs(huol),10. _d 0), huol)
                 stable = 5. _d -1 + sign(5. _d -1 , huol)
                 xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
                 x      = sqrt(xsq)
                 psimh = -5. _d 0*huol*stable + (1. _d 0-stable)*
     &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +
     &                     2. _d 0*log(5. _d -1*(1. _d 0+xsq)) -
     &                     2. _d 0*atan(x) + pi*.5 _d 0)
                 psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*
     &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))

C--  Update the transfer coefficients
                 rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
                 rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
                 re = rh
C--  Update ustar, tstar, qstar using updated, shifted coefficients.
                 ustar = rd*usm
                 qstar = re*delq
                 tstar = rh*deltap
              enddo

              tau   = rhoa*ustar**2
              tau   = tau*us/usm
              csha  = rhoa*cpair*us*rh*rd
              clha  = rhoa*lath*us*re*rd

              fsha  = csha*deltap
              flha  = clha*delq
              evp   = -flha/lath

C--  Upward long wave radiation
cQQ           mixratio=Qa/(1-Qa)
cQQ           ea=p0*mixratio/(0.62197+mixratio)
cQQ           flwupa=-0.985*stefan*tsf**4
cQQ  &                  *(0.39-0.05*sqrt(ea))
cQQ  &                  *(1-0.6*nc**2)
              if (iceornot.eq.0) then
               flwupa=ocean_emissivity*stefan*tsf**4
               dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3
              elseif (iceornot.eq.2) then
                flwupa=snow_emissivity*stefan*tsf**4
                dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3
              else
                flwupa=ice_emissivity*stefan*tsf**4
                dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3
              endif
cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
c             dflhdT = -clha*ssq*Lath*2.166847 _d -3/(Tsf**2)
              dEvdT  =  clha*ssq*ssq2/(Tsf*Tsf)
              dflhdT = -lath*dEvdT
              dfshdT = -csha
cQQ           dflwupdT= 4.*0.985*stefan*tsf**3
cQQ  &                  *(0.39-0.05*sqrt(ea))
cQQ  &                  *(1-0.6*nc**2)
c total derivative with respect to surface temperature
              df0dT=-dflwupdT+dfshdT+dflhdT

C--  Wind stress at center points
C-   in-lining of function: exf_BulkCdn(umps) = cdrag_1/umps + cdrag_2 + cdrag_3*umps
              bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
              ust = rhoa*bulkf_Cdn*us*uw
              vst = rhoa*bulkf_Cdn*us*vw
#endif /*ALLOW_BULK_FORCE*/

      RETURN
      END
