/[MITgcm]/MITgcm/pkg/seaice/seaice_solve4temp.F
ViewVC logotype

Contents of /MITgcm/pkg/seaice/seaice_solve4temp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.22 - (show annotations) (download)
Mon Jan 30 20:51:43 2012 UTC (12 years, 4 months ago) by jmc
Branch: MAIN
Changes since 1.21: +14 -12 lines 
more editing/cleaning with no effect on results: change sign of F_c, +=upward:
LEGACY defined: F_c +=up is exactly what B(i,j) was in earlier version;
LEGACY undef: consitent with all the F_* (I don't want to be impolite) fluxes
              which are positive upward (and are now documented).
1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_solve4temp.F,v 1.21 2012/01/30 03:38:20 jmc Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5 #ifdef ALLOW_EXF
6 # include "EXF_OPTIONS.h"
7 #endif
8
9 CBOP
10 C !ROUTINE: SEAICE_SOLVE4TEMP
11 C !INTERFACE:
12 SUBROUTINE SEAICE_SOLVE4TEMP(
13 I UG, HICE_ACTUAL, HSNOW_ACTUAL,
14 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
15 I F_lh_max,
16 #endif
17 U TSURF,
18 O F_ia, IcePenetSW,
19 O FWsublim,
20 I bi, bj, myTime, myIter, myThid )
21
22 C !DESCRIPTION: \bv
23 C *==========================================================*
24 C | SUBROUTINE SOLVE4TEMP
25 C | o Calculate ice growth rate, surface fluxes and
26 C | temperature of ice surface.
27 C | see Hibler, MWR, 108, 1943-1973, 1980
28 C *==========================================================*
29 C \ev
30
31 C !USES:
32 IMPLICIT NONE
33 C === Global variables ===
34 #include "SIZE.h"
35 #include "GRID.h"
36 #include "EEPARAMS.h"
37 #include "PARAMS.h"
38 #include "FFIELDS.h"
39 #include "SEAICE_SIZE.h"
40 #include "SEAICE_PARAMS.h"
41 #include "SEAICE.h"
42 #ifdef SEAICE_VARIABLE_FREEZING_POINT
43 #include "DYNVARS.h"
44 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
45 #ifdef ALLOW_EXF
46 # include "EXF_FIELDS.h"
47 #endif
48 #ifdef ALLOW_AUTODIFF_TAMC
49 # include "tamc.h"
50 #endif
51
52 C !INPUT PARAMETERS:
53 C UG :: atmospheric wind speed (m/s)
54 C HICE_ACTUAL :: actual ice thickness
55 C HSNOW_ACTUAL :: actual snow thickness
56 C TSURF :: surface temperature of ice/snow in Kelvin
57 C bi,bj :: tile indices
58 C myTime :: current time in simulation
59 C myIter :: iteration number in simulation
60 C myThid :: my Thread Id number
61 C !OUTPUT PARAMETERS:
62 C TSURF :: updated surface temperature of ice/snow in Kelvin
63 C F_ia :: upward seaice/snow surface heat flux to atmosphere (W/m^2)
64 C IcePenetSW :: short wave heat flux transmitted through ice (+=upward)
65 C FWsublim :: fresh water (mass) flux due to sublimation (+=up)(kg/m^2/s)
66 _RL UG (1:sNx,1:sNy)
67 _RL HICE_ACTUAL (1:sNx,1:sNy)
68 _RL HSNOW_ACTUAL(1:sNx,1:sNy)
69 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
70 _RL F_lh_max (1:sNx,1:sNy)
71 #endif
72 _RL TSURF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
73 _RL F_ia (1:sNx,1:sNy)
74 _RL IcePenetSW (1:sNx,1:sNy)
75 _RL FWsublim (1:sNx,1:sNy)
76 INTEGER bi, bj
77 _RL myTime
78 INTEGER myIter, myThid
79 CEOP
80
81 #if defined(ALLOW_ATM_TEMP) && defined(ALLOW_DOWNWARD_RADIATION)
82 C !LOCAL VARIABLES:
83 C === Local variables ===
84 C i, j :: Loop counters
85 C kSrf :: vertical index of surface layer
86 INTEGER i, j
87 #ifdef SEAICE_VARIABLE_FREEZING_POINT
88 INTEGER kSrf
89 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
90 INTEGER ITER
91 C TB :: ocean temperature in contact with ice (=seawater freezing point) (K)
92 _RL TB (1:sNx,1:sNy)
93 _RL D1, D1I
94 _RL D3(1:sNx,1:sNy)
95 _RL TMELT, XKI, XKS, HCUT, XIO
96 _RL SurfMeltTemp
97 C effConduct :: effective conductivity of combined ice and snow
98 _RL effConduct(1:sNx,1:sNy)
99 C lhSublim :: latent heat of sublimation (SEAICE_lhEvap + SEAICE_lhFusion)
100 _RL lhSublim
101 C t1,t2,t3,t4 :: powers of temperature
102 _RL t1, t2, t3, t4
103
104 C Constants to calculate Saturation Vapor Pressure
105 _RL dqhice_dTice
106 #ifdef SEAICE_SOLVE4TEMP_LEGACY
107 _RL C1, C2, C3, C4, C5, QS1
108 #else /* SEAICE_SOLVE4TEMP_LEGACY */
109 _RL dFiDTs1
110 _RL lnTEN
111 _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
112 C specific humidity at ice surface variables
113 _RL mm_pi,mm_log10pi
114 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
115
116 C F_c :: conductive heat flux through seaice+snow (+=upward)
117 C F_lh :: latent heat flux (sublimation) (+=upward)
118 #ifdef SEAICE_SOLVE4TEMP_LEGACY
119 C A1 :: part of atmos surface flux (+=downward) independent of tsurf
120 C A2 :: part of atmos surface flux (+=upward) which depends on tsurf
121 C A3 :: derivative of (A2-F_c) versus tsurf
122 _RL A2 (1:sNx,1:sNy)
123 _RL A3 (1:sNx,1:sNy)
124 _RL A1 (1:sNx,1:sNy)
125 #else /* SEAICE_SOLVE4TEMP_LEGACY */
126 C F_lwu :: upward long-wave surface heat flux (+=upward)
127 C F_sens :: sensible surface heat flux (+=upward)
128 _RL F_lwu (1:sNx,1:sNy)
129 _RL F_sens (1:sNx,1:sNy)
130 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
131 _RL F_lh (1:sNx,1:sNy)
132 _RL F_c (1:sNx,1:sNy)
133 _RL qhice (1:sNx,1:sNy)
134 _RL absorbedSW (1:sNx,1:sNy)
135 _RL penetSWFrac
136
137 C local copies of global variables
138 _RL tsurfLoc (1:sNx,1:sNy)
139 _RL atempLoc (1:sNx,1:sNy)
140 _RL lwdownLoc (1:sNx,1:sNy)
141 _RL ALB (1:sNx,1:sNy)
142 _RL ALB_ICE (1:sNx,1:sNy)
143 _RL ALB_SNOW (1:sNx,1:sNy)
144 C iceOrNot :: this is HICE_ACTUAL.GT.0.
145 LOGICAL iceOrNot(1:sNx,1:sNy)
146 #ifdef SEAICE_DEBUG
147 C F_io_net :: upward conductive heat flux through seaice+snow
148 C F_ia_net :: net heat flux divergence at the sea ice/snow surface:
149 C includes ice conductive fluxes and atmospheric fluxes (W/m^2)
150 _RL F_io_net (1:sNx,1:sNy)
151 _RL F_ia_net (1:sNx,1:sNy)
152 #endif /* SEAICE_DEBUG */
153
154 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
155
156 #ifdef ALLOW_AUTODIFF_TAMC
157 CADJ INIT comlev1_solve4temp = COMMON, sNx*sNy*NMAX_TICE
158 #endif /* ALLOW_AUTODIFF_TAMC */
159
160 #ifdef SEAICE_SOLVE4TEMP_LEGACY
161 C MAYKUT CONSTANTS FOR SAT. VAP. PRESSURE TEMP. POLYNOMIAL
162 C1= 2.7798202 _d -06
163 C2= -2.6913393 _d -03
164 C3= 0.97920849 _d +00
165 C4= -158.63779 _d +00
166 C5= 9653.1925 _d +00
167
168 QS1=0.622 _d +00/1013.0 _d +00
169 #else /* SEAICE_SOLVE4TEMP_LEGACY */
170 lnTEN = LOG(10.0 _d 0)
171 aa1 = 2663.5 _d 0
172 aa2 = 12.537 _d 0
173 bb1 = 0.622 _d 0
174 bb2 = 1.0 _d 0 - bb1
175 Ppascals = 100000. _d 0
176 C cc0 = TEN ** aa2
177 cc0 = EXP(aa2*lnTEN)
178 cc1 = cc0*aa1*bb1*Ppascals*lnTEN
179 cc2 = cc0*bb2
180 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
181
182 #ifdef SEAICE_VARIABLE_FREEZING_POINT
183 kSrf = 1
184 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
185
186 C SENSIBLE HEAT CONSTANT
187 D1=SEAICE_dalton*SEAICE_cpAir*SEAICE_rhoAir
188
189 C ICE LATENT HEAT CONSTANT
190 lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
191 D1I=SEAICE_dalton*lhSublim*SEAICE_rhoAir
192
193 C MELTING TEMPERATURE OF ICE
194 #ifdef SEAICE_SOLVE4TEMP_LEGACY
195 TMELT = 273.16 _d +00
196 SurfMeltTemp = 273.159 _d +00
197 #else /* SEAICE_SOLVE4TEMP_LEGACY */
198 TMELT = celsius2K
199 SurfMeltTemp = TMELT
200 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
201
202 C ICE CONDUCTIVITY
203 XKI=SEAICE_iceConduct
204
205 C SNOW CONDUCTIVITY
206 XKS=SEAICE_snowConduct
207
208 C CUTOFF SNOW THICKNESS
209 C Snow-Thickness above HCUT: SW optically thick snow (=> snow-albedo).
210 C Snow-Thickness below HCUT: linear transition to ice-albedo
211 #ifdef SEAICE_SOLVE4TEMP_LEGACY
212 HCUT = SEAICE_snowThick
213 #else /* SEAICE_SOLVE4TEMP_LEGACY */
214 HCUT = 0. _d 0
215 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
216
217 C PENETRATION SHORTWAVE RADIATION FACTOR
218 XIO=SEAICE_shortwave
219
220 C Initialize variables
221 DO J=1,sNy
222 DO I=1,sNx
223 C HICE_ACTUAL is modified in this routine, but at the same time
224 C used to decided where there is ice, therefore we save this information
225 C here in a separate array
226 iceOrNot (I,J) = HICE_ACTUAL(I,J) .GT. 0. _d 0
227 IcePenetSW(I,J) = 0. _d 0
228 absorbedSW(I,J) = 0. _d 0
229 qhice (I,J) = 0. _d 0
230 F_ia (I,J) = 0. _d 0
231 c F_io_net (I,J) = 0. _d 0
232 c F_ia_net (I,J) = 0. _d 0
233 F_lh (I,J) = 0. _d 0
234
235 C Reset the snow/ice surface to TMELT and bound the atmospheric temperature
236 #ifdef SEAICE_SOLVE4TEMP_LEGACY
237 A1(I,J) = 0.0 _d 0
238 A2(I,J) = 0.0 _d 0
239 A3(I,J) = 0.0 _d 0
240 tsurfLoc (I,J) = MIN(273.16 _d 0 + MAX_TICE, TSURF(I,J,bi,bj))
241 atempLoc (I,J) = MAX(273.16 _d 0 + MIN_ATEMP,ATEMP(I,J,bi,bj))
242 lwdownLoc(I,J) = MAX( MIN_LWDOWN, LWDOWN(I,J,bi,bj) )
243 #else /* SEAICE_SOLVE4TEMP_LEGACY */
244 F_lwu (I,J) = 0. _d 0
245 F_sens (I,J) = 0. _d 0
246 tsurfLoc (I,J) = TSURF(I,J,bi,bj)
247 atempLoc (I,J) = MAX( celsius2K+MIN_ATEMP, ATEMP(I,J,bi,bj) )
248 lwdownLoc(I,J) = LWDOWN(I,J,bi,bj)
249 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
250
251 C FREEZING TEMPERATURE OF SEAWATER
252 #ifdef SEAICE_VARIABLE_FREEZING_POINT
253 C Use a variable seawater freezing point
254 TB(I,J) = -0.0575 _d 0*salt(I,J,kSrf,bi,bj) + 0.0901 _d 0
255 & + celsius2K
256 #else
257 C Use a constant freezing temperature (SEAICE_VARIABLE_FREEZING_POINT undef)
258 #ifdef SEAICE_SOLVE4TEMP_LEGACY
259 TB(I,J) = 271.2 _d 0
260 #else /* SEAICE_SOLVE4TEMP_LEGACY */
261 TB(I,J) = celsius2K + SEAICE_freeze
262 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
263 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
264 IF(HSNOW_ACTUAL(I,J).GT.0.0) THEN
265 C Stefan-Boltzmann constant times emissivity
266 D3(I,J)=SEAICE_snow_emiss*SEAICE_boltzmann
267 #ifdef EXF_LWDOWN_WITH_EMISSIVITY
268 C This is now [(1-emiss)*lwdown - lwdown]
269 lwdownloc(I,J) = SEAICE_snow_emiss*lwdownloc(I,J)
270 #else /* use the old hard wired inconsistent value */
271 lwdownloc(I,J) = 0.97 _d 0*lwdownloc(I,J)
272 #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
273 ELSE
274 C Stefan-Boltzmann constant times emissivity
275 D3(I,J)=SEAICE_ice_emiss*SEAICE_boltzmann
276 #ifdef EXF_LWDOWN_WITH_EMISSIVITY
277 C This is now [(1-emiss)*lwdown - lwdown]
278 lwdownloc(I,J) = SEAICE_ice_emiss*lwdownloc(I,J)
279 #else /* use the old hard wired inconsistent value */
280 lwdownloc(I,J) = 0.97 _d 0*lwdownloc(I,J)
281 #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
282 ENDIF
283 ENDDO
284 ENDDO
285
286 DO J=1,sNy
287 DO I=1,sNx
288
289 C DECIDE ON ALBEDO
290 IF ( iceOrNot(I,J) ) THEN
291
292 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 ) THEN
293 IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
294 ALB_ICE (I,J) = SEAICE_wetIceAlb_south
295 ALB_SNOW(I,J) = SEAICE_wetSnowAlb_south
296 ELSE ! no surface melting
297 ALB_ICE (I,J) = SEAICE_dryIceAlb_south
298 ALB_SNOW(I,J) = SEAICE_drySnowAlb_south
299 ENDIF
300 ELSE !/ Northern Hemisphere
301 IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
302 ALB_ICE (I,J) = SEAICE_wetIceAlb
303 ALB_SNOW(I,J) = SEAICE_wetSnowAlb
304 ELSE ! no surface melting
305 ALB_ICE (I,J) = SEAICE_dryIceAlb
306 ALB_SNOW(I,J) = SEAICE_drySnowAlb
307 ENDIF
308 ENDIF !/ Albedo for snow and ice
309
310 C If actual snow thickness exceeds the cutoff thickness, use snow albedo
311 IF (HSNOW_ACTUAL(I,J) .GT. HCUT) THEN
312 ALB(I,J) = ALB_SNOW(I,J)
313 ELSEIF ( HCUT.LE.ZERO ) THEN
314 ALB(I,J) = ALB_ICE(I,J)
315 ELSE
316 C otherwise, use linear transition between ice and snow albedo
317 ALB(I,J) = MIN( ALB_ICE(I,J) + HSNOW_ACTUAL(I,J)/HCUT
318 & *(ALB_SNOW(I,J) -ALB_ICE(I,J))
319 & , ALB_SNOW(I,J) )
320 ENDIF
321
322 C Determine the fraction of shortwave radiative flux remaining
323 C at ocean interface after scattering through the snow and ice.
324 C If snow is present, no radiation penetrates through snow+ice
325 IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
326 penetSWFrac = 0.0 _d 0
327 ELSE
328 penetSWFrac = XIO*EXP(-1.5 _d 0 * HICE_ACTUAL(I,J))
329 ENDIF
330 C The shortwave radiative flux leaving ocean beneath ice (+=up).
331 IcePenetSW(I,J) = -(1.0 _d 0 - ALB(I,J))
332 & *penetSWFrac * SWDOWN(I,J,bi,bj)
333 C The shortwave radiative flux convergence in the seaice.
334 absorbedSW(I,J) = (1.0 _d 0 - ALB(I,J))
335 & *(1.0 _d 0 - penetSWFrac)* SWDOWN(I,J,bi,bj)
336
337 #ifdef SEAICE_SOLVE4TEMP_LEGACY
338 C Now determine fixed (relative to tsurf) forcing term in heat budget
339 A1(I,J) = absorbedSW(I,J) + lwdownLoc(I,J)
340 & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
341 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
342
343 C The effective conductivity of the two-layer
344 C snow/ice system.
345 #ifdef SEAICE_SOLVE4TEMP_LEGACY
346 effConduct(I,J)=
347 & XKS/(HSNOW_ACTUAL(I,J)/HICE_ACTUAL(I,J) +
348 & XKS/XKI)/HICE_ACTUAL(I,J)
349 #else /* SEAICE_SOLVE4TEMP_LEGACY */
350 C Set a mininum sea ice thickness of 5 cm to bound
351 C the magnitude of conductive heat fluxes.
352 Cif * now taken care of by SEAICE_hice_reg in seaice_growth
353 c hice_tmp = max(HICE_ACTUAL(I,J),5. _d -2)
354 effConduct(I,J) = XKI * XKS /
355 & (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
356 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
357
358 #ifdef SEAICE_DEBUG
359 IF ( (I .EQ. SEAICE_debugPointI) .AND.
360 & (J .EQ. SEAICE_debugPointJ) ) THEN
361 print '(A,i6)','-----------------------------------'
362 print '(A,i6)','ibi merged initialization ', myIter
363 print '(A,i6,4(1x,D24.15))',
364 & 'ibi iter, TSL, TS ',myIter,
365 & tsurfLoc(I,J), TSURF(I,J,bi,bj)
366 print '(A,i6,4(1x,D24.15))',
367 & 'ibi iter, TMELT ',myIter,TMELT
368 print '(A,i6,4(1x,D24.15))',
369 & 'ibi iter, HIA, EFKCON ',myIter,
370 & HICE_ACTUAL(I,J), effConduct(I,J)
371 print '(A,i6,4(1x,D24.15))',
372 & 'ibi iter, HSNOW ',myIter,
373 & HSNOW_ACTUAL(I,J), ALB(I,J)
374 print '(A,i6)','-----------------------------------'
375 print '(A,i6)','ibi energy balance iterat ', myIter
376 ENDIF
377 #endif /* SEAICE_DEBUG */
378
379 ENDIF !/* iceOrNot */
380 ENDDO !/* i */
381 ENDDO !/* j */
382
383 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
384 DO ITER=1,IMAX_TICE
385 DO J=1,sNy
386 DO I=1,sNx
387 #ifdef ALLOW_AUTODIFF_TAMC
388 iicekey = I + sNx*(J-1) + (ITER-1)*sNx*sNy
389 CADJ STORE tsurfloc(i,j) = comlev1_solve4temp,
390 CADJ & key = iicekey, byte = isbyte
391 #endif /* ALLOW_AUTODIFF_TAMC */
392
393 IF ( iceOrNot(I,J) ) THEN
394
395 t1 = tsurfLoc(I,J)
396 t2 = t1*t1
397 t3 = t2*t1
398 t4 = t2*t2
399
400 C Calculate the specific humidity in the BL above the snow/ice
401 #ifdef SEAICE_SOLVE4TEMP_LEGACY
402 C Use the Maykut polynomial
403 qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
404 dqhice_dTice = 0. _d 0
405
406 #else /* SEAICE_SOLVE4TEMP_LEGACY */
407 C Use an approximation which is more accurate at low temperatures
408
409 C log 10 of the sat vap pressure
410 mm_log10pi = -aa1 / t1 + aa2
411 C The saturation vapor pressure (SVP) in the surface
412 C boundary layer (BL) above the snow/ice.
413 C mm_pi = TEN **(mm_log10pi)
414 C The following form does the same, but is faster
415 mm_pi = EXP(mm_log10pi*lnTEN)
416 qhice(I,J) = bb1*mm_pi / (Ppascals - (1.0 _d 0 - bb1) *
417 & mm_pi)
418 C A constant for SVP derivative w.r.t TICE
419 C cc3t = TEN **(aa1 / t1)
420 C The following form does the same, but is faster
421 cc3t = EXP(aa1 / t1 * lnTEN)
422 C d(qh)/d(TICE)
423 dqhice_dTice = cc1*cc3t/((cc2-cc3t*Ppascals)**2 *t2)
424 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
425
426 C Calculate the flux terms based on the updated tsurfLoc
427 F_c(I,J) = effConduct(I,J)*(TB(I,J)-tsurfLoc(I,J))
428 F_lh(I,J) = D1I*UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
429 #ifdef SEAICE_SOLVE4TEMP_LEGACY
430 A2(I,J) = D1*UG(I,J)*t1+D1I*UG(I,J)*qhice(I,J)+D3(I,J)*t4
431 A3(I,J) = 4.0 _d 0*D3(I,J)*t3 + effConduct(I,J)+D1*UG(I,J)
432 & + D1I*UG(I,J)*dqhice_dTice
433 #else /* SEAICE_SOLVE4TEMP_LEGACY */
434 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
435 C if the latent heat flux implied by tsurfLoc exceeds
436 C F_lh_max, cap F_lh and decouple the flux magnitude from TICE
437 IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
438 F_lh(I,J) = F_lh_max(I,J)
439 dqhice_dTice = ZERO
440 ENDIF
441 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
442
443 C d(F_ia)/d(TICE)
444 dFiDTs1 = 4.0 _d 0*D3(I,J)*t3 + effConduct(I,J) + D1*UG(I,J)
445 & + D1I*UG(I,J)*dqhice_dTice
446
447 F_lwu(I,J) = t4 * D3(I,J)
448 F_sens(I,J)= D1 * UG(I,J) * (t1 - atempLoc(I,J))
449 F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
450 & -F_c(I,J) + F_sens(I,J) + F_lh(I,J)
451 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
452
453 #ifdef SEAICE_DEBUG
454 IF ( (I .EQ. SEAICE_debugPointI) .AND.
455 & (J .EQ. SEAICE_debugPointJ) ) THEN
456 print '(A,i6,4(1x,D24.15))',
457 & 'ice-iter qhICE, ', ITER,qhIce(I,J)
458 #ifdef SEAICE_SOLVE4TEMP_LEGACY
459 print '(A,i6,4(1x,D24.15))',
460 & 'ice-iter A1 A2 B ',ITER,A1(I,J),A2(I,J),F_c(I,J)
461 print '(A,i6,4(1x,D24.15))',
462 & 'ice-iter A3 (-A1+A2) ',ITER,A3(I,J),-A1(I,J)+A2(I,J)
463 #else /* SEAICE_SOLVE4TEMP_LEGACY */
464 print '(A,i6,4(1x,D24.15))',
465 & 'ice-iter dFiDTs1 F_ia ', ITER, dFiDTs1, F_ia(I,J)
466 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
467 ENDIF
468 #endif /* SEAICE_DEBUG */
469
470 C Update tsurfLoc
471 #ifdef SEAICE_SOLVE4TEMP_LEGACY
472 C update tsurf as solution of : Fc = A2 - A1 + A3 *delta.tsurf
473 tsurfLoc(I,J)=tsurfLoc(I,J)
474 & +(A1(I,J)-A2(I,J)+F_c(I,J))/A3(I,J)
475 tsurfLoc(I,J) = MAX( 273.16 _d 0+MIN_TICE, tsurfLoc(I,J) )
476
477 #else /* SEAICE_SOLVE4TEMP_LEGACY */
478 tsurfLoc(I,J) = tsurfLoc(I,J) - F_ia(I,J) / dFiDTs1
479
480 C If the search leads to tsurfLoc < 50 Kelvin, restart the search
481 C at tsurfLoc = TMELT. Note that one solution to the energy balance problem
482 C is an extremely low temperature - a temperature far below realistic values.
483 IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) THEN
484 tsurfLoc(I,J) = TMELT
485 ENDIF
486 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
487 tsurfLoc(I,J) = MIN( tsurfLoc(I,J), TMELT )
488
489 #ifdef SEAICE_DEBUG
490 IF ( (I .EQ. SEAICE_debugPointI) .AND.
491 & (J .EQ. SEAICE_debugPointJ) ) THEN
492 print '(A,i6,4(1x,D24.15))',
493 & 'ice-iter tsurfLc,|dif|', ITER,
494 & tsurfLoc(I,J),
495 & LOG10(ABS(tsurfLoc(I,J) - t1))
496 ENDIF
497 #endif /* SEAICE_DEBUG */
498
499 ENDIF !/* iceOrNot */
500 ENDDO !/* i */
501 ENDDO !/* j */
502 ENDDO !/* Iterations */
503 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
504
505 DO J=1,sNy
506 DO I=1,sNx
507 IF ( iceOrNot(I,J) ) THEN
508
509 C Save updated tsurf and finalize the flux terms
510 TSURF(I,J,bi,bj) = tsurfLoc(I,J)
511
512 #ifdef SEAICE_SOLVE4TEMP_LEGACY
513 F_ia(I,J)=-A1(I,J)+A2(I,J)
514 #else /* SEAICE_SOLVE4TEMP_LEGACY */
515 C Recalculate the fluxes based on the (possibly) adjusted TSURF
516 t1 = tsurfLoc(I,J)
517 t2 = t1*t1
518 t3 = t2*t1
519 t4 = t2*t2
520
521 C log 10 of the sat vap pressure
522 mm_log10pi = -aa1 / t1 + aa2
523 C saturation vapor pressure
524 C mm_pi = TEN **(mm_log10pi)
525 C The following form does the same, but is faster
526 mm_pi = EXP(mm_log10pi*lnTEN)
527 C over ice specific humidity
528 qhice(I,J) = bb1*mm_pi/(Ppascals- (1.0 _d 0 - bb1) * mm_pi)
529 F_lh(I,J) = D1I * UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
530 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
531 IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
532 F_lh(I,J) = F_lh_max(I,J)
533 ENDIF
534 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
535
536 F_c(I,J) = effConduct(I,J) * (TB(I,J) - t1)
537 F_lwu(I,J) = t4 * D3(I,J)
538 F_sens(I,J) = D1 * UG(I,J) * (t1 - atempLoc(I,J))
539
540 C The flux between the ice/snow surface and the atmosphere.
541 C (excludes upward conductive fluxes)
542 F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
543 & + F_sens(I,J) + F_lh(I,J)
544 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
545
546 #ifdef SEAICE_MODIFY_GROWTH_ADJ
547 Cgf no additional dependency through solver, snow, etc.
548 IF ( SEAICEadjMODE.GE.2 ) THEN
549 CALL ZERO_ADJ_1D( 1, TSURF(I,J,bi,bj), myThid)
550 t1 = TSURF(I,J,bi,bj)
551 t2 = t1*t1
552 t3 = t2*t1
553 t4 = t2*t2
554 qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
555
556 A1(I,J)=0.3 _d 0 *SWDOWN(I,J,bi,bj)+lwdownLoc(I,J)
557 & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
558 A2(I,J)= D1*UG(I,J)*t1+D1I*UG(I,J)*qhice(I,J)+D3(I,J)*t4
559
560 F_ia(I,J)=-A1(I,J)+A2(I,J)
561 IcePenetSW(I,J)= 0. _d 0
562 ENDIF
563 #endif /* SEAICE_MODIFY_GROWTH_ADJ */
564
565 C Fresh water flux (kg/m^2/s) from latent heat of sublimation.
566 C F_lh is positive upward (sea ice looses heat) and FWsublim
567 C is also positive upward (atmosphere gains freshwater)
568 FWsublim(I,J) = F_lh(I,J)/lhSublim
569
570 #ifdef SEAICE_DEBUG
571 C Caclulate the net ice-ocean and ice-atmosphere fluxes
572 IF (F_c(I,J) .GT. 0.0 _d 0) THEN
573 F_io_net(I,J) = F_c(I,J)
574 F_ia_net(I,J) = 0.0 _d 0
575 ELSE
576 F_io_net(I,J) = 0.0 _d 0
577 F_ia_net(I,J) = F_ia(I,J)
578 ENDIF !/* conductive fluxes up or down */
579
580 IF ( (I .EQ. SEAICE_debugPointI) .AND.
581 & (J .EQ. SEAICE_debugPointJ) ) THEN
582 print '(A)','----------------------------------------'
583 print '(A,i6)','ibi complete ', myIter
584 print '(A,4(1x,D24.15))',
585 & 'ibi T(SURF, surfLoc,atmos) ',
586 & TSURF(I,J,bi,bj), tsurfLoc(I,J),atempLoc(I,J)
587 print '(A,4(1x,D24.15))',
588 & 'ibi LWL ', lwdownLoc(I,J)
589 print '(A,4(1x,D24.15))',
590 & 'ibi QSW(Total, Penetrating)',
591 & SWDOWN(I,J,bi,bj), IcePenetSW(I,J)
592 print '(A,4(1x,D24.15))',
593 & 'ibi qh(ATM ICE) ',
594 & AQH(I,J,bi,bj),qhice(I,J)
595 #ifndef SEAICE_SOLVE4TEMP_LEGACY
596 print '(A,4(1x,D24.15))',
597 & 'ibi F(lwd,swi,lwu) ',
598 & -lwdownLoc(I,J), -absorbedSW(I,J), F_lwu(I,J)
599 print '(A,4(1x,D24.15))',
600 & 'ibi F(c,lh,sens) ',
601 & F_c(I,J), F_lh(I,J), F_sens(I,J)
602 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
603 IF (F_lh_max(I,J) .GT. ZERO) THEN
604 print '(A,4(1x,D24.15))',
605 & 'ibi F_lh_max, F_lh/lhmax) ',
606 & F_lh_max(I,J), F_lh(I,J)/ F_lh_max(I,J)
607 ELSE
608 print '(A,4(1x,D24.15))',
609 & 'ibi F_lh_max = ZERO! '
610 ENDIF
611 print '(A,4(1x,D24.15))',
612 & 'ibi FWsub, FWsubm*dT/rhoI ',
613 & FWsublim(I,J),
614 & FWsublim(I,J)*SEAICE_deltaTtherm/SEAICE_rhoICE
615 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
616 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
617 print '(A,4(1x,D24.15))',
618 & 'ibi F_ia, F_ia_net, F_c ',
619 #ifdef SEAICE_SOLVE4TEMP_LEGACY
620 & -A1(I,J)+A2(I,J), -A1(I,J)+A2(I,J)-F_c(I,J), F_c(I,J)
621 #else /* SEAICE_SOLVE4TEMP_LEGACY */
622 & F_ia(I,J), F_ia_net(I,J), F_c(I,J)
623 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
624 print '(A)','----------------------------------------'
625 ENDIF
626 #endif /* SEAICE_DEBUG */
627
628 ENDIF !/* iceOrNot */
629 ENDDO !/* i */
630 ENDDO !/* j */
631
632 #endif /* ALLOW_ATM_TEMP && ALLOW_DOWNWARD_RADIATION */
633 RETURN
634 END
  ViewVC Help
Powered by ViewVC 1.1.22