pkg/thsice/thsice_get_exf.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_get_exf.F,v 1.7 2007/04/27 15:51:28 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 ==
#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      :: Thread no. that called this routine.
      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

      _RL     aln

      _RL     hsLocal, hlLocal
      integer iter
      _RL     delq
      _RL     deltap
      _RL     hqol
      _RL     htol
      _RL     huol
      _RL     psimh
      _RL     psixh
      _RL     qstar
      _RL     rd
      _RL     re
      _RL     rdn
      _RL     rh
      _RL     ssq
      _RL     stable
      _RL     tstar
      _RL     t0
      _RL     ustar
      _RL     uzn
      _RL     shn
      _RL     xsq
      _RL     x
      _RL     tau
      _RL     tmpbulk
      _RL     czol

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
#ifdef ALLOW_AUTODIFF_TAMC
      integer ikey_1
      integer ikey_2
#endif

C     == external functions ==

c     _RL       exf_BulkqSat
c     external  exf_BulkqSat
      _RL       exf_BulkCdn
      external  exf_BulkCdn
      _RL       exf_BulkRhn
      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 :
      czol = hu*karman*gravity_mks

C     copy a few variables to names used in bulkf_formula_lay
      Tsf    =  tsfCel+cen2kel
      Ts2    = Tsf*Tsf
      IF ( iceornot.EQ.0 ) THEN
        lath  = flamb
        dEvdT = cvapor_exp
      ELSE
        lath  = flamb+flami
        dEvdT = cvapor_exp_ice
      ENDIF

Cph   This statement cannot be a PARAMETER statement in the header,
Cph   but must come here; it is not fortran77 standard
      aln = log(ht/zref)

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.
C             t0     : virtual temperature (K)
C             ssq    : sea surface humidity (kg/kg)
C             deltap : potential temperature diff (K)

          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 = cvapor_fac_ice/exp(cvapor_exp_ice/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)
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
#endif
           tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
           rdn    = sqrt(tmpbulk)
           ustar  = rdn*sh(i,j,bi,bj)
           tmpbulk= exf_BulkRhn(stable)
           tstar  = tmpbulk*deltap
           qstar  = cdalton*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   = czol*(tstar/t0 +
     &           qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
     &           ustar**2
            huol   = max(huol,zolmin)
            stable = exf_half + sign(exf_half, huol)
            htol   = huol*ht/hu
            hqol   = huol*hq/hu

C     Evaluate all stability functions assuming hq = ht.
            xsq    = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
            x     = sqrt(xsq)
            psimh  = -psim_fac*huol*stable +
     &           (exf_one - stable)*
     &           (log((exf_one + x*(exf_two + x))*
     &           (exf_one + xsq)/8.) - exf_two*atan(x) +
     &           pi*exf_half)
            xsq    = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
            psixh  = -psim_fac*htol*stable + (exf_one - stable)*
     &           exf_two*log((exf_one + xsq)/exf_two)

C     Shift wind speed using old coefficient
ccc   rd   = rdn/(exf_one + rdn/karman*
ccc   &                 (log(hu/zref) - psimh) )
            rd     = rdn/(exf_one - rdn/karman*psimh )
            shn    = sh(i,j,bi,bj)*rd/rdn
            uzn    = max(shn, umin)

C     Update the transfer coefficients at 10 meters
C     and neutral stability.

            tmpbulk= exf_BulkCdn(uzn)
            rdn    = sqrt(tmpbulk)

C     Shift all coefficients to the measurement height
C     and stability.
c     rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))
            rd     = rdn/(exf_one - rdn/karman*psimh)
            tmpbulk= exf_BulkRhn(stable)
            rh     = tmpbulk/( exf_one +
     &           tmpbulk/karman*(aln - psixh) )
            re     = cdalton/( exf_one +
     &           cdalton/karman*(aln - psixh) )

C     Update ustar, tstar, qstar using updated, shifted
C     coefficients.
            ustar  = rd*sh(i,j,bi,bj)
            qstar  = re*delq
            tstar  = rh*deltap
           enddo

           tau    = atmrho*ustar**2
           tau    = tau*us(i,j,bi,bj)/sh(i,j,bi,bj)

           evapLoc       = -tau*qstar/ustar
           hlLocal = -lath*evapLoc
           hsLocal = atmcp*tau*tstar/ustar
#ifndef EXF_READ_EVAP
cdm             evap(i,j,bi,bj)    = tau*qstar/ustar
cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
C jmc: do not reset evap which contains evaporation over ice-free ocean fraction
c          evap(i,j,bi,bj)    = -recip_rhonil*evapLoc
#endif

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

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
           flwNet_dwn  = lwdown(i,j,bi,bj) - flwup
#else
      STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
#endif
           flxExceptSw = flwNet_dwn + hsLocal + hlLocal

          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 */
#endif /* ALLOW_EXF */

      RETURN
      END
1	mlosch	1.8	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_get_exf.F,v 1.7 2007/04/27 15:51:28 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			#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			C myThid :: Thread no. that called this routine.
51			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	mlosch	1.1
68			_RL aln
69
70	jmc	1.4	_RL hsLocal, hlLocal
71	mlosch	1.1	integer iter
72			_RL delq
73			_RL deltap
74			_RL hqol
75			_RL htol
76			_RL huol
77			_RL psimh
78			_RL psixh
79			_RL qstar
80			_RL rd
81			_RL re
82			_RL rdn
83			_RL rh
84			_RL ssq
85			_RL stable
86			_RL tstar
87			_RL t0
88			_RL ustar
89			_RL uzn
90			_RL shn
91			_RL xsq
92			_RL x
93			_RL tau
94			_RL tmpbulk
95	mlosch	1.8	_RL czol
96	mlosch	1.1
97			C additional variables that are copied from bulkf_formula_lay:
98			C upward LW at surface (W m-2)
99			_RL flwup
100			C net (downward) LW at surface (W m-2)
101			_RL flwNet_dwn
102	jmc	1.3	C gradients of latent/sensible net upward heat flux
103	mlosch	1.1	C w/ respect to temperature
104			_RL dflhdT, dfshdT, dflwupdT
105			C emissivities, called emittance in exf
106			_RL emiss
107	jmc	1.3	C Tsf :: surface temperature [K]
108			C Ts2 :: surface temperature square [K^2]
109			_RL Tsf
110	mlosch	1.1	_RL Ts2
111	jmc	1.3	C latent heat of evaporation or sublimation [J/kg]
112			_RL lath
113	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
114			integer ikey_1
115			integer ikey_2
116			#endif
117	mlosch	1.1
118	jmc	1.3	C == external functions ==
119	mlosch	1.1
120	jmc	1.3	c _RL exf_BulkqSat
121			c external exf_BulkqSat
122	mlosch	1.1	_RL exf_BulkCdn
123			external exf_BulkCdn
124			_RL exf_BulkRhn
125			external exf_BulkRhn
126
127	jmc	1.3	C == end of interface ==
128	mlosch	1.1
129	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
130			act1 = bi - myBxLo(myThid)
131			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
132			act2 = bj - myByLo(myThid)
133			max2 = myByHi(myThid) - myByLo(myThid) + 1
134			act3 = myThid - 1
135			max3 = nTx*nTy
136			act4 = ikey_dynamics - 1
137
138			ikey_1 = i
139			& + sNx*(j-1)
140			& + sNxsNyact1
141			& + sNxsNymax1*act2
142			& + sNxsNymax1max2act3
143			& + sNxsNymax1max2max3*act4
144			#endif
145
146	mlosch	1.8	C-- Set surface parameters :
147			czol = hukarmangravity_mks
148
149	mlosch	1.1	C copy a few variables to names used in bulkf_formula_lay
150	jmc	1.3	Tsf = tsfCel+cen2kel
151			Ts2 = Tsf*Tsf
152			IF ( iceornot.EQ.0 ) THEN
153			lath = flamb
154			dEvdT = cvapor_exp
155			ELSE
156			lath = flamb+flami
157			dEvdT = cvapor_exp_ice
158			ENDIF
159	mlosch	1.1
160	jmc	1.3	Cph This statement cannot be a PARAMETER statement in the header,
161			Cph but must come here; it is not fortran77 standard
162	mlosch	1.1	aln = log(ht/zref)
163
164	jmc	1.3	C-- Use atmospheric state to compute surface fluxes.
165	mlosch	1.1
166	jmc	1.3	C-- Compute the turbulent surface fluxes.
167	mlosch	1.1
168	jmc	1.3	C Initial guess: z/l=0.0; hu=ht=hq=z
169			C Iterations: converge on z/l and hence the fluxes.
170			C t0 : virtual temperature (K)
171			C ssq : sea surface humidity (kg/kg)
172			C deltap : potential temperature diff (K)
173	mlosch	1.1
174			if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
175			t0 = atemp(i,j,bi,bj)*
176			& (exf_one + humid_fac*aqh(i,j,bi,bj))
177	jmc	1.3	c tmpbulk= exf_BulkqSat(Tsf)
178			c ssq = saltsat*tmpbulk/atmrho
179			tmpbulk = cvapor_fac_ice/exp(cvapor_exp_ice/Tsf)
180			ssq = tmpbulk/atmrho
181			deltap = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
182	mlosch	1.1	delq = aqh(i,j,bi,bj) - ssq
183			stable = exf_half + sign(exf_half, deltap)
184	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
185			CADJ STORE sh(i,j,bi,bj) = comlev1_exf_1, key = ikey_1
186			#endif
187	mlosch	1.1	tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
188			rdn = sqrt(tmpbulk)
189			ustar = rdn*sh(i,j,bi,bj)
190			tmpbulk= exf_BulkRhn(stable)
191	jmc	1.3	tstar = tmpbulk*deltap
192			qstar = cdalton*delq
193	mlosch	1.1
194			do iter = 1,niter_bulk
195
196	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
197			ikey_2 = iter
198			& + niter_bulk*(i-1)
199			& + niter_bulksNx(j-1)
200			& + niter_bulksNxsNy*act1
201			& + niter_bulksNxsNymax1act2
202			& + niter_bulksNxsNymax1max2*act3
203			& + niter_bulksNxsNymax1max2max3act4
204
205			CADJ STORE rdn = comlev1_exf_2, key = ikey_2
206			CADJ STORE ustar = comlev1_exf_2, key = ikey_2
207			CADJ STORE qstar = comlev1_exf_2, key = ikey_2
208			CADJ STORE tstar = comlev1_exf_2, key = ikey_2
209			CADJ STORE sh(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
210			#endif
211
212	mlosch	1.1	huol = czol*(tstar/t0 +
213			& qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
214			& ustar**2
215			huol = max(huol,zolmin)
216			stable = exf_half + sign(exf_half, huol)
217			htol = huol*ht/hu
218			hqol = huol*hq/hu
219	jmc	1.3
220			C Evaluate all stability functions assuming hq = ht.
221	mlosch	1.1	xsq = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
222			x = sqrt(xsq)
223			psimh = -psim_fachuolstable +
224			& (exf_one - stable)*
225			& (log((exf_one + x(exf_two + x))
226			& (exf_one + xsq)/8.) - exf_two*atan(x) +
227			& pi*exf_half)
228			xsq = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
229			psixh = -psim_fachtolstable + (exf_one - stable)*
230			& exf_two*log((exf_one + xsq)/exf_two)
231	jmc	1.3
232			C Shift wind speed using old coefficient
233	mlosch	1.1	ccc rd = rdn/(exf_one + rdn/karman*
234			ccc & (log(hu/zref) - psimh) )
235			rd = rdn/(exf_one - rdn/karman*psimh )
236			shn = sh(i,j,bi,bj)*rd/rdn
237			uzn = max(shn, umin)
238	jmc	1.3
239			C Update the transfer coefficients at 10 meters
240			C and neutral stability.
241
242	mlosch	1.1	tmpbulk= exf_BulkCdn(uzn)
243			rdn = sqrt(tmpbulk)
244	jmc	1.3
245			C Shift all coefficients to the measurement height
246			C and stability.
247	mlosch	1.1	c rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))
248			rd = rdn/(exf_one - rdn/karman*psimh)
249			tmpbulk= exf_BulkRhn(stable)
250	jmc	1.3	rh = tmpbulk/( exf_one +
251	mlosch	1.1	& tmpbulk/karman*(aln - psixh) )
252	jmc	1.3	re = cdalton/( exf_one +
253	mlosch	1.1	& cdalton/karman*(aln - psixh) )
254	jmc	1.3
255			C Update ustar, tstar, qstar using updated, shifted
256			C coefficients.
257	mlosch	1.1	ustar = rd*sh(i,j,bi,bj)
258	jmc	1.3	qstar = re*delq
259	mlosch	1.1	tstar = rh*deltap
260			enddo
261
262	jmc	1.3	tau = atmrhoustar*2
263			tau = tau*us(i,j,bi,bj)/sh(i,j,bi,bj)
264	mlosch	1.1
265	jmc	1.3	evapLoc = -tau*qstar/ustar
266	jmc	1.4	hlLocal = -lath*evapLoc
267			hsLocal = atmcptautstar/ustar
268	mlosch	1.1	#ifndef EXF_READ_EVAP
269			cdm evap(i,j,bi,bj) = tau*qstar/ustar
270			cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
271	jmc	1.4	C jmc: do not reset evap which contains evaporation over ice-free ocean fraction
272			c evap(i,j,bi,bj) = -recip_rhonil*evapLoc
273	mlosch	1.1	#endif
274
275			C--- surf.Temp derivative of turbulent Fluxes
276	jmc	1.3	dEvdT = (taure/ustar)ssq*dEvdT/Ts2
277	mlosch	1.1	dflhdT = -lath*dEvdT
278	jmc	1.3	dfshdT = -atmcptaurh/ustar
279	mlosch	1.1
280			C--- Upward long wave radiation
281			IF ( iceornot.EQ.0 ) THEN
282			emiss = ocean_emissivity
283			ELSEIF (iceornot.EQ.2) THEN
284			emiss = snow_emissivity
285			ELSE
286			emiss = ice_emissivity
287			ENDIF
288			flwup = emissstefanBoltzmannTs2*Ts2
289	jmc	1.3	dflwupdT = emissstefanBoltzmannTs2Tsf 4. _d 0
290
291	mlosch	1.1	C-- Total derivative with respect to surface temperature
292			df0dT = -dflwupdT+dfshdT+dflhdT
293	jmc	1.3
294	jmc	1.7	#ifdef ALLOW_DOWNWARD_RADIATION
295	mlosch	1.1	flwNet_dwn = lwdown(i,j,bi,bj) - flwup
296	jmc	1.7	#else
297			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
298			#endif
299	jmc	1.4	flxExceptSw = flwNet_dwn + hsLocal + hlLocal
300	mlosch	1.1
301			endif
302
303			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
304
305	jmc	1.7	#else /* ALLOW_ATM_TEMP */
306			STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'
307			#endif /* ALLOW_ATM_TEMP */
308	mlosch	1.1	#endif /* ALLOW_EXF */
309
310			RETURN
311			END