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

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

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

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