/[MITgcm]/MITgcm/pkg/thsice/thsice_calc_thickn.F
ViewVC logotype

Contents of /MITgcm/pkg/thsice/thsice_calc_thickn.F

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


Revision 1.13 - (show annotations) (download)
Mon Aug 27 20:55:28 2007 UTC (16 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59g
Changes since 1.12: +31 -15 lines
prevent (or reduce) lateral melting if sea-ice fraction is too close to iceMaskMin

1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_calc_thickn.F,v 1.12 2007/06/03 22:17:48 jmc Exp $
2 C $Name: $
3
4 #include "THSICE_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: THSICE_CALC_THICKN
8 C !INTERFACE:
9 SUBROUTINE THSICE_CALC_THICKN(
10 I bi, bj, siLo, siHi, sjLo, sjHi,
11 I iMin,iMax, jMin,jMax, dBugFlag,
12 I iceMask, tFrz, tOce, v2oc,
13 I snowP, prcAtm, sHeat, flxCnB,
14 U icFrac, hIce, hSnow, tSrf, qIc1, qIc2,
15 U frwAtm, fzMlOc, flx2oc,
16 O frw2oc, fsalt,
17 I myTime, myIter, myThid )
18 C !DESCRIPTION: \bv
19 C *==========================================================*
20 C | S/R THSICE_CALC_THICKN
21 C | o Calculate ice & snow thickness changes
22 C *==========================================================*
23 C \ev
24 C ADAPTED FROM:
25 C LANL CICE.v2.0.2
26 C-----------------------------------------------------------------------
27 C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
28 C.. See Bitz, C. M. and W. H. Lipscomb, 1999: "An energy-conserving
29 C.. thermodynamic sea ice model for climate study." J. Geophys.
30 C.. Res., 104, 15669 - 15677.
31 C.. Winton, M., 1999: "A reformulated three-layer sea ice model."
32 C.. Submitted to J. Atmos. Ocean. Technol.
33 C.. authors Elizabeth C. Hunke and William Lipscomb
34 C.. Fluid Dynamics Group, Los Alamos National Laboratory
35 C-----------------------------------------------------------------------
36 Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
37 C.. Compute temperature change using Winton model with 2 ice layers, of
38 C.. which only the top layer has a variable heat capacity.
39
40 C---------------------------------
41 C parameters that control the partitioning between lateral (ice area) and
42 C vertical (ice thickness) ice volume changes.
43 C a) surface melting and bottom melting (runtime parameter: fracEnMelt):
44 C frace is the fraction of available heat that is used for
45 C lateral melting (and 1-frace reduces the thickness ) when
46 C o hi < hThinIce & frac > lowIcFrac2 : frace=1 (lateral melting only)
47 C o hThinIce < hi < hThickIce & frac > lowIcFrac1 : frace=fracEnMelt
48 C o hi > hThickIce or frac < lowIcFrac1 : frace=0 (thinning only)
49 C b) ocean freezing (and ice forming):
50 C - conductive heat flux (below sea-ice) always increases thickness.
51 C - under sea-ice, freezing potential (x iceFraction) is used to increase ice
52 C thickness or ice fraction (lateral growth), according to:
53 C o hi < hThinIce : use freezing potential to grow ice vertically;
54 C o hThinIce < hi < hThickIce : use partition factor fracEnFreez for lateral growth
55 c and (1-fracEnFreez) to increase thickness.
56 C o hi > hThickIce : use all freezing potential to grow ice laterally
57 C (up to areaMax)
58 C - over open ocean, use freezing potential [x(1-iceFraction)] to grow ice laterally
59 C - lateral growth forms ice of the same or =hNewIceMax thickness, the less of the 2.
60 C - starts to form sea-ice over fraction iceMaskMin, as minimum ice-volume is reached
61 C---------------------------------
62 C !USES:
63 IMPLICIT NONE
64
65 C == Global variables ===
66 #include "EEPARAMS.h"
67 #include "THSICE_SIZE.h"
68 #include "THSICE_PARAMS.h"
69 #ifdef ALLOW_AUTODIFF_TAMC
70 # include "SIZE.h"
71 # include "tamc.h"
72 # include "tamc_keys.h"
73 #endif
74
75 C !INPUT/OUTPUT PARAMETERS:
76 C == Routine Arguments ==
77 C siLo,siHi :: size of input/output array: 1rst dim. lower,higher bounds
78 C sjLo,sjHi :: size of input/output array: 2nd dim. lower,higher bounds
79 C bi,bj :: tile indices
80 C iMin,iMax :: computation domain: 1rst index range
81 C jMin,jMax :: computation domain: 2nd index range
82 C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
83 C--- Input:
84 C iceMask :: sea-ice fractional mask [0-1]
85 C tFrz (Tf) :: sea-water freezing temperature [oC] (function of S)
86 C tOce (oceTs) :: surface level oceanic temperature [oC]
87 C v2oc (oceV2s) :: square of ocean surface-level velocity [m2/s2]
88 C snowP (snowPr) :: snow precipitation [kg/m2/s]
89 C prcAtm (=) :: total precip from the atmosphere [kg/m2/s]
90 C sHeat(sHeating):: surf heating flux left to melt snow or ice (= Atmos-conduction)
91 C flxCnB (=) :: heat flux conducted through the ice to bottom surface
92 C--- Modified (input&output):
93 C icFrac(iceFrac):: fraction of grid area covered in ice
94 C hIce (hi) :: ice height [m]
95 C hSnow (hs) :: snow height [m]
96 C tSrf (Tsf) :: surface (ice or snow) temperature
97 C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level
98 C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level
99 C frwAtm (evpAtm):: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
100 C fzMlOc (frzmlt):: ocean mixed-layer freezing/melting potential [W/m2]
101 C flx2oc (qleft) :: net heat flux to ocean [W/m2] (> 0 downward)
102 C--- Output
103 C frw2oc (fresh) :: Total fresh water flux to ocean [kg/m2/s] (> 0 downward)
104 C fsalt (=) :: salt flux to ocean [g/m2/s] (> 0 downward)
105 C--- Input:
106 C myTime :: current Time of simulation [s]
107 C myIter :: current Iteration number in simulation
108 C myThid :: my Thread Id number
109 INTEGER siLo, siHi, sjLo, sjHi
110 INTEGER bi,bj
111 INTEGER iMin, iMax
112 INTEGER jMin, jMax
113 LOGICAL dBugFlag
114 _RL iceMask(siLo:siHi,sjLo:sjHi)
115 _RL tFrz (siLo:siHi,sjLo:sjHi)
116 _RL tOce (siLo:siHi,sjLo:sjHi)
117 _RL v2oc (siLo:siHi,sjLo:sjHi)
118 _RL snowP (siLo:siHi,sjLo:sjHi)
119 _RL prcAtm (siLo:siHi,sjLo:sjHi)
120 _RL sHeat (siLo:siHi,sjLo:sjHi)
121 _RL flxCnB (siLo:siHi,sjLo:sjHi)
122 _RL icFrac (siLo:siHi,sjLo:sjHi)
123 _RL hIce (siLo:siHi,sjLo:sjHi)
124 _RL hSnow (siLo:siHi,sjLo:sjHi)
125 _RL tSrf (siLo:siHi,sjLo:sjHi)
126 _RL qIc1 (siLo:siHi,sjLo:sjHi)
127 _RL qIc2 (siLo:siHi,sjLo:sjHi)
128 _RL frwAtm (siLo:siHi,sjLo:sjHi)
129 _RL fzMlOc (siLo:siHi,sjLo:sjHi)
130 _RL flx2oc (siLo:siHi,sjLo:sjHi)
131 _RL frw2oc (siLo:siHi,sjLo:sjHi)
132 _RL fsalt (siLo:siHi,sjLo:sjHi)
133 _RL myTime
134 INTEGER myIter
135 INTEGER myThid
136 CEOP
137
138 #ifdef ALLOW_THSICE
139
140 C !LOCAL VARIABLES:
141 C--- local copy of input/output argument list variables (see description above)
142 _RL frzmlt
143 _RL Tf
144 _RL oceTs, oceV2s, snowPr
145 _RL sHeating
146 c _RL flxCnB
147 c _RL evpAtm
148 _RL iceFrac
149 _RL hi
150 _RL hs
151 _RL Tsf
152 _RL qicen(nlyr)
153 _RL qleft
154 _RL fresh
155 c _RL fsalt
156 C == Local Variables ==
157 INTEGER i,j,k ! loop indices
158 _RL rec_nlyr ! reciprocal of number of ice layers (real value)
159 C evap :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
160 C Fbot :: oceanic heat flux used to melt/form ice [W/m2]
161 _RL evap
162 _RL Fbot
163 _RL etop ! energy for top melting (J m-2)
164 _RL ebot ! energy for bottom melting (J m-2)
165 _RL etope ! energy (from top) for lateral melting (J m-2)
166 _RL ebote ! energy (from bottom) for lateral melting (J m-2)
167 _RL extend ! total energy for lateral melting (J m-2)
168 _RL hnew(nlyr) ! new ice layer thickness (m)
169 _RL hlyr ! individual ice layer thickness (m)
170 _RL dhi ! change in ice thickness
171 _RL dhs ! change in snow thickness
172 _RL rq ! rho * q for a layer
173 _RL rqh ! rho * q * h for a layer
174 _RL qbot ! enthalpy for new ice at bottom surf (J/kg)
175 _RL dt ! timestep
176 _RL esurp ! surplus energy from melting (J m-2)
177 _RL mwater0 ! fresh water mass gained/lost (kg/m^2)
178 _RL msalt0 ! salt gained/lost (kg/m^2)
179 _RL freshe ! fresh water gain from extension melting
180 _RL salte ! salt gained from extension melting
181
182 _RL ustar, cpchr
183 _RL chi, chs
184 _RL frace, rs, hq
185 C lowIcFrac1 :: ice-fraction lower limit above which partial (lowIcFrac1)
186 C lowIcFrac2 :: or full (lowIcFrac2) lateral melting is allowed.
187 _RL lowIcFrac1, lowIcFrac2
188
189 C- define grid-point location where to print debugging values
190 #include "THSICE_DEBUG.h"
191
192 1010 FORMAT(A,I3,3F8.3)
193 1020 FORMAT(A,1P4E11.3)
194
195 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
196
197 #ifdef ALLOW_AUTODIFF_TAMC
198 act1 = bi - myBxLo(myThid)
199 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
200 act2 = bj - myByLo(myThid)
201 max2 = myByHi(myThid) - myByLo(myThid) + 1
202 act3 = myThid - 1
203 max3 = nTx*nTy
204 act4 = ikey_dynamics - 1
205 #endif /* ALLOW_AUTODIFF_TAMC */
206
207 rec_nlyr = nlyr
208 rec_nlyr = 1. _d 0 / rec_nlyr
209 dt = thSIce_deltaT
210
211 C for now, use hard coded threshold (iceMaskMin +1.% and +10.%)
212 lowIcFrac1 = iceMaskMin*1.01 _d 0
213 lowIcFrac2 = iceMaskMin*1.10 _d 0
214
215 DO j = jMin, jMax
216 DO i = iMin, iMax
217 #ifdef ALLOW_AUTODIFF_TAMC
218 ikey_1 = i
219 & + sNx*(j-1)
220 & + sNx*sNy*act1
221 & + sNx*sNy*max1*act2
222 & + sNx*sNy*max1*max2*act3
223 & + sNx*sNy*max1*max2*max3*act4
224 #endif /* ALLOW_AUTODIFF_TAMC */
225 C--
226 #ifdef ALLOW_AUTODIFF_TAMC
227 CADJ STORE frwatm(i,j) = comlev1_thsice_1, key=ikey_1
228 CADJ STORE fzmloc(i,j) = comlev1_thsice_1, key=ikey_1
229 CADJ STORE hice(i,j) = comlev1_thsice_1, key=ikey_1
230 CADJ STORE hsnow(i,j) = comlev1_thsice_1, key=ikey_1
231 CADJ STORE icfrac(i,j) = comlev1_thsice_1, key=ikey_1
232 CADJ STORE qic1(i,j) = comlev1_thsice_1, key=ikey_1
233 CADJ STORE qic2(i,j) = comlev1_thsice_1, key=ikey_1
234 #endif
235
236 IF (iceMask(i,j).GT.0. _d 0) THEN
237 Tf = tFrz(i,j)
238 oceTs = tOce(i,j)
239 oceV2s = v2oc(i,j)
240 snowPr = snowP(i,j)
241 c prcAtm = prcAtm(i,j)
242 sHeating= sHeat(i,j)
243 c flxCnB = flxCnB(i,j)
244 iceFrac = icFrac(i,j)
245 hi = hIce(i,j)
246 hs = hSnow(i,j)
247 Tsf = tSrf(i,j)
248 qicen(1)= qIc1(i,j)
249 qicen(2)= qIc2(i,j)
250 c evpAtm = frwAtm(i,j)
251 frzmlt = fzMlOc(i,j)
252 qleft = flx2oc(i,j)
253 c fresh = frw2oc(i,j)
254 c fsalt = fsalt(i,j)
255 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
256 C initialize energies
257 esurp = 0. _d 0
258
259 evap = frwAtm(i,j)
260
261 C......................................................................
262 C.. Compute growth and/or melting at the top and bottom surfaces.......
263 C......................................................................
264
265 IF (frzmlt.GE. 0. _d 0) THEN
266 C !-----------------------------------------------------------------
267 C ! freezing conditions
268 C !-----------------------------------------------------------------
269 Fbot = frzmlt
270 IF ( iceFrac.LT.iceMaskMax ) THEN
271 IF (hi.GT.hThickIce) THEN
272 C if higher than hThickIce, use all frzmlt energy to grow extra ice
273 Fbot = 0. _d 0
274 ELSEIF (hi.GE.hThinIce) THEN
275 C between hThinIce & hThickIce, use partition factor fracEnFreez
276 Fbot = (1. _d 0 - fracEnFreez)*frzmlt
277 ENDIF
278 ENDIF
279 ELSE
280 C !-----------------------------------------------------------------
281 C ! melting conditions
282 C !-----------------------------------------------------------------
283 ustar = 5. _d -2 !for no currents
284 C frictional velocity between ice and water
285 ustar = SQRT(0.00536 _d 0*oceV2s)
286 ustar=max(5. _d -3,ustar)
287 cpchr =cpWater*rhosw*bMeltCoef
288 Fbot = cpchr*(Tf-oceTs)*ustar ! < 0
289 Fbot = max(Fbot,frzmlt) ! frzmlt < Fbot < 0
290 Fbot = min(Fbot,0. _d 0)
291 ENDIF
292
293 C mass of fresh water and salt initially present in ice
294 mwater0 = rhos*hs + rhoi*hi
295 msalt0 = rhoi*hi*saltIce
296
297 #ifdef ALLOW_DBUG_THSICE
298 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
299 & 'ThSI_CALC_TH: evpAtm, frzmlt, Fbot =', frwAtm(i,j),frzmlt,Fbot
300 #endif
301
302 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
303
304 C Compute energy available for melting/growth.
305
306 IF ( hi.GT.hThickIce .OR. fracEnMelt.EQ.0. _d 0 ) THEN
307 C above certain height (or when no ice fractionation), only melt from top
308 frace = 0. _d 0
309 ELSEIF (hi.LT.hThinIce) THEN
310 C below a certain height, all energy goes to changing ice extent
311 frace = 1. _d 0
312 ELSE
313 frace = fracEnMelt
314 ENDIF
315 C Reduce lateral melting when ice fraction is low : the purpose is to avoid
316 C disappearing of (up to hThinIce thick) sea-ice by over doing lateral melting
317 C (which would bring iceFrac below iceMaskMin).
318 IF ( iceFrac.LE.lowIcFrac1 ) THEN
319 frace = 0. _d 0
320 ELSEIF (iceFrac.LE.lowIcFrac2 ) THEN
321 frace = MIN( frace, fracEnMelt )
322 ENDIF
323
324 c IF (Tsf .EQ. 0. _d 0 .AND. sHeating.GT.0. _d 0) THEN
325 IF ( sHeating.GT.0. _d 0 ) THEN
326 etop = (1. _d 0-frace)*sHeating * dt
327 etope = frace*sHeating * dt
328 ELSE
329 etop = 0. _d 0
330 etope = 0. _d 0
331 C jmc: found few cases where Tsf=0 & sHeating < 0 : add this line to conserv energy:
332 esurp = sHeating * dt
333 ENDIF
334 C-- flux at the base of sea-ice:
335 C conduction H.flx= flxCnB (+ =down); oceanic turbulent H.flx= Fbot (+ =down).
336 C- ==> energy available(+ => melt)= (flxCnB-Fbot)*dt
337 c IF (frzmlt.LT.0. _d 0) THEN
338 c ebot = (1. _d 0-frace)*(flxCnB-Fbot) * dt
339 c ebote = frace*(flxCnB-Fbot) * dt
340 c ELSE
341 c ebot = (flxCnB-Fbot) * dt
342 c ebote = 0. _d 0
343 c ENDIF
344 C- original formulation(above): Loose energy when flxCnB < Fbot < 0
345 ebot = (flxCnB(i,j)-Fbot) * dt
346 IF (ebot.GT.0. _d 0) THEN
347 ebote = frace*ebot
348 ebot = ebot-ebote
349 ELSE
350 ebote = 0. _d 0
351 ENDIF
352 #ifdef ALLOW_DBUG_THSICE
353 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
354 & 'ThSI_CALC_TH: etop,etope,ebot,ebote=', etop,etope,ebot,ebote
355 #endif
356
357 C Initialize layer thicknesses.
358 C Make sure internal ice temperatures do not exceed Tmlt.
359 C If they do, then eliminate the layer. (Dont think this will happen
360 C for reasonable values of i0.)
361
362 hlyr = hi * rec_nlyr
363 DO k = 1, nlyr
364 hnew(k) = hlyr
365 ENDDO
366
367 C Top melt: snow, then ice.
368
369 IF (etop .GT. 0. _d 0) THEN
370 #ifdef ALLOW_AUTODIFF_TAMC
371 CADJ STORE etop = comlev1_thsice_1, key=ikey_1
372 #endif
373 IF (hs. gt. 0. _d 0) THEN
374 rq = rhos * qsnow
375 rqh = rq * hs
376 IF (etop .LT. rqh) THEN
377 hs = hs - etop/rq
378 etop = 0. _d 0
379 ELSE
380 etop = etop - rqh
381 hs = 0. _d 0
382 ENDIF
383 ENDIF
384
385 DO k = 1, nlyr
386 #ifdef ALLOW_AUTODIFF_TAMC
387 ikey_2 = k
388 & + nlyr*(i-1)
389 & + nlyr*sNx*(j-1)
390 & + nlyr*sNx*sNy*act1
391 & + nlyr*sNx*sNy*max1*act2
392 & + nlyr*sNx*sNy*max1*max2*act3
393 & + nlyr*sNx*sNy*max1*max2*max3*act4
394 #endif
395 C--
396 #ifdef ALLOW_AUTODIFF_TAMC
397 CADJ STORE etop = comlev1_thsice_2, key=ikey_2
398 CADJ STORE hnew(k) = comlev1_thsice_2, key=ikey_2
399 #endif
400 IF (etop .GT. 0. _d 0) THEN
401 rq = rhoi * qicen(k)
402 rqh = rq * hnew(k)
403 IF (etop .LT. rqh) THEN
404 hnew(k) = hnew(k) - etop / rq
405 etop = 0. _d 0
406 ELSE
407 etop = etop - rqh
408 hnew(k) = 0. _d 0
409 ENDIF
410 ENDIF
411 ENDDO
412 ELSE
413 etop=0. _d 0
414 ENDIF
415 C If ice is gone and melting energy remains
416 c IF (etop .GT. 0. _d 0) THEN
417 c WRITE (6,*) 'QQ All ice melts from top ', i,j
418 c hi=0. _d 0
419 c go to 200
420 c ENDIF
421
422
423 C Bottom melt/growth.
424
425 IF (ebot .LT. 0. _d 0) THEN
426 C Compute enthalpy of new ice growing at bottom surface.
427 qbot = -cpIce *Tf + Lfresh
428 dhi = -ebot / (qbot * rhoi)
429 ebot = 0. _d 0
430 cph k = nlyr
431 #ifdef ALLOW_AUTODIFF_TAMC
432 CADJ STORE hnew(:) = comlev1_thsice_1, key=ikey_1
433 CADJ STORE qicen(:) = comlev1_thsice_1, key=ikey_1
434 #endif
435 qicen(nlyr) = (hnew(nlyr)*qicen(nlyr)+dhi*qbot) /
436 & (hnew(nlyr)+dhi)
437 hnew(nlyr) = hnew(nlyr) + dhi
438 ELSE
439 DO k = nlyr, 1, -1
440 #ifdef ALLOW_AUTODIFF_TAMC
441 ikey_2 = (nlyr-k+1)
442 & + nlyr*(i-1)
443 & + nlyr*sNx*(j-1)
444 & + nlyr*sNx*sNy*act1
445 & + nlyr*sNx*sNy*max1*act2
446 & + nlyr*sNx*sNy*max1*max2*act3
447 & + nlyr*sNx*sNy*max1*max2*max3*act4
448 #endif
449 C--
450 #ifdef ALLOW_AUTODIFF_TAMC
451 CADJ STORE ebot = comlev1_thsice_2, key=ikey_2
452 CADJ STORE hnew(k) = comlev1_thsice_2, key=ikey_2
453 CADJ STORE qicen(k) = comlev1_thsice_2, key=ikey_2
454 #endif
455 IF (ebot.GT.0. _d 0 .AND. hnew(k).GT.0. _d 0) THEN
456 rq = rhoi * qicen(k)
457 rqh = rq * hnew(k)
458 IF (ebot .LT. rqh) THEN
459 hnew(k) = hnew(k) - ebot / rq
460 ebot = 0. _d 0
461 ELSE
462 ebot = ebot - rqh
463 hnew(k) = 0. _d 0
464 ENDIF
465 ENDIF
466 ENDDO
467
468 C If ice melts completely and snow is left, remove the snow with
469 C energy from the mixed layer
470
471 IF (ebot.GT.0. _d 0 .AND. hs.GT.0. _d 0) THEN
472 rq = rhos * qsnow
473 rqh = rq * hs
474 IF (ebot .LT. rqh) THEN
475 hs = hs - ebot / rq
476 ebot = 0. _d 0
477 ELSE
478 ebot = ebot - rqh
479 hs = 0. _d 0
480 ENDIF
481 ENDIF
482 c IF (ebot .GT. 0. _d 0) THEN
483 c IF (dBug) WRITE(6,*) 'All ice (& snow) melts from bottom'
484 c hi=0. _d 0
485 c go to 200
486 c ENDIF
487 ENDIF
488
489 C Compute new total ice thickness.
490 hi = 0. _d 0
491 DO k = 1, nlyr
492 hi = hi + hnew(k)
493 ENDDO
494 #ifdef ALLOW_DBUG_THSICE
495 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
496 & 'ThSI_CALC_TH: etop, ebot, hi, hs =', etop, ebot, hi, hs
497 #endif
498
499 C If hi < hIceMin, melt the ice.
500 IF ( hi.LT.hIceMin .AND. (hi+hs).GT.0. _d 0 ) THEN
501 esurp = esurp - rhos*qsnow*hs
502 DO k = 1, nlyr
503 esurp = esurp - rhoi*qicen(k)*hnew(k)
504 ENDDO
505 hi = 0. _d 0
506 hs = 0. _d 0
507 Tsf=0. _d 0
508 iceFrac = 0. _d 0
509 DO k = 1, nlyr
510 qicen(k) = 0. _d 0
511 ENDDO
512 #ifdef ALLOW_DBUG_THSICE
513 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
514 & 'ThSI_CALC_TH: -1 : esurp=',esurp
515 #endif
516 ENDIF
517
518 C-- do a mass-budget of sea-ice to compute "fresh" = the fresh-water flux
519 C that is returned to the ocean ; needs to be done before snow/evap
520 fresh = (mwater0 - (rhos*hs + rhoi*hi))/dt
521
522 IF ( hi .LE. 0. _d 0 ) THEN
523 C- return snow to the ocean (account for Latent heat of freezing)
524 fresh = fresh + snowPr
525 qleft = qleft - snowPr*Lfresh
526
527 ELSE
528 C- else: hi > 0
529
530 C Let it snow
531
532 hs = hs + dt*snowPr/rhos
533
534 C If ice evap is used to sublimate surface snow/ice or
535 C if no ice pass on to ocean
536 #ifdef ALLOW_AUTODIFF_TAMC
537 CADJ STORE evap = comlev1_thsice_1, key=ikey_1
538 CADJ STORE hs = comlev1_thsice_1, key=ikey_1
539 #endif
540 IF (hs.GT.0. _d 0) THEN
541 IF (evap/rhos *dt.GT.hs) THEN
542 evap=evap-hs*rhos/dt
543 hs=0. _d 0
544 ELSE
545 hs = hs - evap/rhos *dt
546 evap=0. _d 0
547 ENDIF
548 ENDIF
549 #ifdef ALLOW_AUTODIFF_TAMC
550 CADJ STORE evap = comlev1_thsice_1, key=ikey_1
551 CADJ STORE hi = comlev1_thsice_1, key=ikey_1
552 CADJ STORE hnew(:) = comlev1_thsice_1, key=ikey_1
553 CADJ STORE qicen(:) = comlev1_thsice_1, key=ikey_1
554 #endif
555 IF (hi.GT.0. _d 0.AND.evap.GT.0. _d 0) THEN
556 DO k = 1, nlyr
557 #ifdef ALLOW_AUTODIFF_TAMC
558 ikey_2 = k
559 & + nlyr*(i-1)
560 & + nlyr*sNx*(j-1)
561 & + nlyr*sNx*sNy*act1
562 & + nlyr*sNx*sNy*max1*act2
563 & + nlyr*sNx*sNy*max1*max2*act3
564 & + nlyr*sNx*sNy*max1*max2*max3*act4
565 #endif
566 C--
567 #ifdef ALLOW_AUTODIFF_TAMC
568 CADJ STORE evap = comlev1_thsice_2, key=ikey_2
569 CADJ STORE hnew(k) = comlev1_thsice_2, key=ikey_2
570 CADJ STORE qicen(k) = comlev1_thsice_2, key=ikey_2
571 #endif
572 IF (evap .GT. 0. _d 0) THEN
573 C-- original scheme, does not care about ice temp.
574 C- this can produce small error (< 1.W/m2) in the Energy budget
575 c IF (evap/rhoi *dt.GT.hnew(k)) THEN
576 c evap=evap-hnew(k)*rhoi/dt
577 c hnew(k)=0. _d 0
578 c ELSE
579 c hnew(k) = hnew(k) - evap/rhoi *dt
580 c evap=0. _d 0
581 c ENDIF
582 C-- modified scheme. taking into account Ice enthalpy
583 dhi = evap/rhoi*dt
584 IF (dhi.GE.hnew(k)) THEN
585 evap=evap-hnew(k)*rhoi/dt
586 esurp = esurp - hnew(k)*rhoi*(qicen(k)-Lfresh)
587 hnew(k)=0. _d 0
588 ELSE
589 #ifdef ALLOW_AUTODIFF_TAMC
590 CADJ STORE hnew(k) = comlev1_thsice_2, key=ikey_2
591 #endif
592 hq = hnew(k)*qicen(k)-dhi*Lfresh
593 hnew(k) = hnew(k) - dhi
594 qicen(k)=hq/hnew(k)
595 evap=0. _d 0
596 ENDIF
597 C-------
598 ENDIF
599 ENDDO
600 ENDIF
601 c IF (evap .GT. 0. _d 0) THEN
602 c WRITE (6,*) 'BB All ice sublimates', i,j
603 c hi=0. _d 0
604 c go to 200
605 c ENDIF
606
607 C Compute new total ice thickness.
608
609 hi = 0. _d 0
610 DO k = 1, nlyr
611 hi = hi + hnew(k)
612 ENDDO
613
614 C If hi < hIceMin, melt the ice.
615 IF ( hi.GT.0. _d 0 .AND. hi.LT.hIceMin ) THEN
616 fresh = fresh + (rhos*hs + rhoi*hi)/dt
617 esurp = esurp - rhos*qsnow*hs
618 DO k = 1, nlyr
619 esurp = esurp - rhoi*qicen(k)*hnew(k)
620 ENDDO
621 hi = 0. _d 0
622 hs = 0. _d 0
623 Tsf=0. _d 0
624 iceFrac = 0. _d 0
625 DO k = 1, nlyr
626 qicen(k) = 0. _d 0
627 ENDDO
628 #ifdef ALLOW_DBUG_THSICE
629 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
630 & 'ThSI_CALC_TH: -2 : esurp,fresh=', esurp, fresh
631 #endif
632 ENDIF
633
634 C- else hi > 0: end
635 ENDIF
636
637 IF ( hi .GT. 0. _d 0 ) THEN
638
639 C If there is enough snow to lower the ice/snow interface below
640 C freeboard, convert enough snow to ice to bring the interface back
641 C to sea-level. Adjust enthalpy of top ice layer accordingly.
642
643 #ifdef ALLOW_AUTODIFF_TAMC
644 CADJ STORE hi = comlev1_thsice_1, key=ikey_1
645 CADJ STORE hs = comlev1_thsice_1, key=ikey_1
646 CADJ STORE hnew(:) = comlev1_thsice_1, key=ikey_1
647 CADJ STORE qicen(:) = comlev1_thsice_1, key=ikey_1
648 #endif
649 IF ( hs .GT. hi*rhoiw/rhos ) THEN
650 cBB WRITE (6,*) 'Freeboard adjusts'
651 dhi = (hs * rhos - hi * rhoiw) / rhosw
652 dhs = dhi * rhoi / rhos
653 #ifdef ALLOW_AUTODIFF_TAMC
654 CADJ STORE hnew(:) = comlev1_thsice_1, key=ikey_1
655 #endif
656 rqh = rhoi*qicen(1)*hnew(1) + rhos*qsnow*dhs
657 hnew(1) = hnew(1) + dhi
658 qicen(1) = rqh / (rhoi*hnew(1))
659 hi = hi + dhi
660 hs = hs - dhs
661 ENDIF
662
663
664 C limit ice height
665 C- NOTE: this part does not conserve Energy ;
666 C but surplus of fresh water and salt are taken into account.
667 IF (hi.GT.hiMax) THEN
668 cBB print*,'BBerr, hi>hiMax',i,j,hi
669 chi=hi-hiMax
670 DO k=1,nlyr
671 hnew(k)=hnew(k)-chi/2. _d 0
672 ENDDO
673 fresh = fresh + chi*rhoi/dt
674 ENDIF
675 IF (hs.GT.hsMax) THEN
676 c print*,'BBerr, hs>hsMax',i,j,hs
677 chs=hs-hsMax
678 hs=hsMax
679 fresh = fresh + chs*rhos/dt
680 ENDIF
681
682 C Compute new total ice thickness.
683
684 hi = 0. _d 0
685 DO k = 1, nlyr
686 hi = hi + hnew(k)
687 ENDDO
688
689 #ifdef ALLOW_DBUG_THSICE
690 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
691 & 'ThSI_CALC_TH: b-Winton: hnew, qice =', hnew, qicen
692 #endif
693
694 hlyr = hi * rec_nlyr
695 CALL THSICE_RESHAPE_LAYERS(
696 U qicen,
697 I hlyr, hnew, myThid )
698
699 #ifdef ALLOW_DBUG_THSICE
700 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
701 & 'ThSI_CALC_TH: iceFrac,hi, qtot, hs =', iceFrac,hi,
702 & (qicen(1)+qicen(2))*0.5, hs
703 #endif
704
705 C- if hi > 0 : end
706 ENDIF
707 200 CONTINUE
708
709 C- Compute surplus energy left over from melting.
710
711 IF (hi.LE.0. _d 0) iceFrac=0. _d 0
712
713 C.. heat fluxes left over for ocean
714 qleft = qleft + (Fbot+(esurp+etop+ebot)/dt)
715 #ifdef ALLOW_DBUG_THSICE
716 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
717 & 'ThSI_CALC_TH: [esurp,etop+ebot]/dt =',esurp/dt,etop/dt,ebot/dt
718 #endif
719
720 C-- Evaporation left to the ocean :
721 fresh = fresh - evap
722 C- Correct Atmos. fluxes for this different latent heat:
723 C evap was computed over freezing surf.(Tsf<0), latent heat = Lvap+Lfresh
724 C but should be Lvap only for the fraction "evap" that is left to the ocean.
725 qleft = qleft + evap*Lfresh
726
727 C fresh and salt fluxes
728 c fresh = (mwater0 - (rhos*(hs) + rhoi*(hi)))/dt-evap
729 c fsalt = (msalt0 - rhoi*hi*saltIce)/35. _d 0/dt ! for same units as fresh
730 C note (jmc): fresh is computed from a sea-ice mass budget that already
731 C contains, at this point, snow & evaporation (of snow & ice)
732 C but are not meant to be part of ice/ocean fresh-water flux.
733 C fix: a) like below or b) by making the budget before snow/evap is added
734 c fresh = (mwater0 - (rhos*(hs) + rhoi*(hi)))/dt
735 c & + snow(i,j,bi,bj)*rhos - frwAtm
736 fsalt(i,j) = (msalt0 - rhoi*hi*saltIce)/dt
737
738 #ifdef ALLOW_DBUG_THSICE
739 IF (dBug(i,j,bi,bj) ) THEN
740 WRITE(6,1020)'ThSI_CALC_TH:dH2O,Ev[kg],fresh,fsalt',
741 & (mwater0-(rhos*hs+rhoi*hi))/dt,evap,fresh,fsalt(i,j)
742 WRITE(6,1020)'ThSI_CALC_TH: Qleft,Fbot,extend/dt =',
743 & Qleft,Fbot,(etope+ebote)/dt
744 ENDIF
745 #endif
746
747 C-- add remaining liquid Precip (rain+RunOff) directly to ocean:
748 fresh = fresh + (prcAtm(i,j)-snowPr)
749
750 C-- note: at this point, iceFrac has not been changed (unless reset to zero)
751 C and it can only be reduced by lateral melting in the following part:
752
753 C calculate extent changes
754 extend=etope+ebote
755 IF (iceFrac.GT.0. _d 0.AND.extend.GT.0. _d 0) THEN
756 rq = rhoi * 0.5 _d 0*(qicen(1)+qicen(2))
757 rs = rhos * qsnow
758 rqh = rq * hi + rs * hs
759 freshe=(rhos*hs+rhoi*hi)/dt
760 salte=(rhoi*hi*saltIce)/dt
761 IF ( extend.LT.rqh ) THEN
762 iceFrac=(1. _d 0-extend/rqh)*iceFrac
763 ENDIF
764 IF ( extend.LT.rqh .AND. iceFrac.GE.iceMaskMin ) THEN
765 fresh=fresh+extend/rqh*freshe
766 fsalt(i,j)=fsalt(i,j)+extend/rqh*salte
767 ELSE
768 iceFrac=0. _d 0
769 hi=0. _d 0
770 hs=0. _d 0
771 qleft=qleft+(extend-rqh)/dt
772 fresh=fresh+freshe
773 fsalt(i,j)=fsalt(i,j)+salte
774 ENDIF
775 ELSEIF (extend.GT.0. _d 0) THEN
776 qleft=qleft+extend/dt
777 ENDIF
778
779 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
780 C-- Update output variables :
781
782 C- Diagnostic of Atmos. fresh water flux (E-P) over sea-ice :
783 C substract precip from Evap (<- stored in frwAtm array)
784 frwAtm(i,j) = frwAtm(i,j) - prcAtm(i,j)
785
786 C- update Mixed-Layer Freezing potential heat flux by substracting the
787 C part which has already been accounted for (Fbot):
788 fzMlOc(i,j) = frzmlt - Fbot*iceMask(i,j)
789
790 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
791 #ifdef ALLOW_DBUG_THSICE
792 IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
793 & 'ThSI_CALC_TH: iceFrac,flx2oc,fsalt,frw2oc=',
794 & iceFrac, qleft, fsalt(i,j), fresh
795 #endif
796 #ifdef CHECK_ENERGY_CONSERV
797 CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 0,
798 I iceMask(i,j), iceFrac, hi, hs, qicen,
799 I qleft, fresh, fsalt,
800 I myTime, myIter, myThid )
801 #endif /* CHECK_ENERGY_CONSERV */
802 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
803 C-- Update Sea-Ice state output:
804 icFrac(i,j) = iceFrac
805 hIce(i,j) = hi
806 hSnow(i,j ) = hs
807 tSrf(i,j) = Tsf
808 qIc1(i,j) = qicen(1)
809 qIc2(i,j) = qicen(2)
810 C-- Update oceanic flux output:
811 flx2oc(i,j) = qleft
812 frw2oc(i,j) = fresh
813 c fsalt(i,j) = fsalt
814 ENDIF
815 ENDDO
816 ENDDO
817 #endif /* ALLOW_THSICE */
818
819 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
820
821 RETURN
822 END

  ViewVC Help
Powered by ViewVC 1.1.22