#include "BULK_FORCE_OPTIONS.h"

CBOP
C     !ROUTINE: BULKF_FORMULA_AIM
C     !INTERFACE:
      SUBROUTINE BULKF_FORMULA_AIM(
     I                   Tsurf, SLRD,
     I                   T1, T0, Q0, Vsurf,
     O                   SHF, EVAP, SLRU,
     O                   dEvp, sFlx,
     I                   iceornot, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R BULKF_FORMULA_AIM
C     | o compute surface flux over ocean and sea-ice,
C     |   using AIM surface flux formulation
C     *==========================================================*
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     Resolution parameters

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "BULKF_PARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C--   Input:
C    FMASK  :: fractional land-sea mask        (2-dim)
C    Tsurf  :: surface temperature        (2-dim)
C    SSR    :: sfc sw radiation (net flux)     (2-dim)
C    SLRD   :: sfc lw radiation (downward flux)(2-dim)
C    T1     :: near-surface air temperature (from Pot.temp)
C    T0     :: near-surface air temperature    (2-dim)
C    Q0     :: near-surface sp. humidity [g/kg](2-dim)
C    Vsurf  :: surface wind speed        [m/s] (2-dim,input)
C--   Output:
C    SHF    :: sensible heat flux              (2-dim)
C    EVAP   :: evaporation [g/(m^2 s)]         (2-dim)
C    SLRU   :: sfc lw radiation (upward flux)  (2-dim)
C    Shf0   :: sensible heat flux over freezing surf.
C    dShf   :: sensible heat flux derivative relative to surf. temp
C    Evp0   :: evaporation computed over freezing surface (Ts=0.oC)
C    dEvp   :: evaporation derivative relative to surf. temp
C    Slr0   :: upward long wave radiation over freezing surf.
C    dSlr   :: upward long wave rad. derivative relative to surf. temp
C    sFlx   :: net heat flux (+=down) except SW, function of surf. temp Ts:
C              0: Flux(Ts=0.oC) ; 1: Flux(Ts^n) ; 2: d.Flux/d.Ts(Ts^n)
C    TSFC   :: surface temperature (clim.)     (2-dim)
C    TSKIN  :: skin surface temperature        (2-dim)
C--   Input:
C    iceornot :: 0=open water, 1=ice cover
C    myThid :: Thread number for this instance of the routine
C--
      INTEGER NGP
      PARAMETER ( NGP = 1 )
c     _RL  PSA(NGP), FMASK(NGP), EMISloc
      _RL  Tsurf(NGP)
c     _RL  SSR(NGP)
      _RL  SLRD(NGP)
      _RL  T1(NGP), T0(NGP), Q0(NGP), Vsurf(NGP)

      _RL  SHF(NGP), EVAP(NGP), SLRU(NGP)
      _RL  dEvp(NGP), sFlx(NGP,0:2)
c     _RL  Shf0(NGP), dShf(NGP), Evp0(NGP), dEvp(NGP)
c     _RL  Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
c     _RL  TSFC(NGP), TSKIN(NGP)

      INTEGER iceornot
      INTEGER myThid
CEOP

#ifdef ALLOW_FORMULA_AIM

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C      FWIND0 = ratio of near-sfc wind to lowest-level wind
C      CHS    = heat exchange coefficient over sea
C      VGUST  = wind speed for sub-grid-scale gusts
C      DTHETA = Potential temp. gradient for stability correction
C      dTstab = potential temp. increment for stability function derivative
C      FSTAB  = Amplitude of stability correction (fraction)
C       P0    = reference pressure                 [Pa=N/m2]
C       GG    = gravity accel.                     [m/s2]
C       RD    = gas constant for dry air           [J/kg/K]
C       CP    = specific heat at constant pressure [J/kg/K]
C       ALHC  = latent heat of condensation        [J/g]
C       ALHF  = latent heat of freezing            [J/g]
C       SBC   = Stefan-Boltzmann constant
C      EMISloc :: longwave surface emissivity
c     _RL FWIND0, CHS, VGUST, DTHETA, dTstab, FSTAB
      _RL P0, ALHC, ALHF, RD, CP, SBC, EMISloc
      EQUIVALENCE ( ocean_emissivity , EMISloc )
      EQUIVALENCE ( Lvap   , ALHC )
      EQUIVALENCE ( Lfresh , ALHF )
      EQUIVALENCE ( Rgas   , RD )
      EQUIVALENCE ( cpair  , CP )
      EQUIVALENCE ( stefan , SBC )

C-- Local variables:
C    PSA     :: norm. surface pressure [p/p0]   (2-dim)
C    DENVV   :: surface flux (sens,lat.) coeff. (=Rho*|V|) [kg/m2/s]
C    CDENVV  :: surf. heat flux (sens.,lat.) coeff including stability effect
C    ALHevp  :: Latent Heat of evaporation
      _RL PSA(NGP)
      _RL DENVV(NGP), PRD
      _RL  Shf0(NGP), dShf(NGP), Evp0(NGP)
      _RL  Slr0(NGP), dSlr(NGP)
      _RL  TSFC(NGP), TSKIN(NGP)
      _RL CDENVV(NGP), RDTH, FSSICE
      _RL ALHevp, Fstb0, dTstb, dFstb
      _RL QSAT0(NGP,2)
      _RL QDUMMY(1), RDUMMY(1), TS2
      INTEGER J
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      PSA(1) = 1. _d 0
      P0 = 1. _d +5
C---

      ALHevp = ALHC
C     Evap of snow/ice: account for Latent Heat of freezing :
c     IF ( aim_energPrecip .OR. useThSIce ) ALHevp = ALHC + ALHF
      IF ( iceornot.GE.1 ) ALHevp = ALHC + ALHF

C     1.4 Density * wind speed (including gustiness factor)

      PRD = P0/RD
c     VG2 = VGUST*VGUST
c     factWind2 = FWIND0*FWIND0

      DO J=1,NGP
c       SPEED0(J)=SQRT(factWind2*Vsurf2(J)+VG2)
c       DENVV(J)=(PRD*PSA(J)/T0(J))*SPEED0(J)
C-- assuming input file "WspeedFile" contains the time-average "SPEED0"
C     from AIM output (aimPhytave: fields # 15 ; aimDiag: WINDS ) :
        DENVV(J)=(PRD*PSA(J)/T0(J))*Vsurf(J)
      ENDDO

C     1.5 Define effective skin temperature to compensate for
C         non-linearity of heat/moisture fluxes during the daily cycle

      DO J=1,NGP
        TSKIN(J) = Tsurf(J) + celsius2K
        TSFC(J)=273.16 _d 0
      ENDDO

C--   2. Computation of fluxes over land and sea

C     2.1 Stability correction

      RDTH = FSTAB/DTHETA

      DO J=1,NGP
        FSSICE=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH
        CDENVV(J)=CHS*DENVV(J)*FSSICE
      ENDDO

      IF ( dTstab.GT.0. _d 0 ) THEN
C-    account for stability function derivative relative to Tsurf:
C note: to avoid discontinuity in the derivative (because of min,max), compute
C   the derivative using the discrete form: F(Ts+dTstab)-F(Ts-dTstab)/2.dTstab
       DO J=1,NGP
        Fstb0 = 1.+MIN(DTHETA,MAX(-DTHETA,TSFC(J) -T1(J)))*RDTH
        Shf0(J) = CHS*DENVV(J)*Fstb0
        dTstb = ( DTHETA+dTstab-ABS(TSKIN(J)-T1(J)) )/dTstab
        dFstb = RDTH*MIN(1. _d 0, MAX(0. _d 0, dTstb*0.5 _d 0))
        dShf(J) = CHS*DENVV(J)*dFstb
       ENDDO
      ENDIF

C     2.2 Evaporation

      CALL BULKF_SH2RH_AIM( 2, NGP, TSKIN, PSA, 1. _d 0,
     &                      QDUMMY, dEvp,  QSAT0(1,1), myThid )
      CALL BULKF_SH2RH_AIM( 0, NGP, TSFC,  PSA, 1. _d 0,
     &                      QDUMMY, RDUMMY,QSAT0(1,2), myThid )

      IF ( dTstab.GT.0. _d 0 ) THEN
C-    account for stability function derivative relative to Tsurf:
       DO J=1,NGP
        EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
        Evp0(J) =   Shf0(J)*(QSAT0(J,2)-Q0(J))
        dEvp(J) = CDENVV(J)*dEvp(J)
     &            + dShf(J)*(QSAT0(J,1)-Q0(J))
       ENDDO
      ELSE
       DO J=1,NGP
        EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
        Evp0(J) = CDENVV(J)*(QSAT0(J,2)-Q0(J))
        dEvp(J) = CDENVV(J)*dEvp(J)
       ENDDO
      ENDIF

C     2.3 Sensible heat flux

      IF ( dTstab.GT.0. _d 0 ) THEN
C-    account for stability function derivative relative to Tsurf:
       DO J=1,NGP
        SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
        Shf0(J) =   Shf0(J)*CP*(TSFC(J) -T0(J))
        dShf(J) = CDENVV(J)*CP
     &            + dShf(J)*CP*(TSKIN(J)-T0(J))
        dShf(J) = MAX( dShf(J), 0. _d 0 )
C--   do not allow negative derivative vs Ts of Sensible+Latent H.flux:
C     a) quiet unrealistic ;
C     b) garantee positive deriv. of total H.flux (needed for implicit solver)
        dEvp(J) = MAX( dEvp(J), -dShf(J)/ALHevp )
       ENDDO
      ELSE
       DO J=1,NGP
        SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
        Shf0(J) = CDENVV(J)*CP*(TSFC(J) -T0(J))
        dShf(J) = CDENVV(J)*CP
       ENDDO
      ENDIF

C     2.4 Emission of lw radiation from the surface

      DO J=1,NGP
        TS2     = TSFC(J)*TSFC(J)
        Slr0(J) = SBC*TS2*TS2
        TS2     = TSKIN(J)*TSKIN(J)
        SLRU(J) = SBC*TS2*TS2
        dSlr(J)  = 4. _d 0 *SBC*TS2*TSKIN(J)
      ENDDO

C--   Compute net surface heat flux and its derivative ./. surf. temp.
      DO J=1,NGP
        sFlx(J,0)= ( SLRD(J) - EMISloc*Slr0(J) )
     &           - ( Shf0(J) + ALHevp*Evp0(J) )
        sFlx(J,1)= ( SLRD(J) - EMISloc*SLRU(J) )
     &           - ( SHF(J)  + ALHevp*EVAP(J) )
        sFlx(J,2)=            -EMISloc*dSlr(J)
     &           - ( dShf(J) + ALHevp*dEvp(J) )
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_FORMULA_AIM */

      RETURN
      END
