#include "CHEAPAML_OPTIONS.h"
#undef ALLOW_THSICE

CBOP
C     !ROUTINE: CHEAPAML_LANL_FLUX
C     !INTERFACE:
      SUBROUTINE CHEAPAML_LANL_FLUX(
     I                    i,j,bi,bj,
     O                    fsha, flha, evp, xolw, ssqt, q100 )

C     !DESCRIPTION:
C     ==================================================================
C     SUBROUTINE cheapaml_LANL_flux
C     ==================================================================
C     o compute surface fluxes using LANL algorithms

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
c#include "FFIELDS.h"
c#ifdef ALLOW_THSICE
c#include "THSICE_VARS.h"
c#endif
#include "CHEAPAML.h"

C     !INPUT PARAMETERS:
      INTEGER i,j,bi,bj
C     !OUTPUT PARAMETERS:
      _RL fsha, flha, evp, xolw, ssqt, q100
CEOP

C       Output:
C       ustress, vstress - wind stress
C       fsha - sensible heat flux
C       flha - latent heat flux
C       xolw - oceanic upwelled long wave radiation
C       ssqt - sat. specific humidity at atm layer top

C       Input:
C       uwind, vwind  - mean wind speed (m/s)
C       Tair  - mean air temperature (K)  at height ht (m)
C       theta(k=1) - sea surface temperature (C)
C       Qair - Specific humidity kg/kg
C       Solar - short wave net solar flux at surface (W/m^2)
C       Tr - relaxation profile for temperature on boundaries (C)
C       qr - relaxation profile for specific humidity (kg/kg)
C       i,j,bi,bj - indices of data

C     !LOCAL VARIABLES:
C     iceornot :: variables to include seaice effect
      INTEGER iceornot
      _RL deltaTm
      _RL uss,usm,uw,vw
      _RL cheapaml_BulkCdn
      _RL to
      _RL t
      _RL t0,QaR
      _RL ssq, q
      _RL deltap, delq, pt, psx100, z100ol
      _RL rdn,ren,rhn,zice,zref
      _RL rd,re,rh,tta,ttas,toa,ttt
      _RL ustar,tstar,qstar,ht,hu,hq
      _RL aln,cdalton,czol,psim_fac
      _RL huol,stable,xsq,x,psimh,psixh
      _RL clha, csha
      INTEGER niter_bulk,iter

C useful values
C hardwire atmospheric relative humidity at 80%
        QaR=0.8 _d 0
C factor to compute rainfall from specific humidity
C inverse of time step
        deltaTm=1. _d 0/deltaT
C reference values to compute turbulent flux
              ht=zt
              hq=zq
              hu=zu
              zref = zt
              zice=.0005 _d 0
              aln = log(ht/zref)
C for iterating on turbulence
              niter_bulk = 5
              cdalton = 0.0346000 _d 0
              czol = zref*xkar*gravity
              psim_fac=5. _d 0

C     determine wind stress
        IF(.NOT.useStressOption)THEN

             if (maskC(i,j,1,bi,bj).ne.0. _d 0) then
#ifdef ALLOW_THSICE
               if (ICEMASK(i,j,bi,bj).gt.0. _d 0) then
                if (snowheight(i,j,bi,bj).gt.3. _d -1) then
                   iceornot=2
                 else
                   iceornot=1
                 endif
               else
                 iceornot=0
               endif
#else
               iceornot=0
#endif
                       uw=uwind(i,j,bi,bj)
                       vw=vwind(i,j,bi,bj)
                       uss=sqrt(uw**2+vw**2)
                       usm=max(uss,1. _d 0)
                  cheapaml_BulkCdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
                       ustress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*uw
                       vstress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*vw
             else
               usm=0. _d 0
               ustress(i,j,bi,bj) = 0. _d 0
               vstress(i,j,bi,bj) = 0. _d 0
                endif
C wind stress computed
                ENDIF
C diabatic and freshwater flux forcing
        to=theta(i,j,1,bi,bj)
        t=Tair(i,j,bi,bj)
        toa=to+Celsius2K
        tta=t+Celsius2K
        ttas=tta+gamma_blk*zref
        ttt=tta-( cheaphgrid(i,j,bi,bj)- zref)*gamma_blk
        pt=p0*(1-gamma_blk*cheaphgrid(i,j,bi,bj)/ttas)
     &     **(gravity/gamma_blk/gasR)

C specific humidities
              ssq= ssq0*exp( lath*(ssq1-ssq2/toa) ) / p0
              ssqt = ssq0*exp( lath*(ssq1-ssq2/ttt) ) / pt
C     saturation no more at the top:
              ssqt = 0.7 _d 0*ssq

            if (useFreshWaterFlux) then
            q=qair(i,j,bi,bj)
            else
            q=QaR * ssq
            endif

C adjust temperature from reference height to formula height
            deltap = t  - to + gamma_blk*(zref-ht)
            delq   = q - ssq
            ttas   = tta+gamma_blk*(zref-ht)
            t0     = ttas*(1. _d 0 + humid_fac*q)

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+q))
                 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

                        usm=max(uss,0.5 _d 0)
                csha   = rhoa*cpair*usm*rh*rd
                clha   = rhoa*lath*usm*re*rd

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

C the sensible and latent heat fluxes, fsha and flha,
C are computed so that positive values are downward.
C the convention for cheapaml is upward fluxes are positive,
C so they must be multiplied by -1
        fsha=-fsha
        flha=-flha

C oceanic upwelled long wave
        xolw=stefan*(toa)**4
C compute specific humidity at 100m
                 huol   = czol/ustar**2 *(tstar/t0 +
     &                    qstar/(1. _d 0/humid_fac+q))
                 huol   = sign( min(abs(huol),10. _d 0), huol)
                 stable = 5. _d -1 + sign(5. _d -1 , huol)
                 z100ol   = 100. _d 0 *xkar*gravity/ustar**2 *(tstar/t0
     &                    + qstar/(1. _d 0/humid_fac+q))
                 xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*z100ol)),1. _d 0)
                 x      = sqrt(xsq)
                 psx100  = -5. _d 0*z100ol*stable + (1. _d 0-stable)*
     &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
                 q100=ssq+qstar*(dlog(100. _d 0/zice)-psx100)

      RETURN
      END
