pkg/thsice/thsice_get_exf.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_get_exf.F,v 1.10 2007/05/09 12:46:10 jmc 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 ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#ifdef ALLOW_EXF
# 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
      _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     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     copy a few variables to names used in bulkf_formula_lay
      Tsf    =  tsfCel+cen2kel
      Ts2    = Tsf*Tsf
C     wind speed
      ws     = wspeed(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.
        IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN

C--   air - surface difference of temperature & humidity
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

          IF ( useStabilityFct_overIce ) THEN
C--   Compute the turbulent surface fluxes (function of stability).

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

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

           t0     = atemp(i,j,bi,bj)*
     &            (exf_one + humid_fac*aqh(i,j,bi,bj))
           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

          ELSE
C--   Compute the turbulent surface fluxes using fixed transfert Coeffs
C      with no stability dependence ( useStabilityFct_overIce = false )

           evapLoc =      -atmrho*exf_iceCe*wsm*delq
           hlLocal = -lath*evapLoc
           hsLocal = atmcp*atmrho*exf_iceCh*wsm*deltap
C---  surf.Temp derivative of turbulent Fluxes
           dEvdT  = (atmrho*exf_iceCe*wsm)*(ssq*qsat_exp/Ts2)
           dflhdT = -lath*dEvdT
           dfshdT = -atmcp*atmrho*exf_iceCh*wsm

          ENDIF
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         flwNet_dwn  =       lwdown(i,j,bi,bj) - flwup
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.11	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_get_exf.F,v 1.10 2007/05/09 12:46:10 jmc 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	jmc	1.10	#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31	mlosch	1.1	#ifdef ALLOW_EXF
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	_RL czol
76			_RL ws ! wind speed [m/s] (unlimited)
77			_RL wsm ! limited wind speed [m/s] (> umin)
78			_RL t0 ! virtual temperature [K]
79			_RL ustar ! friction velocity [m/s]
80			_RL tstar ! turbulent temperature scale [K]
81			_RL qstar ! turbulent humidity scale [kg/kg]
82			_RL ssq
83			_RL rd ! = sqrt(Cd) [-]
84			_RL re ! = Ce / sqrt(Cd) [-]
85			_RL rh ! = Ch / sqrt(Cd) [-]
86			_RL rdn, ren, rhn ! neutral, zref (=10m) values of rd, re, rh
87			_RL usn, usm ! neutral, zref (=10m) wind-speed (+limited)
88			_RL stable ! = 1 if stable ; = 0 if unstable
89			_RL huol ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
90			_RL htol ! stability parameter at zth [-]
91			_RL hqol
92			_RL x ! stability function [-]
93			_RL xsq ! = x^2 [-]
94			_RL psimh ! momentum stability function
95			_RL psixh ! latent & sensib. stability function
96			_RL zwln ! = log(zwd/zref)
97			_RL ztln ! = log(zth/zref)
98			_RL tau ! surface stress coef = rhoA * Ws * sqrt(Cd)
99			_RL tmpbulk
100	mlosch	1.1
101			C additional variables that are copied from bulkf_formula_lay:
102			C upward LW at surface (W m-2)
103			_RL flwup
104			C net (downward) LW at surface (W m-2)
105			_RL flwNet_dwn
106	jmc	1.3	C gradients of latent/sensible net upward heat flux
107	mlosch	1.1	C w/ respect to temperature
108			_RL dflhdT, dfshdT, dflwupdT
109			C emissivities, called emittance in exf
110			_RL emiss
111	jmc	1.3	C Tsf :: surface temperature [K]
112			C Ts2 :: surface temperature square [K^2]
113			_RL Tsf
114	mlosch	1.1	_RL Ts2
115	jmc	1.3	C latent heat of evaporation or sublimation [J/kg]
116			_RL lath
117	jmc	1.9	_RL qsat_fac, qsat_exp
118	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
119	jmc	1.9	INTEGER ikey_1
120			INTEGER ikey_2
121	heimbach	1.6	#endif
122	mlosch	1.1
123	jmc	1.3	C == external functions ==
124	mlosch	1.1
125	jmc	1.3	c _RL exf_BulkqSat
126			c external exf_BulkqSat
127	jmc	1.9	c _RL exf_BulkCdn
128			c external exf_BulkCdn
129			c _RL exf_BulkRhn
130			c external exf_BulkRhn
131	mlosch	1.1
132	jmc	1.3	C == end of interface ==
133	mlosch	1.1
134	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
135			act1 = bi - myBxLo(myThid)
136			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
137			act2 = bj - myByLo(myThid)
138			max2 = myByHi(myThid) - myByLo(myThid) + 1
139			act3 = myThid - 1
140			max3 = nTx*nTy
141			act4 = ikey_dynamics - 1
142
143			ikey_1 = i
144			& + sNx*(j-1)
145			& + sNxsNyact1
146			& + sNxsNymax1*act2
147			& + sNxsNymax1max2act3
148			& + sNxsNymax1max2max3*act4
149			#endif
150
151	mlosch	1.1	C copy a few variables to names used in bulkf_formula_lay
152	jmc	1.3	Tsf = tsfCel+cen2kel
153			Ts2 = Tsf*Tsf
154	jmc	1.9	C wind speed
155	jmc	1.11	ws = wspeed(i,j,bi,bj)
156	jmc	1.9	#ifdef ALLOW_AUTODIFF_TAMC
157			CADJ STORE sh(i,j,bi,bj) = comlev1_exf_1, key = ikey_1
158			#endif
159			wsm = sh(i,j,bi,bj)
160	jmc	1.3	IF ( iceornot.EQ.0 ) THEN
161			lath = flamb
162	jmc	1.9	qsat_fac = cvapor_fac
163			qsat_exp = cvapor_exp
164	jmc	1.3	ELSE
165			lath = flamb+flami
166	jmc	1.9	qsat_fac = cvapor_fac_ice
167			qsat_exp = cvapor_exp_ice
168	jmc	1.3	ENDIF
169	mlosch	1.1
170	jmc	1.3	C-- Use atmospheric state to compute surface fluxes.
171	jmc	1.10	IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
172	mlosch	1.1
173	jmc	1.10	C-- air - surface difference of temperature & humidity
174			c tmpbulk= exf_BulkqSat(Tsf)
175			c ssq = saltsat*tmpbulk/atmrho
176			tmpbulk= qsat_fac*EXP(-qsat_exp/Tsf)
177			ssq = tmpbulk/atmrho
178			deltap = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
179			delq = aqh(i,j,bi,bj) - ssq
180
181			IF ( useStabilityFct_overIce ) THEN
182			C-- Compute the turbulent surface fluxes (function of stability).
183
184			C-- Set surface parameters :
185			zwln = LOG(hu/zref)
186			ztln = LOG(ht/zref)
187			czol = hukarmangravity_mks
188	mlosch	1.1
189	jmc	1.3	C Initial guess: z/l=0.0; hu=ht=hq=z
190			C Iterations: converge on z/l and hence the fluxes.
191	mlosch	1.1
192			t0 = atemp(i,j,bi,bj)*
193			& (exf_one + humid_fac*aqh(i,j,bi,bj))
194	jmc	1.9	stable = exf_half + SIGN(exf_half, deltap)
195			c tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
196			tmpbulk= cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
197			rdn = SQRT(tmpbulk)
198			C-- initial guess for exchange other coefficients:
199			c rhn = exf_BulkRhn(stable)
200			rhn = (exf_one-stable)cstanton_1 + stablecstanton_2
201			ren = cDalton
202			C-- calculate turbulent scales
203			ustar = rdn*wsm
204			tstar = rhn*deltap
205			qstar = ren*delq
206	mlosch	1.1
207	jmc	1.9	DO iter = 1,niter_bulk
208	mlosch	1.1
209	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
210			ikey_2 = iter
211			& + niter_bulk*(i-1)
212			& + niter_bulksNx(j-1)
213			& + niter_bulksNxsNy*act1
214			& + niter_bulksNxsNymax1act2
215			& + niter_bulksNxsNymax1max2*act3
216			& + niter_bulksNxsNymax1max2max3act4
217
218			CADJ STORE rdn = comlev1_exf_2, key = ikey_2
219			CADJ STORE ustar = comlev1_exf_2, key = ikey_2
220			CADJ STORE qstar = comlev1_exf_2, key = ikey_2
221			CADJ STORE tstar = comlev1_exf_2, key = ikey_2
222			CADJ STORE sh(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
223			#endif
224
225	jmc	1.9	huol = (tstar/t0 +
226			& qstar/(exf_one/humid_fac+aqh(i,j,bi,bj))
227			& )czol/(ustarustar)
228			C- Large&Pond_1981 code (zolmin default = -100):
229			c huol = MAX(huol,zolmin)
230			C- Large&Yeager_2004 code:
231			huol = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
232	mlosch	1.1	htol = huol*ht/hu
233			hqol = huol*hq/hu
234	jmc	1.9	stable = exf_half + SIGN(exf_half, huol)
235	jmc	1.3
236			C Evaluate all stability functions assuming hq = ht.
237	jmc	1.9	C- Large&Pond_1981 code:
238			c xsq = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)
239			C- Large&Yeager_2004 code:
240			xsq = SQRT( ABS(exf_one - huol*16. _d 0) )
241			x = SQRT(xsq)
242			psimh = -psim_fachuolstable
243			& + (exf_one-stable)
244			& ( LOG( (exf_one + exf_twox + xsq)
245			& (exf_one+xsq)0.125 _d 0 )
246			& -exf_twoATAN(x) + exf_halfpi )
247			C- Large&Pond_1981 code:
248			c xsq = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
249			C- Large&Yeager_2004 code:
250			xsq = SQRT( ABS(exf_one - htol*16. _d 0) )
251			psixh = -psim_fachtolstable
252			& + (exf_one-stable)
253			& exf_twoLOG( exf_half*(exf_one+xsq) )
254	jmc	1.3
255			C Shift wind speed using old coefficient
256	jmc	1.9	usn = ws/( exf_one + rdn*(zwln-psimh)/karman )
257			usm = MAX(usn, umin)
258
259			C- Update the 10m, neutral stability transfer coefficients
260			c tmpbulk= exf_BulkCdn(usm)
261			tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm
262			rdn = SQRT(tmpbulk)
263			c rhn = exf_BulkRhn(stable)
264			rhn = (exf_one-stable)cstanton_1 + stablecstanton_2
265
266			C Shift all coefficients to the measurement height and stability.
267			rd = rdn/( exf_one + rdn*(zwln-psimh)/karman )
268			rh = rhn/( exf_one + rhn*(ztln-psixh)/karman )
269			re = ren/( exf_one + ren*(ztln-psixh)/karman )
270
271			C Update ustar, tstar, qstar using updated, shifted coefficients.
272			ustar = rd*wsm
273	jmc	1.3	qstar = re*delq
274	mlosch	1.1	tstar = rh*deltap
275	jmc	1.9	ENDDO
276	mlosch	1.1
277	jmc	1.9	tau = atmrhordws
278	mlosch	1.1
279	jmc	1.9	evapLoc = -tau*qstar
280	jmc	1.4	hlLocal = -lath*evapLoc
281	jmc	1.9	hsLocal = atmcptautstar
282			c ustress = taurdUwindSpeed
283			c vstress = taurdVwindSpeed
284	mlosch	1.1
285			C--- surf.Temp derivative of turbulent Fluxes
286	jmc	1.9	dEvdT = (taure)ssq*qsat_exp/Ts2
287	mlosch	1.1	dflhdT = -lath*dEvdT
288	jmc	1.9	dfshdT = -atmcptaurh
289	mlosch	1.1
290	jmc	1.10	ELSE
291			C-- Compute the turbulent surface fluxes using fixed transfert Coeffs
292			C with no stability dependence ( useStabilityFct_overIce = false )
293
294			evapLoc = -atmrhoexf_iceCewsm*delq
295			hlLocal = -lath*evapLoc
296			hsLocal = atmcpatmrhoexf_iceChwsmdeltap
297			C--- surf.Temp derivative of turbulent Fluxes
298			dEvdT = (atmrhoexf_iceCewsm)(ssqqsat_exp/Ts2)
299			dflhdT = -lath*dEvdT
300			dfshdT = -atmcpatmrhoexf_iceCh*wsm
301
302			ENDIF
303	mlosch	1.1	C--- Upward long wave radiation
304	jmc	1.10	IF ( iceornot.EQ.0 ) THEN
305	mlosch	1.1	emiss = ocean_emissivity
306	jmc	1.10	ELSEIF (iceornot.EQ.2) THEN
307	mlosch	1.1	emiss = snow_emissivity
308	jmc	1.10	ELSE
309	mlosch	1.1	emiss = ice_emissivity
310	jmc	1.10	ENDIF
311			flwup = emissstefanBoltzmannTs2*Ts2
312			dflwupdT = emissstefanBoltzmannTs2Tsf 4. _d 0
313	jmc	1.3
314	mlosch	1.1	C-- Total derivative with respect to surface temperature
315	jmc	1.10	df0dT = -dflwupdT+dfshdT+dflhdT
316	jmc	1.3
317	jmc	1.7	#ifdef ALLOW_DOWNWARD_RADIATION
318	jmc	1.10	c flwNet_dwn = lwdown(i,j,bi,bj) - flwup
319			C- assume long-wave albedo = 1 - emissivity :
320			flwNet_dwn = emiss*lwdown(i,j,bi,bj) - flwup
321	jmc	1.7	#else
322			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
323			#endif
324	jmc	1.10	flxExceptSw = flwNet_dwn + hsLocal + hlLocal
325	mlosch	1.1
326	jmc	1.10	ELSE
327			flxExceptSw = 0. _d 0
328			df0dT = 0. _d 0
329			evapLoc = 0. _d 0
330			dEvdT = 0. _d 0
331			ENDIF
332	mlosch	1.1
333			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
334
335	jmc	1.7	#else /* ALLOW_ATM_TEMP */
336			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'
337			#endif /* ALLOW_ATM_TEMP */
338	jmc	1.9	#ifdef EXF_READ_EVAP
339			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: EXF_READ_EVAP defined'
340			#endif /* EXF_READ_EVAP */
341	mlosch	1.1	#endif /* ALLOW_EXF */
342
343			RETURN
344			END