C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_aim.F,v 1.2 2006/06/22 14:10:29 jmc Exp $
C $Name: $
#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