forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
thsice_get_exf.F
531 lines (498 loc) · 18.6 KB
/
thsice_get_exf.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
#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 bi, bj, it2,
I iMin,iMax, jMin,jMax,
I icFlag, hSnow1, tsfCel,
O flxExcSw, dFlxdT, evapLoc, dEvdT,
I myTime, myIter, 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"
# include "THSICE_SIZE.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C bi,bj :: tile indices
C it :: solv4temp iteration
C iMin,iMax :: computation domain: 1rst index range
C jMin,jMax :: computation domain: 2nd index range
C icFlag :: True= get fluxes at this location ; False= do nothing
C hSnow1 :: snow height [m]
C tsfCel :: surface (ice or snow) temperature (oC)
C flxExcSw :: net (downward) surface heat flux, except short-wave [W/m2]
C dFlxdT :: 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 myTime :: current Time of simulation [s]
C myIter :: current Iteration number in simulation
C myThid :: my Thread Id number
INTEGER bi, bj
INTEGER it2
INTEGER iMin, iMax
INTEGER jMin, jMax
_RL icFlag (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hSnow1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tsfCel (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxExcSw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dFlxdT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL evapLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dEvdT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_EXF
#ifdef ALLOW_ATM_TEMP
#ifdef ALLOW_DOWNWARD_RADIATION
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
_RL hlLocal
INTEGER iter
INTEGER i, j
_RL czol
_RL wsm ! limited wind speed [m/s] (> umin)
_RL t0 ! virtual temperature [K]
C copied from exf_bulkformulae:
C these need to be 2D-arrays for vectorizing code
C turbulent temperature scale [K]
_RL tstar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C turbulent humidity scale [kg/kg]
_RL qstar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C friction velocity [m/s]
_RL ustar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C neutral, zref (=10m) values of rd
_RL rdn (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL rd (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = sqrt(Cd) [-]
_RL rh (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = Ch / sqrt(Cd) [-]
_RL re (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = Ce / sqrt(Cd) [-]
C specific humidity difference [kg/kg]
_RL delq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL deltap(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef EXF_CALC_ATMRHO
C local atmospheric density [kg/m^3]
_RL atmrho_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif
C
_RL ssq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ren, rhn ! neutral, zref (=10m) values of re, rh
_RL usn, usm ! neutral, zref (=10m) wind-speed (+limited)
_RL stable ! = 1 if stable ; = 0 if unstable
C stability parameter at zwd [-] (=z/Monin-Obuklov length)
_RL huol
_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
_RL dfshdT
_RL 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
_RL qsat_exp
#ifdef ALLOW_DBUG_THSICE
LOGICAL dBugFlag
INTEGER stdUnit
#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 ==
C- Define grid-point location where to print debugging values
#include "THSICE_DEBUG.h"
#ifdef ALLOW_DBUG_THSICE
dBugFlag = debugLevel.GE.debLevC
stdUnit = standardMessageUnit
#endif
C-- Set surface parameters :
zwln = LOG(hu/zref)
ztln = LOG(ht/zref)
czol = hu*karman*gravity_mks
ren = cDalton
C more abbreviations
lath = flamb+flami
qsat_fac = cvapor_fac_ice
qsat_exp = cvapor_exp_ice
C initialisation of local arrays
DO j = 1-OLy,sNy+OLy
DO i = 1-OLx,sNx+OLx
tstar(i,j) = 0. _d 0
qstar(i,j) = 0. _d 0
ustar(i,j) = 0. _d 0
rdn(i,j) = 0. _d 0
rd(i,j) = 0. _d 0
rh(i,j) = 0. _d 0
re(i,j) = 0. _d 0
delq(i,j) = 0. _d 0
deltap(i,j) = 0. _d 0
ssq(i,j) = 0. _d 0
ENDDO
ENDDO
C
DO j=jMin,jMax
DO i=iMin,iMax
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) .AND. (icFlag(i,j).GT.0. _d 0) )
& WRITE(stdUnit,'(A,2I4,2I2,2F12.6)')
& 'ThSI_GET_EXF: i,j,atemp,lwd=',
& i,j,bi,bj, atemp(i,j,bi,bj),lwdown(i,j,bi,bj)
#endif
#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*(it2-1)
& + sNx*sNy*MaxTsf*act1
& + sNx*sNy*MaxTsf*max1*act2
& + sNx*sNy*MaxTsf*max1*max2*act3
& + sNx*sNy*MaxTsf*max1*max2*max3*act4
#endif
C-- Use atmospheric state to compute surface fluxes.
IF ( (icFlag(i,j).GT.0. _d 0) .AND.
& (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
IF ( hSnow1(i,j).GT.3. _d -1 ) THEN
emiss = snow_emissivity
ELSE
emiss = ice_emissivity
ENDIF
C copy a few variables to names used in bulkf_formula_lay
Tsf = tsfCel(i,j)+cen2kel
Ts2 = Tsf*Tsf
C wind speed
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sh(i,j,bi,bj) = comlev1_thsice_3, key = ikey_1
#endif
wsm = sh(i,j,bi,bj)
C-- air - surface difference of temperature & humidity
c tmpbulk= exf_BulkqSat(Tsf)
c ssq(i,j) = saltsat*tmpbulk/atmrho
tmpbulk = qsat_fac*EXP(-qsat_exp/Tsf)
#ifdef EXF_CALC_ATMRHO
atmrho_loc(i,j) = apressure(i,j,bi,bj) /
& (287.04 _d 0*atemp(i,j,bi,bj)
& *(1. _d 0 + humid_fac*aqh(i,j,bi,bj)))
ssq(i,j) = tmpbulk/atmrho_loc(i,j)
#else
ssq(i,j) = tmpbulk/atmrho
#endif
deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
delq(i,j) = aqh(i,j,bi,bj) - ssq(i,j)
C Do the part of the output variables that do not depend
C on the ice here to save a few re-computations
C This is not yet dEvdT, but just a cheap way to save a 2D-field
C for ssq and recomputing Ts2 lateron
dEvdT(i,j) = ssq(i,j)*qsat_exp/Ts2
flwup = emiss*stefanBoltzmann*Ts2*Ts2
dflwupdT = emiss*stefanBoltzmann*Ts2*Tsf * 4. _d 0
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
C-- This is not yet the total derivative with respect to surface temperature
dFlxdT(i,j) = -dflwupdT
C-- This is not yet the Net downward radiation excluding shortwave
flxExcSw(i,j) = flwNet_dwn
ENDIF
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
IF ( useStabilityFct_overIce ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
#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*(it2-1)
& + sNx*sNy*MaxTsf*act1
& + sNx*sNy*MaxTsf*max1*act2
& + sNx*sNy*MaxTsf*max1*max2*act3
& + sNx*sNy*MaxTsf*max1*max2*max3*act4
C--
CADJ STORE sh(i,j,bi,bj) = comlev1_thsice_3, key = ikey_1
#endif
IF ( (icFlag(i,j).GT.0. _d 0) .AND.
& (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
C-- Compute the turbulent surface fluxes (function of stability).
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(i,j))
c tmpbulk = exf_BulkCdn(sh(i,j,bi,bj))
wsm = sh(i,j,bi,bj)
tmpbulk = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
IF (tmpbulk.NE.0.) THEN
rdn(i,j) = SQRT(tmpbulk)
ELSE
rdn(i,j) = 0. _d 0
ENDIF
C-- initial guess for exchange other coefficients:
c rhn = exf_BulkRhn(stable)
rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
C-- calculate turbulent scales
ustar(i,j) = rdn(i,j)*wsm
tstar(i,j) = rhn*deltap(i,j)
qstar(i,j) = ren*delq(i,j)
ENDIF
ENDDO
ENDDO
C start iteration
DO iter = 1,niter_bulk
DO j=jMin,jMax
DO i=iMin,iMax
IF ( (icFlag(i,j).GT.0. _d 0) .AND.
& (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
ikey_2 = iter
& + niter_bulk*(i-1)
& + niter_bulk*sNx*(j-1)
& + niter_bulk*sNx*sNy*(it2-1)
& + niter_bulk*sNx*sNy*MaxTsf*act1
& + niter_bulk*sNx*sNy*MaxTsf*max1*act2
& + niter_bulk*sNx*sNy*MaxTsf*max1*max2*act3
& + niter_bulk*sNx*sNy*MaxTsf*max1*max2*max3*act4
CADJ STORE rdn(i,j) = comlev1_thsice_5, key = ikey_2
CADJ STORE ustar(i,j) = comlev1_thsice_5, key = ikey_2
CADJ STORE qstar(i,j) = comlev1_thsice_5, key = ikey_2
CADJ STORE tstar(i,j) = comlev1_thsice_5, key = ikey_2
CADJ STORE sh(i,j,bi,bj) = comlev1_thsice_5, key = ikey_2
#endif
t0 = atemp(i,j,bi,bj)*
& (exf_one + humid_fac*aqh(i,j,bi,bj))
huol = (tstar(i,j)/t0 +
& qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
& )*czol/(ustar(i,j)*ustar(i,j))
#ifdef ALLOW_BULK_LARGEYEAGER04
C- Large&Yeager_2004 code:
huol = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
#else
C- Large&Pond_1981 code (zolmin default = -100):
huol = MAX(huol,zolmin)
#endif /* ALLOW_BULK_LARGEYEAGER04 */
htol = huol*ht/hu
hqol = huol*hq/hu
stable = exf_half + SIGN(exf_half, huol)
C Evaluate all stability functions assuming hq = ht.
#ifdef ALLOW_BULK_LARGEYEAGER04
C- Large&Yeager_2004 code:
xsq = SQRT( ABS(exf_one - huol*16. _d 0) )
#else
C- Large&Pond_1981 code:
xsq = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)
#endif /* ALLOW_BULK_LARGEYEAGER04 */
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 )
#ifdef ALLOW_BULK_LARGEYEAGER04
C- Large&Yeager_2004 code:
xsq = SQRT( ABS(exf_one - htol*16. _d 0) )
#else
C- Large&Pond_1981 code:
xsq = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
#endif /* ALLOW_BULK_LARGEYEAGER04 */
psixh = -psim_fac*htol*stable
& + (exf_one-stable)
& *exf_two*LOG( exf_half*(exf_one+xsq) )
C Shift wind speed using old coefficient
#ifdef ALLOW_BULK_LARGEYEAGER04
C-- Large&Yeager04:
usn = wspeed(i,j,bi,bj)
& /( exf_one + rdn(i,j)*(zwln-psimh)/karman )
#else
C-- Large&Pond1981:
usn = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
#endif /* ALLOW_BULK_LARGEYEAGER04 */
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(i,j) = 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.
#ifdef ALLOW_BULK_LARGEYEAGER04
rd(i,j)= rdn(i,j)/( exf_one + rdn(i,j)*(zwln-psimh)/karman )
#else
rd(i,j)= rdn(i,j)/( exf_one - rdn(i,j)/karman*psimh )
#endif /* ALLOW_BULK_LARGEYEAGER04 */
rh(i,j)= rhn/( exf_one + rhn*(ztln-psixh)/karman )
re(i,j)= ren/( exf_one + ren*(ztln-psixh)/karman )
C Update ustar, tstar, qstar using updated, shifted coefficients.
ustar(i,j) = rd(i,j)*sh(i,j,bi,bj)
qstar(i,j) = re(i,j)*delq(i,j)
tstar(i,j) = rh(i,j)*deltap(i,j)
ENDIF
C end i/j-loops
ENDDO
ENDDO
C end iteration loop
ENDDO
DO j=jMin,jMax
DO i=iMin,iMax
IF ( (icFlag(i,j).GT.0. _d 0) .AND.
& (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
#ifdef EXF_CALC_ATMRHO
tau = atmrho_loc(i,j)*rd(i,j)*wspeed(i,j,bi,bj)
#else
tau = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
#endif
evapLoc(i,j) = -tau*qstar(i,j)
hlLocal = -lath*evapLoc(i,j)
hsLocal = atmcp*tau*tstar(i,j)
c ustress = tau*rd(i,j)*UwindSpeed
c vstress = tau*rd(i,j)*VwindSpeed
C--- surf.Temp derivative of turbulent Fluxes
C complete computation of dEvdT
dEvdT(i,j) = (tau*re(i,j))*dEvdT(i,j)
dflhdT = -lath*dEvdT(i,j)
dfshdT = -atmcp*tau*rh(i,j)
C-- Update total derivative with respect to surface temperature
dFlxdT(i,j) = dFlxdT(i,j) + dfshdT + dflhdT
C-- Update net downward radiation excluding shortwave
flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
ENDIF
ENDDO
ENDDO
ELSE
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Compute the turbulent surface fluxes using fixed transfert Coeffs
C with no stability dependence ( useStabilityFct_overIce = false )
DO j=jMin,jMax
DO i=iMin,iMax
IF ( (icFlag(i,j).GT.0. _d 0) .AND.
& (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
wsm = sh(i,j,bi,bj)
#ifdef EXF_CALC_ATMRHO
tau = atmrho_loc(i,j)*exf_iceCe*wsm
#else
tau = atmrho*exf_iceCe*wsm
#endif
evapLoc(i,j) = -tau*delq(i,j)
hlLocal = -lath*evapLoc(i,j)
#ifdef EXF_CALC_ATMRHO
hsLocal = atmcp*atmrho_loc(i,j)
& *exf_iceCh*wsm*deltap(i,j)
#else
hsLocal = atmcp*atmrho*exf_iceCh*wsm*deltap(i,j)
#endif
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,4F12.6)')
& 'ThSI_GET_EXF: wsm,hl,hs,Lw=',
& wsm,hlLocal,hsLocal,flxExcSw(i,j)
#endif
C--- surf.Temp derivative of turbulent Fluxes
C complete computation of dEvdT
dEvdT(i,j) = tau*dEvdT(i,j)
dflhdT = -lath*dEvdT(i,j)
#ifdef EXF_CALC_ATMRHO
dfshdT = -atmcp*atmrho_loc(i,j)*exf_iceCh*wsm
#else
dfshdT = -atmcp*atmrho*exf_iceCh*wsm
#endif
C-- Update total derivative with respect to surface temperature
dFlxdT(i,j) = dFlxdT(i,j) + dfshdT + dflhdT
C-- Update net downward radiation excluding shortwave
flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,4F12.6)')
& 'ThSI_GET_EXF: flx,dFlxdT,evap,dEvdT',
& flxExcSw(i,j), dFlxdT(i,j), evapLoc(i,j),dEvdT(i,j)
#endif
ENDIF
ENDDO
ENDDO
C endif useStabilityFct_overIce
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
DO j=jMin,jMax
DO i=iMin,iMax
IF ( (icFlag(i,j).GT.0. _d 0) .AND.
& (atemp(i,j,bi,bj).LE.0. _d 0) ) THEN
C-- in case atemp is zero:
flxExcSw(i,j) = 0. _d 0
dFlxdT (i,j) = 0. _d 0
evapLoc (i,j) = 0. _d 0
dEvdT (i,j) = 0. _d 0
ENDIF
ENDDO
ENDDO
#else /* ALLOW_DOWNWARD_RADIATION */
STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
#endif /* ALLOW_DOWNWARD_RADIATION */
#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