pkg/thsice/thsice_get_exf.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_get_exf.F,v 1.8 2007/05/03 21:51:12 mlosch Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"
#ifdef ALLOW_EXF
#include "EXF_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: THSICE_GET_EXF
C     !INTERFACE:
      SUBROUTINE THSICE_GET_EXF(
     I                         iceornot, tsfCel,
     O                         flxExceptSw, df0dT, evapLoc, dEvdT,
     I                         i,j,bi,bj,myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R  THSICE_GET_EXF
C     *==========================================================*
C     | Interface S/R : get Surface Fluxes from pkg EXF
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     == Global data ==
#ifdef ALLOW_EXF
# include "SIZE.h"
# include "EEPARAMS.h"
# include "PARAMS.h"
# include "EXF_CONSTANTS.h"
# include "EXF_PARAM.h"
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     iceornot    :: 0=open water, 1=ice cover
C     tsfCel      :: surface (ice or snow) temperature (oC)
C     flxExceptSw :: net (downward) surface heat flux, except short-wave [W/m2]
C     df0dT       :: deriv of flx with respect to Tsf    [W/m/K]
C     evapLoc     :: surface evaporation (>0 if evaporate) [kg/m2/s]
C     dEvdT       :: deriv of evap. with respect to Tsf  [kg/m2/s/K]
C     i,j, bi,bj  :: current grid point indices
C     myThid      :: My Thread Id number
      INTEGER i,j, bi,bj
      INTEGER myThid
      INTEGER iceornot
      _RL  tsfCel
      _RL  flxExceptSw
      _RL  df0dT
      _RL  evapLoc
      _RL  dEvdT
CEOP

#ifdef ALLOW_EXF
#ifdef ALLOW_ATM_TEMP

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C     === Local variables ===
C     hsLocal, hlLocal :: sensible & latent heat flux over sea-ice
C             t0       :: virtual temperature (K)
C             ssq      :: saturation specific humidity (kg/kg)
C             deltap   :: potential temperature diff (K)

      _RL     hsLocal, hlLocal
      INTEGER iter
      _RL     delq
      _RL     deltap
c-----------------------
      _RL czol
      _RL ws                 ! wind speed [m/s] (unlimited)
      _RL wsm                ! limited wind speed [m/s] (> umin)
      _RL t0                 ! virtual temperature [K]
      _RL ustar              ! friction velocity [m/s]
      _RL tstar              ! turbulent temperature scale [K]
      _RL qstar              ! turbulent humidity scale  [kg/kg]
      _RL ssq
      _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 usn, usm           ! neutral, zref (=10m) wind-speed (+limited)
      _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 hqol
      _RL x                  ! stability function  [-]
      _RL xsq                ! = x^2               [-]
      _RL psimh              ! momentum stability function
      _RL psixh              ! latent & sensib. stability function
      _RL zwln               ! = log(zwd/zref)
      _RL ztln               ! = log(zth/zref)
      _RL tau                ! surface stress coef = rhoA * Ws * sqrt(Cd)
      _RL tmpbulk
c-----------------------

C     additional variables that are copied from bulkf_formula_lay:
C     upward LW at surface (W m-2)
      _RL  flwup
C     net (downward) LW at surface (W m-2)
      _RL  flwNet_dwn
C     gradients of latent/sensible net upward heat flux
C     w/ respect to temperature
      _RL dflhdT, dfshdT, dflwupdT
C     emissivities, called emittance in exf
      _RL     emiss
C     Tsf    :: surface temperature [K]
C     Ts2    :: surface temperature square [K^2]
      _RL Tsf
      _RL Ts2
C     latent heat of evaporation or sublimation [J/kg]
      _RL lath
      _RL qsat_fac, qsat_exp
#ifdef ALLOW_AUTODIFF_TAMC
      INTEGER ikey_1
      INTEGER ikey_2
#endif

C     == external functions ==

c     _RL       exf_BulkqSat
c     external  exf_BulkqSat
c     _RL       exf_BulkCdn
c     external  exf_BulkCdn
c     _RL       exf_BulkRhn
c     external  exf_BulkRhn

C     == end of interface ==

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1

          ikey_1 = i
     &         + sNx*(j-1)
     &         + sNx*sNy*act1
     &         + sNx*sNy*max1*act2
     &         + sNx*sNy*max1*max2*act3
     &         + sNx*sNy*max1*max2*max3*act4
#endif

C-- Set surface parameters :
      zwln = LOG(hu/zref)
      ztln = LOG(ht/zref)
      czol = hu*karman*gravity_mks

C     copy a few variables to names used in bulkf_formula_lay
      Tsf    =  tsfCel+cen2kel
      Ts2    = Tsf*Tsf
C     wind speed
      ws     = us(i,j,bi,bj)
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
#endif
      wsm    = sh(i,j,bi,bj)
      IF ( iceornot.EQ.0 ) THEN
        lath  = flamb
        qsat_fac = cvapor_fac
        qsat_exp = cvapor_exp
      ELSE
        lath  = flamb+flami
        qsat_fac = cvapor_fac_ice
        qsat_exp = cvapor_exp_ice
      ENDIF

C--   Use atmospheric state to compute surface fluxes.

C--   Compute the turbulent surface fluxes.

C             Initial guess: z/l=0.0; hu=ht=hq=z
C             Iterations:    converge on z/l and hence the fluxes.

          IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
           t0     = atemp(i,j,bi,bj)*
     &            (exf_one + humid_fac*aqh(i,j,bi,bj))
c          tmpbulk= exf_BulkqSat(Tsf)
c          ssq    = saltsat*tmpbulk/atmrho
           tmpbulk= qsat_fac*EXP(-qsat_exp/Tsf)
           ssq    = tmpbulk/atmrho
           deltap = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
           delq   = aqh(i,j,bi,bj) - ssq
           stable = exf_half + SIGN(exf_half, deltap)
c          tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
           tmpbulk= cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
           rdn    = SQRT(tmpbulk)
C--  initial guess for exchange other coefficients:
c          rhn    = exf_BulkRhn(stable)
           rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2
           ren    = cDalton
C--  calculate turbulent scales
           ustar  = rdn*wsm
           tstar  = rhn*deltap
           qstar  = ren*delq

           DO iter = 1,niter_bulk

#ifdef ALLOW_AUTODIFF_TAMC
                   ikey_2 = iter
     &                  + niter_bulk*(i-1)
     &                  + niter_bulk*sNx*(j-1)
     &                  + niter_bulk*sNx*sNy*act1
     &                  + niter_bulk*sNx*sNy*max1*act2
     &                  + niter_bulk*sNx*sNy*max1*max2*act3
     &                  + niter_bulk*sNx*sNy*max1*max2*max3*act4

CADJ STORE rdn    = comlev1_exf_2, key = ikey_2
CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
CADJ STORE qstar  = comlev1_exf_2, key = ikey_2
CADJ STORE tstar  = comlev1_exf_2, key = ikey_2
CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2
#endif

            huol   = (tstar/t0 +
     &                qstar/(exf_one/humid_fac+aqh(i,j,bi,bj))
     &               )*czol/(ustar*ustar)
C-    Large&Pond_1981 code (zolmin default = -100):
c           huol   = MAX(huol,zolmin)
C-    Large&Yeager_2004 code:
            huol   = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
            htol   = huol*ht/hu
            hqol   = huol*hq/hu
            stable = exf_half + SIGN(exf_half, huol)

C     Evaluate all stability functions assuming hq = ht.
C-    Large&Pond_1981 code:
c           xsq    = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)
C-    Large&Yeager_2004 code:
            xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
            x      = SQRT(xsq)
            psimh  = -psim_fac*huol*stable
     &             + (exf_one-stable)
     &              *( LOG( (exf_one + exf_two*x + xsq)
     &                     *(exf_one+xsq)*0.125 _d 0 )
     &                 -exf_two*ATAN(x) + exf_half*pi )
C-    Large&Pond_1981 code:
c           xsq    = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
C-    Large&Yeager_2004 code:
            xsq    = SQRT( ABS(exf_one - htol*16. _d 0) )
            psixh  = -psim_fac*htol*stable
     &             + (exf_one-stable)
     &              *exf_two*LOG( exf_half*(exf_one+xsq) )

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

C-    Update the 10m, neutral stability transfer coefficients
c           tmpbulk= exf_BulkCdn(usm)
            tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm
            rdn    = SQRT(tmpbulk)
c           rhn    = exf_BulkRhn(stable)
            rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2

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

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

           tau     = atmrho*rd*ws

           evapLoc = -tau*qstar
           hlLocal = -lath*evapLoc
           hsLocal = atmcp*tau*tstar
c          ustress = tau*rd*UwindSpeed
c          vstress = tau*rd*VwindSpeed

C---  surf.Temp derivative of turbulent Fluxes
           dEvdT  =  (tau*re)*ssq*qsat_exp/Ts2
           dflhdT = -lath*dEvdT
           dfshdT = -atmcp*tau*rh

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

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

#ifdef ALLOW_DOWNWARD_RADIATION
C-    assume long-wave albedo = 1 - emissivity
           flwNet_dwn  = emiss*lwdown(i,j,bi,bj) - flwup
#else
      STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
#endif
           flxExceptSw = flwNet_dwn + hsLocal + hlLocal

          ELSE
           flxExceptSw = 0. _d 0
           df0dT   = 0. _d 0
           evapLoc = 0. _d 0
           dEvdT   = 0. _d 0
          ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#else /* ALLOW_ATM_TEMP */
      STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'
#endif /* ALLOW_ATM_TEMP */
#ifdef EXF_READ_EVAP
      STOP 'ABNORMAL END: S/R THSICE_GET_EXF: EXF_READ_EVAP defined'
#endif /* EXF_READ_EVAP */
#endif /* ALLOW_EXF */

      RETURN
      END
1	jmc	1.9	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_get_exf.F,v 1.8 2007/05/03 21:51:12 mlosch Exp $
2	mlosch	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5			#ifdef ALLOW_EXF
6			#include "EXF_OPTIONS.h"
7			#endif
8
9			CBOP
10			C !ROUTINE: THSICE_GET_EXF
11			C !INTERFACE:
12			SUBROUTINE THSICE_GET_EXF(
13	jmc	1.3	I iceornot, tsfCel,
14	mlosch	1.1	O flxExceptSw, df0dT, evapLoc, dEvdT,
15			I i,j,bi,bj,myThid )
16			C !DESCRIPTION: \bv
17			C ==========================================================
18			C \| S/R THSICE_GET_EXF
19			C ==========================================================
20			C \| Interface S/R : get Surface Fluxes from pkg EXF
21			C ==========================================================
22			C \ev
23
24			C !USES:
25			IMPLICIT NONE
26
27			C == Global data ==
28			#ifdef ALLOW_EXF
29			# include "SIZE.h"
30			# include "EEPARAMS.h"
31			# include "PARAMS.h"
32	jmc	1.5	# include "EXF_CONSTANTS.h"
33			# include "EXF_PARAM.h"
34			# include "EXF_FIELDS.h"
35	mlosch	1.1	#endif
36	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
37			# include "tamc.h"
38			# include "tamc_keys.h"
39			#endif
40	mlosch	1.1
41			C !INPUT/OUTPUT PARAMETERS:
42			C === Routine arguments ===
43			C iceornot :: 0=open water, 1=ice cover
44	jmc	1.3	C tsfCel :: surface (ice or snow) temperature (oC)
45	mlosch	1.1	C flxExceptSw :: net (downward) surface heat flux, except short-wave [W/m2]
46			C df0dT :: deriv of flx with respect to Tsf [W/m/K]
47			C evapLoc :: surface evaporation (>0 if evaporate) [kg/m2/s]
48			C dEvdT :: deriv of evap. with respect to Tsf [kg/m2/s/K]
49			C i,j, bi,bj :: current grid point indices
50	jmc	1.9	C myThid :: My Thread Id number
51	mlosch	1.1	INTEGER i,j, bi,bj
52			INTEGER myThid
53			INTEGER iceornot
54	jmc	1.3	_RL tsfCel
55	mlosch	1.1	_RL flxExceptSw
56			_RL df0dT
57			_RL evapLoc
58			_RL dEvdT
59			CEOP
60
61			#ifdef ALLOW_EXF
62	jmc	1.7	#ifdef ALLOW_ATM_TEMP
63	mlosch	1.1
64			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
65			C === Local variables ===
66	jmc	1.4	C hsLocal, hlLocal :: sensible & latent heat flux over sea-ice
67	jmc	1.9	C t0 :: virtual temperature (K)
68			C ssq :: saturation specific humidity (kg/kg)
69			C deltap :: potential temperature diff (K)
70	mlosch	1.1
71	jmc	1.4	_RL hsLocal, hlLocal
72	jmc	1.9	INTEGER iter
73	mlosch	1.1	_RL delq
74			_RL deltap
75	jmc	1.9	c-----------------------
76			_RL czol
77			_RL ws ! wind speed [m/s] (unlimited)
78			_RL wsm ! limited wind speed [m/s] (> umin)
79			_RL t0 ! virtual temperature [K]
80			_RL ustar ! friction velocity [m/s]
81			_RL tstar ! turbulent temperature scale [K]
82			_RL qstar ! turbulent humidity scale [kg/kg]
83			_RL ssq
84			_RL rd ! = sqrt(Cd) [-]
85			_RL re ! = Ce / sqrt(Cd) [-]
86			_RL rh ! = Ch / sqrt(Cd) [-]
87			_RL rdn, ren, rhn ! neutral, zref (=10m) values of rd, re, rh
88			_RL usn, usm ! neutral, zref (=10m) wind-speed (+limited)
89			_RL stable ! = 1 if stable ; = 0 if unstable
90			_RL huol ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
91			_RL htol ! stability parameter at zth [-]
92			_RL hqol
93			_RL x ! stability function [-]
94			_RL xsq ! = x^2 [-]
95			_RL psimh ! momentum stability function
96			_RL psixh ! latent & sensib. stability function
97			_RL zwln ! = log(zwd/zref)
98			_RL ztln ! = log(zth/zref)
99			_RL tau ! surface stress coef = rhoA * Ws * sqrt(Cd)
100			_RL tmpbulk
101			c-----------------------
102	mlosch	1.1
103			C additional variables that are copied from bulkf_formula_lay:
104			C upward LW at surface (W m-2)
105			_RL flwup
106			C net (downward) LW at surface (W m-2)
107			_RL flwNet_dwn
108	jmc	1.3	C gradients of latent/sensible net upward heat flux
109	mlosch	1.1	C w/ respect to temperature
110			_RL dflhdT, dfshdT, dflwupdT
111			C emissivities, called emittance in exf
112			_RL emiss
113	jmc	1.3	C Tsf :: surface temperature [K]
114			C Ts2 :: surface temperature square [K^2]
115			_RL Tsf
116	mlosch	1.1	_RL Ts2
117	jmc	1.3	C latent heat of evaporation or sublimation [J/kg]
118			_RL lath
119	jmc	1.9	_RL qsat_fac, qsat_exp
120	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
121	jmc	1.9	INTEGER ikey_1
122			INTEGER ikey_2
123	heimbach	1.6	#endif
124	mlosch	1.1
125	jmc	1.3	C == external functions ==
126	mlosch	1.1
127	jmc	1.3	c _RL exf_BulkqSat
128			c external exf_BulkqSat
129	jmc	1.9	c _RL exf_BulkCdn
130			c external exf_BulkCdn
131			c _RL exf_BulkRhn
132			c external exf_BulkRhn
133	mlosch	1.1
134	jmc	1.3	C == end of interface ==
135	mlosch	1.1
136	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
137			act1 = bi - myBxLo(myThid)
138			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
139			act2 = bj - myByLo(myThid)
140			max2 = myByHi(myThid) - myByLo(myThid) + 1
141			act3 = myThid - 1
142			max3 = nTx*nTy
143			act4 = ikey_dynamics - 1
144
145			ikey_1 = i
146			& + sNx*(j-1)
147			& + sNxsNyact1
148			& + sNxsNymax1*act2
149			& + sNxsNymax1max2act3
150			& + sNxsNymax1max2max3*act4
151			#endif
152
153	mlosch	1.8	C-- Set surface parameters :
154	jmc	1.9	zwln = LOG(hu/zref)
155			ztln = LOG(ht/zref)
156	mlosch	1.8	czol = hukarmangravity_mks
157
158	mlosch	1.1	C copy a few variables to names used in bulkf_formula_lay
159	jmc	1.3	Tsf = tsfCel+cen2kel
160			Ts2 = Tsf*Tsf
161	jmc	1.9	C wind speed
162			ws = us(i,j,bi,bj)
163			#ifdef ALLOW_AUTODIFF_TAMC
164			CADJ STORE sh(i,j,bi,bj) = comlev1_exf_1, key = ikey_1
165			#endif
166			wsm = sh(i,j,bi,bj)
167	jmc	1.3	IF ( iceornot.EQ.0 ) THEN
168			lath = flamb
169	jmc	1.9	qsat_fac = cvapor_fac
170			qsat_exp = cvapor_exp
171	jmc	1.3	ELSE
172			lath = flamb+flami
173	jmc	1.9	qsat_fac = cvapor_fac_ice
174			qsat_exp = cvapor_exp_ice
175	jmc	1.3	ENDIF
176	mlosch	1.1
177	jmc	1.3	C-- Use atmospheric state to compute surface fluxes.
178	mlosch	1.1
179	jmc	1.3	C-- Compute the turbulent surface fluxes.
180	mlosch	1.1
181	jmc	1.3	C Initial guess: z/l=0.0; hu=ht=hq=z
182			C Iterations: converge on z/l and hence the fluxes.
183	mlosch	1.1
184	jmc	1.9	IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
185	mlosch	1.1	t0 = atemp(i,j,bi,bj)*
186			& (exf_one + humid_fac*aqh(i,j,bi,bj))
187	jmc	1.3	c tmpbulk= exf_BulkqSat(Tsf)
188			c ssq = saltsat*tmpbulk/atmrho
189	jmc	1.9	tmpbulk= qsat_fac*EXP(-qsat_exp/Tsf)
190	jmc	1.3	ssq = tmpbulk/atmrho
191	jmc	1.9	deltap = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
192	mlosch	1.1	delq = aqh(i,j,bi,bj) - ssq
193	jmc	1.9	stable = exf_half + SIGN(exf_half, deltap)
194			c tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
195			tmpbulk= cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
196			rdn = SQRT(tmpbulk)
197			C-- initial guess for exchange other coefficients:
198			c rhn = exf_BulkRhn(stable)
199			rhn = (exf_one-stable)cstanton_1 + stablecstanton_2
200			ren = cDalton
201			C-- calculate turbulent scales
202			ustar = rdn*wsm
203			tstar = rhn*deltap
204			qstar = ren*delq
205	mlosch	1.1
206	jmc	1.9	DO iter = 1,niter_bulk
207	mlosch	1.1
208	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
209			ikey_2 = iter
210			& + niter_bulk*(i-1)
211			& + niter_bulksNx(j-1)
212			& + niter_bulksNxsNy*act1
213			& + niter_bulksNxsNymax1act2
214			& + niter_bulksNxsNymax1max2*act3
215			& + niter_bulksNxsNymax1max2max3act4
216
217			CADJ STORE rdn = comlev1_exf_2, key = ikey_2
218			CADJ STORE ustar = comlev1_exf_2, key = ikey_2
219			CADJ STORE qstar = comlev1_exf_2, key = ikey_2
220			CADJ STORE tstar = comlev1_exf_2, key = ikey_2
221			CADJ STORE sh(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
222			#endif
223
224	jmc	1.9	huol = (tstar/t0 +
225			& qstar/(exf_one/humid_fac+aqh(i,j,bi,bj))
226			& )czol/(ustarustar)
227			C- Large&Pond_1981 code (zolmin default = -100):
228			c huol = MAX(huol,zolmin)
229			C- Large&Yeager_2004 code:
230			huol = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
231	mlosch	1.1	htol = huol*ht/hu
232			hqol = huol*hq/hu
233	jmc	1.9	stable = exf_half + SIGN(exf_half, huol)
234	jmc	1.3
235			C Evaluate all stability functions assuming hq = ht.
236	jmc	1.9	C- Large&Pond_1981 code:
237			c xsq = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)
238			C- Large&Yeager_2004 code:
239			xsq = SQRT( ABS(exf_one - huol*16. _d 0) )
240			x = SQRT(xsq)
241			psimh = -psim_fachuolstable
242			& + (exf_one-stable)
243			& ( LOG( (exf_one + exf_twox + xsq)
244			& (exf_one+xsq)0.125 _d 0 )
245			& -exf_twoATAN(x) + exf_halfpi )
246			C- Large&Pond_1981 code:
247			c xsq = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
248			C- Large&Yeager_2004 code:
249			xsq = SQRT( ABS(exf_one - htol*16. _d 0) )
250			psixh = -psim_fachtolstable
251			& + (exf_one-stable)
252			& exf_twoLOG( exf_half*(exf_one+xsq) )
253	jmc	1.3
254			C Shift wind speed using old coefficient
255	jmc	1.9	usn = ws/( exf_one + rdn*(zwln-psimh)/karman )
256			usm = MAX(usn, umin)
257
258			C- Update the 10m, neutral stability transfer coefficients
259			c tmpbulk= exf_BulkCdn(usm)
260			tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm
261			rdn = SQRT(tmpbulk)
262			c rhn = exf_BulkRhn(stable)
263			rhn = (exf_one-stable)cstanton_1 + stablecstanton_2
264
265			C Shift all coefficients to the measurement height and stability.
266			rd = rdn/( exf_one + rdn*(zwln-psimh)/karman )
267			rh = rhn/( exf_one + rhn*(ztln-psixh)/karman )
268			re = ren/( exf_one + ren*(ztln-psixh)/karman )
269
270			C Update ustar, tstar, qstar using updated, shifted coefficients.
271			ustar = rd*wsm
272	jmc	1.3	qstar = re*delq
273	mlosch	1.1	tstar = rh*deltap
274	jmc	1.9	ENDDO
275	mlosch	1.1
276	jmc	1.9	tau = atmrhordws
277	mlosch	1.1
278	jmc	1.9	evapLoc = -tau*qstar
279	jmc	1.4	hlLocal = -lath*evapLoc
280	jmc	1.9	hsLocal = atmcptautstar
281			c ustress = taurdUwindSpeed
282			c vstress = taurdVwindSpeed
283	mlosch	1.1
284			C--- surf.Temp derivative of turbulent Fluxes
285	jmc	1.9	dEvdT = (taure)ssq*qsat_exp/Ts2
286	mlosch	1.1	dflhdT = -lath*dEvdT
287	jmc	1.9	dfshdT = -atmcptaurh
288	mlosch	1.1
289			C--- Upward long wave radiation
290			IF ( iceornot.EQ.0 ) THEN
291			emiss = ocean_emissivity
292			ELSEIF (iceornot.EQ.2) THEN
293			emiss = snow_emissivity
294			ELSE
295			emiss = ice_emissivity
296			ENDIF
297			flwup = emissstefanBoltzmannTs2*Ts2
298	jmc	1.3	dflwupdT = emissstefanBoltzmannTs2Tsf 4. _d 0
299
300	mlosch	1.1	C-- Total derivative with respect to surface temperature
301			df0dT = -dflwupdT+dfshdT+dflhdT
302	jmc	1.3
303	jmc	1.7	#ifdef ALLOW_DOWNWARD_RADIATION
304	jmc	1.9	C- assume long-wave albedo = 1 - emissivity
305			flwNet_dwn = emiss*lwdown(i,j,bi,bj) - flwup
306	jmc	1.7	#else
307			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
308			#endif
309	jmc	1.4	flxExceptSw = flwNet_dwn + hsLocal + hlLocal
310	mlosch	1.1
311	jmc	1.9	ELSE
312			flxExceptSw = 0. _d 0
313			df0dT = 0. _d 0
314			evapLoc = 0. _d 0
315			dEvdT = 0. _d 0
316			ENDIF
317	mlosch	1.1
318			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
319
320	jmc	1.7	#else /* ALLOW_ATM_TEMP */
321			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'
322			#endif /* ALLOW_ATM_TEMP */
323	jmc	1.9	#ifdef EXF_READ_EVAP
324			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: EXF_READ_EVAP defined'
325			#endif /* EXF_READ_EVAP */
326	mlosch	1.1	#endif /* ALLOW_EXF */
327
328			RETURN
329			END