Line data Source code
1 : !=========================================================================
2 : !
3 : ! Update ice and snow internal temperatures and compute
4 : ! thermodynamic growth rates and atmospheric fluxes.
5 : !
6 : ! NOTE: The thermodynamic calculation is split in two for load balancing.
7 : ! First icepack_therm_vertical computes vertical growth rates and coupler
8 : ! fluxes. Then icepack_therm_itd does thermodynamic calculations not
9 : ! needed for coupling.
10 : !
11 : ! authors: William H. Lipscomb, LANL
12 : ! C. M. Bitz, UW
13 : ! Elizabeth C. Hunke, LANL
14 : !
15 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
16 : ! 2004: Block structure added by William Lipscomb
17 : ! 2006: Streamlined for efficiency by Elizabeth Hunke
18 : ! Converted to free source form (F90)
19 :
20 : module icepack_therm_vertical
21 :
22 : use icepack_kinds
23 : use icepack_parameters, only: c0, c1, p001, p5, puny
24 : use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin
25 : use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity
26 : use icepack_parameters, only: ktherm, heat_capacity, calc_Tsfc
27 : use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair
28 : use icepack_parameters, only: rfracmin, rfracmax, pndaspect, dpscale, frzpnd
29 : use icepack_parameters, only: phi_i_mushy
30 :
31 : use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso
32 : use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
33 : use icepack_tracers, only: n_aero, n_iso
34 :
35 : use icepack_therm_shared, only: ferrmax, l_brine
36 : use icepack_therm_shared, only: calculate_tin_from_qin, Tmin
37 : use icepack_therm_shared, only: hi_min
38 : use icepack_therm_bl99, only: temperature_changes
39 : use icepack_therm_0layer, only: zerolayer_temperature
40 : use icepack_therm_mushy, only: temperature_changes_salinity
41 :
42 : use icepack_warnings, only: warnstr, icepack_warnings_add
43 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
44 :
45 : use icepack_mushy_physics, only: icepack_mushy_temperature_mush
46 : use icepack_mushy_physics, only: liquidus_temperature_mush
47 : use icepack_mushy_physics, only: enthalpy_mush, enthalpy_of_melting
48 :
49 : use icepack_aerosol, only: update_aerosol
50 : use icepack_isotope, only: update_isotope
51 : use icepack_atmo, only: neutral_drag_coeffs, icepack_atm_boundary
52 : use icepack_age, only: increment_age
53 : use icepack_firstyear, only: update_FYarea
54 : use icepack_flux, only: set_sfcflux, merge_fluxes
55 : use icepack_meltpond_cesm, only: compute_ponds_cesm
56 : use icepack_meltpond_lvl, only: compute_ponds_lvl
57 : use icepack_meltpond_topo, only: compute_ponds_topo
58 :
59 : implicit none
60 :
61 : private
62 : public :: frzmlt_bottom_lateral, &
63 : thermo_vertical, &
64 : icepack_step_therm1
65 :
66 : !=======================================================================
67 :
68 : contains
69 :
70 : !=======================================================================
71 : !
72 : ! Driver for updating ice and snow internal temperatures and
73 : ! computing thermodynamic growth rates and atmospheric fluxes.
74 : !
75 : ! authors: William H. Lipscomb, LANL
76 : ! C. M. Bitz, UW
77 :
78 4360794 : subroutine thermo_vertical (nilyr, nslyr, &
79 : dt, aicen, &
80 : vicen, vsnon, &
81 4360794 : Tsf, zSin, &
82 4360794 : zqin, zqsn, &
83 : apond, hpond, &
84 : tr_pond_topo,&
85 : flw, potT, &
86 : Qa, rhoa, &
87 : fsnow, fpond, &
88 : fbot, Tbot, &
89 : Tsnice, sss, &
90 : lhcoef, shcoef, &
91 : fswsfc, fswint, &
92 4360794 : Sswabs, Iswabs, &
93 : fsurfn, fcondtopn, &
94 : fcondbotn, &
95 : fsensn, flatn, &
96 : flwoutn, evapn, &
97 : evapsn, evapin, &
98 : freshn, fsaltn, &
99 : fhocnn, meltt, &
100 : melts, meltb, &
101 : congel, snoice, &
102 : mlt_onset, frz_onset, &
103 : yday, dsnow, &
104 : prescribed_ice)
105 :
106 : integer (kind=int_kind), intent(in) :: &
107 : nilyr , & ! number of ice layers
108 : nslyr ! number of snow layers
109 :
110 : real (kind=dbl_kind), intent(in) :: &
111 : dt ! time step
112 :
113 : ! ice state variables
114 : real (kind=dbl_kind), intent(inout) :: &
115 : aicen , & ! concentration of ice
116 : vicen , & ! volume per unit area of ice (m)
117 : vsnon ! volume per unit area of snow (m)
118 :
119 : ! tracers
120 : real (kind=dbl_kind), intent(inout) :: &
121 : Tsf , & ! ice/snow top surface temp, same as Tsfcn (deg C)
122 : apond , & ! melt pond area fraction
123 : hpond ! melt pond depth (m)
124 : ! iage ! ice age (s)
125 :
126 : logical (kind=log_kind), intent(in) :: &
127 : tr_pond_topo ! if .true., use melt pond tracer
128 :
129 : logical (kind=log_kind), intent(in), optional :: &
130 : prescribed_ice ! if .true., use prescribed ice instead of computed
131 :
132 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
133 : zqsn , & ! snow layer enthalpy, zqsn < 0 (J m-3)
134 : zqin , & ! ice layer enthalpy, zqin < 0 (J m-3)
135 : zSin ! internal ice layer salinities
136 :
137 : ! input from atmosphere
138 : real (kind=dbl_kind), &
139 : intent(in) :: &
140 : flw , & ! incoming longwave radiation (W/m^2)
141 : potT , & ! air potential temperature (K)
142 : Qa , & ! specific humidity (kg/kg)
143 : rhoa , & ! air density (kg/m^3)
144 : fsnow , & ! snowfall rate (kg m-2 s-1)
145 : shcoef , & ! transfer coefficient for sensible heat
146 : lhcoef ! transfer coefficient for latent heat
147 :
148 : real (kind=dbl_kind), &
149 : intent(inout) :: &
150 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
151 : fswint , & ! SW absorbed in ice interior, below surface (W m-2)
152 : fpond ! fresh water flux to ponds (kg/m^2/s)
153 :
154 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
155 : Sswabs , & ! SW radiation absorbed in snow layers (W m-2)
156 : Iswabs ! SW radiation absorbed in ice layers (W m-2)
157 :
158 : ! input from ocean
159 : real (kind=dbl_kind), intent(in) :: &
160 : fbot , & ! ice-ocean heat flux at bottom surface (W/m^2)
161 : Tbot , & ! ice bottom surface temperature (deg C)
162 : sss ! ocean salinity
163 :
164 : ! coupler fluxes to atmosphere
165 : real (kind=dbl_kind), intent(out):: &
166 : flwoutn , & ! outgoing longwave radiation (W/m^2)
167 : evapn , & ! evaporative water flux (kg/m^2/s)
168 : evapsn , & ! evaporative water flux over snow (kg/m^2/s)
169 : evapin ! evaporative water flux over ice (kg/m^2/s)
170 :
171 : ! Note: these are intent out if calc_Tsfc = T, otherwise intent in
172 : real (kind=dbl_kind), intent(inout):: &
173 : fsensn , & ! sensible heat flux (W/m^2)
174 : flatn , & ! latent heat flux (W/m^2)
175 : fsurfn , & ! net flux to top surface, excluding fcondtopn
176 : fcondtopn, & ! downward cond flux at top surface (W m-2)
177 : fcondbotn ! downward cond flux at bottom surface (W m-2)
178 :
179 : ! coupler fluxes to ocean
180 : real (kind=dbl_kind), intent(out):: &
181 : freshn , & ! fresh water flux to ocean (kg/m^2/s)
182 : fsaltn , & ! salt flux to ocean (kg/m^2/s)
183 : fhocnn ! net heat flux to ocean (W/m^2)
184 :
185 : ! diagnostic fields
186 : real (kind=dbl_kind), &
187 : intent(inout):: &
188 : Tsnice , & ! snow ice interface temperature (deg C)
189 : meltt , & ! top ice melt (m/step-->cm/day)
190 : melts , & ! snow melt (m/step-->cm/day)
191 : meltb , & ! basal ice melt (m/step-->cm/day)
192 : congel , & ! basal ice growth (m/step-->cm/day)
193 : snoice , & ! snow-ice formation (m/step-->cm/day)
194 : dsnow , & ! change in snow thickness (m/step-->cm/day)
195 : mlt_onset, & ! day of year that sfc melting begins
196 : frz_onset ! day of year that freezing begins (congel or frazil)
197 :
198 : real (kind=dbl_kind), intent(in) :: &
199 : yday ! day of year
200 :
201 : ! local variables
202 :
203 : integer (kind=int_kind) :: &
204 : k ! ice layer index
205 :
206 : real (kind=dbl_kind) :: &
207 1486692 : dhi , & ! change in ice thickness
208 1486692 : dhs ! change in snow thickness
209 :
210 : ! 2D state variables (thickness, temperature)
211 :
212 : real (kind=dbl_kind) :: &
213 1486692 : hilyr , & ! ice layer thickness
214 1486692 : hslyr , & ! snow layer thickness
215 1486692 : hin , & ! ice thickness (m)
216 1486692 : hsn , & ! snow thickness (m)
217 1486692 : hsn_new , & ! thickness of new snow (m)
218 1486692 : worki , & ! local work array
219 1486692 : works ! local work array
220 :
221 : real (kind=dbl_kind), dimension (nilyr) :: &
222 19971408 : zTin ! internal ice layer temperatures
223 :
224 : real (kind=dbl_kind), dimension (nslyr) :: &
225 7656394 : zTsn ! internal snow layer temperatures
226 :
227 : ! other 2D flux and energy variables
228 :
229 : real (kind=dbl_kind) :: &
230 1486692 : einit , & ! initial energy of melting (J m-2)
231 1486692 : efinal , & ! final energy of melting (J m-2)
232 1486692 : einter ! intermediate energy
233 :
234 : real (kind=dbl_kind) :: &
235 1486692 : fadvocn ! advective heat flux to ocean
236 :
237 : character(len=*),parameter :: subname='(thermo_vertical)'
238 :
239 : !-----------------------------------------------------------------
240 : ! Initialize
241 : !-----------------------------------------------------------------
242 :
243 4360794 : flwoutn = c0
244 4360794 : evapn = c0
245 4360794 : evapsn = c0
246 4360794 : evapin = c0
247 4360794 : freshn = c0
248 4360794 : fsaltn = c0
249 4360794 : fhocnn = c0
250 4360794 : fadvocn = c0
251 :
252 4360794 : meltt = c0
253 4360794 : meltb = c0
254 4360794 : melts = c0
255 4360794 : congel = c0
256 4360794 : snoice = c0
257 4360794 : dsnow = c0
258 9630500 : zTsn(:) = c0
259 33066714 : zTin(:) = c0
260 :
261 4360794 : if (calc_Tsfc) then
262 4360794 : fsensn = c0
263 4360794 : flatn = c0
264 4360794 : fsurfn = c0
265 4360794 : fcondtopn = c0
266 : endif
267 :
268 : !-----------------------------------------------------------------
269 : ! Compute variables needed for vertical thermo calculation
270 : !-----------------------------------------------------------------
271 :
272 : call init_vertical_profile (nilyr, nslyr, &
273 : aicen, &
274 : vicen, vsnon, &
275 : hin, hilyr, &
276 : hsn, hslyr, &
277 0 : zqin, zTin, &
278 0 : zqsn, zTsn, &
279 0 : zSin, &
280 4360794 : einit )
281 4360794 : if (icepack_warnings_aborted(subname)) return
282 :
283 : ! Save initial ice and snow thickness (for fresh and fsalt)
284 4360794 : worki = hin
285 4360794 : works = hsn
286 :
287 : !-----------------------------------------------------------------
288 : ! Compute new surface temperature and internal ice and snow
289 : ! temperatures.
290 : !-----------------------------------------------------------------
291 :
292 4360794 : if (heat_capacity) then ! usual case
293 :
294 4057521 : if (ktherm == 2) then
295 :
296 : call temperature_changes_salinity(dt, &
297 : nilyr, nslyr, &
298 : rhoa, flw, &
299 : potT, Qa, &
300 : shcoef, lhcoef, &
301 : fswsfc, fswint, &
302 0 : Sswabs, Iswabs, &
303 : hilyr, hslyr, &
304 : apond, hpond, &
305 0 : zqin, zTin, &
306 0 : zqsn, zTsn, &
307 0 : zSin, &
308 : Tsf, Tbot, &
309 : sss, &
310 : fsensn, flatn, &
311 : flwoutn, fsurfn, &
312 : fcondtopn, fcondbotn, &
313 3172532 : fadvocn, snoice)
314 3172532 : if (icepack_warnings_aborted(subname)) return
315 :
316 : else ! ktherm
317 :
318 : call temperature_changes(dt, &
319 : nilyr, nslyr, &
320 : rhoa, flw, &
321 : potT, Qa, &
322 : shcoef, lhcoef, &
323 : fswsfc, fswint, &
324 0 : Sswabs, Iswabs, &
325 : hilyr, hslyr, &
326 0 : zqin, zTin, &
327 0 : zqsn, zTsn, &
328 0 : zSin, &
329 : Tsf, Tbot, &
330 : fsensn, flatn, &
331 : flwoutn, fsurfn, &
332 : fcondtopn, fcondbotn, &
333 884989 : einit )
334 884989 : if (icepack_warnings_aborted(subname)) return
335 :
336 : endif ! ktherm
337 :
338 : else
339 :
340 303273 : if (calc_Tsfc) then
341 :
342 : call zerolayer_temperature(nilyr, nslyr, &
343 : rhoa, flw, &
344 : potT, Qa, &
345 : shcoef, lhcoef, &
346 : fswsfc, &
347 : hilyr, hslyr, &
348 : Tsf, Tbot, &
349 : fsensn, flatn, &
350 : flwoutn, fsurfn, &
351 303273 : fcondtopn, fcondbotn )
352 303273 : if (icepack_warnings_aborted(subname)) return
353 :
354 : else
355 :
356 : !------------------------------------------------------------
357 : ! Set fcondbot = fcondtop for zero layer thermodynamics
358 : ! fcondtop is set in call to set_sfcflux in step_therm1
359 : !------------------------------------------------------------
360 :
361 0 : fcondbotn = fcondtopn ! zero layer
362 :
363 : endif ! calc_Tsfc
364 :
365 : endif ! heat_capacity
366 :
367 : ! intermediate energy for error check
368 :
369 4360794 : einter = c0
370 9630500 : do k = 1, nslyr
371 9630500 : einter = einter + hslyr * zqsn(k)
372 : enddo ! k
373 33066714 : do k = 1, nilyr
374 33066714 : einter = einter + hilyr * zqin(k)
375 : enddo ! k
376 :
377 4360794 : Tsnice = c0
378 4360794 : if ((hslyr+hilyr) > puny) then
379 4360794 : if (hslyr > puny) then
380 3827465 : Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr)
381 : else
382 533329 : Tsnice = Tsf
383 : endif
384 : endif
385 :
386 4360794 : if (icepack_warnings_aborted(subname)) return
387 :
388 : !-----------------------------------------------------------------
389 : ! Compute growth and/or melting at the top and bottom surfaces.
390 : ! Add new snowfall.
391 : ! Repartition ice into equal-thickness layers, conserving energy.
392 : !-----------------------------------------------------------------
393 :
394 : call thickness_changes(nilyr, nslyr, &
395 : dt, yday, &
396 : efinal, &
397 : hin, hilyr, &
398 : hsn, hslyr, &
399 0 : zqin, zqsn, &
400 : fbot, Tbot, &
401 : flatn, fsurfn, &
402 : fcondtopn, fcondbotn, &
403 : fsnow, hsn_new, &
404 : fhocnn, evapn, &
405 : evapsn, evapin, &
406 : meltt, melts, &
407 : meltb, &
408 : congel, snoice, &
409 : mlt_onset, frz_onset, &
410 0 : zSin, sss, &
411 4360794 : dsnow)
412 4360794 : if (icepack_warnings_aborted(subname)) return
413 :
414 : !-----------------------------------------------------------------
415 : ! Check for energy conservation by comparing the change in energy
416 : ! to the net energy input
417 : !-----------------------------------------------------------------
418 :
419 : call conservation_check_vthermo(dt, &
420 : fsurfn, flatn, &
421 : fhocnn, fswint, &
422 : fsnow, einit, &
423 : einter, efinal, &
424 : fcondtopn, fcondbotn, &
425 4360794 : fadvocn, fbot )
426 4360794 : if (icepack_warnings_aborted(subname)) return
427 :
428 : !-----------------------------------------------------------------
429 : ! If prescribed ice, set hi back to old values
430 : !-----------------------------------------------------------------
431 :
432 : #ifdef CESMCOUPLED
433 : if (present(prescribed_ice)) then
434 : if (prescribed_ice) then
435 : hin = worki
436 : fhocnn = c0 ! for diagnostics
437 : endif
438 : endif
439 : #endif
440 :
441 : !-----------------------------------------------------------------
442 : ! Compute fluxes of water and salt from ice to ocean.
443 : ! evapn < 0 => sublimation, evapn > 0 => condensation
444 : ! aerosol flux is accounted for in icepack_aerosol.F90
445 : !-----------------------------------------------------------------
446 :
447 4360794 : dhi = hin - worki
448 4360794 : dhs = hsn - works - hsn_new
449 :
450 4360794 : freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt
451 4360794 : fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
452 4360794 : fhocnn = fhocnn + fadvocn ! for ktherm=2
453 :
454 4360794 : if (hin == c0) then
455 7 : if (tr_pond_topo) fpond = fpond - aicen * apond * hpond
456 : endif
457 :
458 : !-----------------------------------------------------------------
459 : ! Given the vertical thermo state variables, compute the new ice
460 : ! state variables.
461 : !-----------------------------------------------------------------
462 :
463 : call update_state_vthermo(nilyr, nslyr, &
464 : Tbot, Tsf, &
465 : hin, hsn, &
466 0 : zqin, zSin, &
467 0 : zqsn, &
468 : aicen, &
469 4360794 : vicen, vsnon)
470 4360794 : if (icepack_warnings_aborted(subname)) return
471 :
472 : !-----------------------------------------------------------------
473 : ! Reload passive tracer array
474 : !-----------------------------------------------------------------
475 :
476 : end subroutine thermo_vertical
477 :
478 : !=======================================================================
479 : !
480 : ! Compute heat flux to bottom surface.
481 : ! Compute fraction of ice that melts laterally.
482 : !
483 : ! authors C. M. Bitz, UW
484 : ! William H. Lipscomb, LANL
485 : ! Elizabeth C. Hunke, LANL
486 :
487 1370160 : subroutine frzmlt_bottom_lateral (dt, ncat, &
488 : nilyr, nslyr, &
489 : aice, frzmlt, &
490 2740320 : vicen, vsnon, &
491 2740320 : qicen, qsnon, &
492 : sst, Tf, &
493 : ustar_min, &
494 474288 : fbot_xfer_type, &
495 : strocnxT, strocnyT, &
496 : Tbot, fbot, &
497 : rside, Cdn_ocn, &
498 : fside)
499 :
500 : integer (kind=int_kind), intent(in) :: &
501 : ncat , & ! number of thickness categories
502 : nilyr , & ! number of ice layers
503 : nslyr ! number of snow layers
504 :
505 : real (kind=dbl_kind), intent(in) :: &
506 : dt ! time step
507 :
508 : real (kind=dbl_kind), intent(in) :: &
509 : aice , & ! ice concentration
510 : frzmlt , & ! freezing/melting potential (W/m^2)
511 : sst , & ! sea surface temperature (C)
512 : Tf , & ! freezing temperature (C)
513 : ustar_min,& ! minimum friction velocity for ice-ocean heat flux
514 : Cdn_ocn , & ! ocean-ice neutral drag coefficient
515 : strocnxT, & ! ice-ocean stress, x-direction
516 : strocnyT ! ice-ocean stress, y-direction
517 :
518 : character (char_len), intent(in) :: &
519 : fbot_xfer_type ! transfer coefficient type for ice-ocean heat flux
520 :
521 : real (kind=dbl_kind), dimension(:), intent(in) :: &
522 : vicen , & ! ice volume (m)
523 : vsnon ! snow volume (m)
524 :
525 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
526 : qicen , & ! ice layer enthalpy (J m-3)
527 : qsnon ! snow layer enthalpy (J m-3)
528 :
529 : real (kind=dbl_kind), intent(out) :: &
530 : Tbot , & ! ice bottom surface temperature (deg C)
531 : fbot , & ! heat flux to ice bottom (W/m^2)
532 : rside , & ! fraction of ice that melts laterally
533 : fside ! lateral heat flux (W/m^2)
534 :
535 : ! local variables
536 :
537 : integer (kind=int_kind) :: &
538 : n , & ! thickness category index
539 : k ! layer index
540 :
541 : real (kind=dbl_kind) :: &
542 474288 : etot , & ! total energy in column
543 474288 : qavg ! average enthalpy in column (approximate)
544 :
545 : real (kind=dbl_kind) :: &
546 474288 : deltaT , & ! SST - Tbot >= 0
547 474288 : ustar , & ! skin friction velocity for fbot (m/s)
548 474288 : wlat , & ! lateral melt rate (m/s)
549 474288 : xtmp ! temporary variable
550 :
551 : ! Parameters for bottom melting
552 :
553 : real (kind=dbl_kind) :: &
554 474288 : cpchr ! -cp_ocn*rhow*exchange coefficient
555 :
556 : ! Parameters for lateral melting
557 :
558 : real (kind=dbl_kind), parameter :: &
559 : floediam = 300.0_dbl_kind, & ! effective floe diameter (m)
560 : floeshape = 0.66_dbl_kind , & ! constant from Steele (unitless)
561 : m1 = 1.6e-6_dbl_kind , & ! constant from Maykut & Perovich
562 : ! (m/s/deg^(-m2))
563 : m2 = 1.36_dbl_kind ! constant from Maykut & Perovich
564 : ! (unitless)
565 :
566 : character(len=*),parameter :: subname='(frzmlt_bottom_lateral)'
567 :
568 : !-----------------------------------------------------------------
569 : ! Identify grid cells where ice can melt.
570 : !-----------------------------------------------------------------
571 :
572 1370160 : rside = c0
573 1370160 : fside = c0
574 1370160 : Tbot = Tf
575 1370160 : fbot = c0
576 1370160 : wlat = c0
577 :
578 1370160 : if (aice > puny .and. frzmlt < c0) then ! ice can melt
579 :
580 : !-----------------------------------------------------------------
581 : ! Use boundary layer theory for fbot.
582 : ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703.
583 : !-----------------------------------------------------------------
584 :
585 499794 : deltaT = max((sst-Tbot),c0)
586 :
587 : ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
588 499794 : ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow)
589 499794 : ustar = max (ustar,ustar_min)
590 :
591 499794 : if (trim(fbot_xfer_type) == 'Cdn_ocn') then
592 : ! Note: Cdn_ocn has already been used for calculating ustar
593 : ! (formdrag only) --- David Schroeder (CPOM)
594 0 : cpchr = -cp_ocn*rhow*Cdn_ocn
595 : else ! fbot_xfer_type == 'constant'
596 : ! 0.006 = unitless param for basal heat flx ala McPhee and Maykut
597 499794 : cpchr = -cp_ocn*rhow*0.006_dbl_kind
598 : endif
599 :
600 499794 : fbot = cpchr * deltaT * ustar ! < 0
601 499794 : fbot = max (fbot, frzmlt) ! frzmlt < fbot < 0
602 :
603 : !!! uncomment to use all frzmlt for standalone runs
604 : ! fbot = min (c0, frzmlt)
605 :
606 : !-----------------------------------------------------------------
607 : ! Compute rside. See these references:
608 : ! Maykut and Perovich (1987): JGR, 92, 7032-7044
609 : ! Steele (1992): JGR, 97, 17,729-17,738
610 : !-----------------------------------------------------------------
611 :
612 499794 : wlat = m1 * deltaT**m2 ! Maykut & Perovich
613 499794 : rside = wlat*dt*pi/(floeshape*floediam) ! Steele
614 499794 : rside = max(c0,min(rside,c1))
615 :
616 : !-----------------------------------------------------------------
617 : ! Compute heat flux associated with this value of rside.
618 : !-----------------------------------------------------------------
619 :
620 2996636 : do n = 1, ncat
621 :
622 2496842 : etot = c0
623 2496842 : qavg = c0
624 :
625 : ! melting energy/unit area in each column, etot < 0
626 :
627 5958884 : do k = 1, nslyr
628 3462042 : etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
629 5958884 : qavg = qavg + qsnon(k,n)
630 : enddo
631 :
632 18688354 : do k = 1, nilyr
633 16191512 : etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind)
634 18688354 : qavg = qavg + qicen(k,n)
635 : enddo ! nilyr
636 :
637 : ! lateral heat flux, fside < 0
638 2996636 : if (tr_fsd) then ! floe size distribution
639 169535 : fside = fside + wlat*qavg
640 : else ! default floe size
641 2327307 : fside = fside + rside*etot/dt
642 : endif
643 :
644 : enddo ! n
645 :
646 : !-----------------------------------------------------------------
647 : ! Limit bottom and lateral heat fluxes if necessary.
648 : !-----------------------------------------------------------------
649 :
650 499794 : xtmp = frzmlt/(fbot + fside + puny)
651 499794 : xtmp = min(xtmp, c1)
652 499794 : fbot = fbot * xtmp
653 499794 : rside = rside * xtmp
654 499794 : fside = fside * xtmp
655 :
656 : endif
657 :
658 1370160 : end subroutine frzmlt_bottom_lateral
659 :
660 : !=======================================================================
661 : !
662 : ! Given the state variables (vicen, vsnon, zqin, etc.),
663 : ! compute variables needed for the vertical thermodynamics
664 : ! (hin, hsn, zTin, etc.)
665 : !
666 : ! authors William H. Lipscomb, LANL
667 : ! C. M. Bitz, UW
668 :
669 4360794 : subroutine init_vertical_profile(nilyr, nslyr, &
670 : aicen, vicen, &
671 : vsnon, &
672 : hin, hilyr, &
673 : hsn, hslyr, &
674 4360794 : zqin, zTin, &
675 4360794 : zqsn, zTsn, &
676 4360794 : zSin, &
677 : einit )
678 :
679 : integer (kind=int_kind), intent(in) :: &
680 : nilyr , & ! number of ice layers
681 : nslyr ! number of snow layers
682 :
683 : real (kind=dbl_kind), intent(in) :: &
684 : aicen , & ! concentration of ice
685 : vicen , & ! volume per unit area of ice (m)
686 : vsnon ! volume per unit area of snow (m)
687 :
688 : real (kind=dbl_kind), intent(out):: &
689 : hilyr , & ! ice layer thickness
690 : hslyr , & ! snow layer thickness
691 : einit ! initial energy of melting (J m-2)
692 :
693 : real (kind=dbl_kind), intent(out):: &
694 : hin , & ! ice thickness (m)
695 : hsn ! snow thickness (m)
696 :
697 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
698 : zqin , & ! ice layer enthalpy (J m-3)
699 : zTin ! internal ice layer temperatures
700 :
701 : real (kind=dbl_kind), dimension (:), intent(in) :: &
702 : zSin ! internal ice layer salinities
703 :
704 : real (kind=dbl_kind), dimension (:), &
705 : intent(inout) :: &
706 : zqsn , & ! snow enthalpy
707 : zTsn ! snow temperature
708 :
709 : ! local variables
710 : real (kind=dbl_kind), dimension(nilyr) :: &
711 15610614 : Tmlts ! melting temperature
712 :
713 : integer (kind=int_kind) :: &
714 : k ! ice layer index
715 :
716 : real (kind=dbl_kind) :: &
717 1486692 : rnslyr, & ! real(nslyr)
718 1486692 : Tmax ! maximum allowed snow/ice temperature (deg C)
719 :
720 : logical (kind=log_kind) :: & ! for vector-friendly error checks
721 : tsno_high , & ! flag for zTsn > Tmax
722 : tice_high , & ! flag for zTin > Tmlt
723 : tsno_low , & ! flag for zTsn < Tmin
724 : tice_low ! flag for zTin < Tmin
725 :
726 : character(len=*),parameter :: subname='(init_vertical_profile)'
727 :
728 : !-----------------------------------------------------------------
729 : ! Initialize
730 : !-----------------------------------------------------------------
731 :
732 4360794 : rnslyr = real(nslyr,kind=dbl_kind)
733 :
734 4360794 : tsno_high = .false.
735 4360794 : tice_high = .false.
736 4360794 : tsno_low = .false.
737 4360794 : tice_low = .false.
738 :
739 4360794 : einit = c0
740 :
741 : !-----------------------------------------------------------------
742 : ! Surface temperature, ice and snow thickness
743 : ! Initialize internal energy
744 : !-----------------------------------------------------------------
745 :
746 4360794 : hin = vicen / aicen
747 4360794 : hsn = vsnon / aicen
748 4360794 : hilyr = hin / real(nilyr,kind=dbl_kind)
749 4360794 : hslyr = hsn / rnslyr
750 :
751 : !-----------------------------------------------------------------
752 : ! Snow enthalpy and maximum allowed snow temperature
753 : ! If heat_capacity = F, zqsn and zTsn are used only for checking
754 : ! conservation.
755 : !-----------------------------------------------------------------
756 :
757 9630500 : do k = 1, nslyr
758 :
759 : !-----------------------------------------------------------------
760 : ! Tmax based on the idea that dT ~ dq / (rhos*cp_ice)
761 : ! dq ~ q dv / v
762 : ! dv ~ puny = eps11
763 : ! where 'd' denotes an error due to roundoff.
764 : !-----------------------------------------------------------------
765 :
766 5269706 : if (hslyr > hs_min/rnslyr .and. heat_capacity) then
767 : ! zqsn < 0
768 1055746 : Tmax = -zqsn(k)*puny*rnslyr / &
769 3169007 : (rhos*cp_ice*vsnon)
770 : else
771 2100699 : zqsn (k) = -rhos * Lfresh
772 2100699 : Tmax = puny
773 : endif
774 :
775 : !-----------------------------------------------------------------
776 : ! Compute snow temperatures from enthalpies.
777 : ! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0.
778 : !-----------------------------------------------------------------
779 5269706 : zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice
780 :
781 : !-----------------------------------------------------------------
782 : ! Check for zTsn > Tmax (allowing for roundoff error) and zTsn < Tmin.
783 : !-----------------------------------------------------------------
784 9630500 : if (zTsn(k) > Tmax) then
785 0 : tsno_high = .true.
786 5269706 : elseif (zTsn(k) < Tmin) then
787 0 : tsno_low = .true.
788 : endif
789 :
790 : enddo ! nslyr
791 :
792 : !-----------------------------------------------------------------
793 : ! If zTsn is out of bounds, print diagnostics and exit.
794 : !-----------------------------------------------------------------
795 :
796 4360794 : if (tsno_high .and. heat_capacity) then
797 0 : do k = 1, nslyr
798 :
799 0 : if (hslyr > hs_min/rnslyr) then
800 0 : Tmax = -zqsn(k)*puny*rnslyr / &
801 0 : (rhos*cp_ice*vsnon)
802 : else
803 0 : Tmax = puny
804 : endif
805 :
806 0 : if (zTsn(k) > Tmax) then
807 0 : write(warnstr,*) ' '
808 0 : call icepack_warnings_add(warnstr)
809 0 : write(warnstr,*) subname, 'Starting thermo, zTsn > Tmax'
810 0 : call icepack_warnings_add(warnstr)
811 0 : write(warnstr,*) subname, 'zTsn=',zTsn(k)
812 0 : call icepack_warnings_add(warnstr)
813 0 : write(warnstr,*) subname, 'Tmax=',Tmax
814 0 : call icepack_warnings_add(warnstr)
815 0 : write(warnstr,*) subname, 'zqsn',zqsn(k),-Lfresh*rhos,zqsn(k)+Lfresh*rhos
816 0 : call icepack_warnings_add(warnstr)
817 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
818 0 : call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTsn > Tmax" )
819 0 : return
820 : endif
821 :
822 : enddo ! nslyr
823 : endif ! tsno_high
824 :
825 4360794 : if (tsno_low .and. heat_capacity) then
826 0 : do k = 1, nslyr
827 :
828 0 : if (zTsn(k) < Tmin) then ! allowing for roundoff error
829 0 : write(warnstr,*) ' '
830 0 : call icepack_warnings_add(warnstr)
831 0 : write(warnstr,*) subname, 'Starting thermo, zTsn < Tmin'
832 0 : call icepack_warnings_add(warnstr)
833 0 : write(warnstr,*) subname, 'zTsn=', zTsn(k)
834 0 : call icepack_warnings_add(warnstr)
835 0 : write(warnstr,*) subname, 'Tmin=', Tmin
836 0 : call icepack_warnings_add(warnstr)
837 0 : write(warnstr,*) subname, 'zqsn', zqsn(k)
838 0 : call icepack_warnings_add(warnstr)
839 0 : write(warnstr,*) subname, hin
840 0 : call icepack_warnings_add(warnstr)
841 0 : write(warnstr,*) subname, hsn
842 0 : call icepack_warnings_add(warnstr)
843 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
844 0 : call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTsn < Tmin" )
845 0 : return
846 : endif
847 :
848 : enddo ! nslyr
849 : endif ! tsno_low
850 :
851 9630500 : do k = 1, nslyr
852 :
853 5269706 : if (zTsn(k) > c0) then ! correct roundoff error
854 9213 : zTsn(k) = c0
855 9213 : zqsn(k) = -rhos*Lfresh
856 : endif
857 :
858 : !-----------------------------------------------------------------
859 : ! initial energy per unit area of ice/snow, relative to 0 C
860 : !-----------------------------------------------------------------
861 9630500 : einit = einit + hslyr*zqsn(k)
862 :
863 : enddo ! nslyr
864 :
865 33066714 : do k = 1, nilyr
866 :
867 : !---------------------------------------------------------------------
868 : ! Use initial salinity profile for thin ice
869 : !---------------------------------------------------------------------
870 :
871 28705920 : if (ktherm == 1 .and. zSin(k) < min_salin-puny) then
872 0 : write(warnstr,*) ' '
873 0 : call icepack_warnings_add(warnstr)
874 0 : write(warnstr,*) subname, 'Starting zSin < min_salin, layer', k
875 0 : call icepack_warnings_add(warnstr)
876 0 : write(warnstr,*) subname, 'zSin =', zSin(k)
877 0 : call icepack_warnings_add(warnstr)
878 0 : write(warnstr,*) subname, 'min_salin =', min_salin
879 0 : call icepack_warnings_add(warnstr)
880 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
881 0 : call icepack_warnings_add(subname//" init_vertical_profile: Starting zSin < min_salin, layer" )
882 0 : return
883 : endif
884 :
885 28705920 : if (ktherm == 2) then
886 22207724 : Tmlts(k) = liquidus_temperature_mush(zSin(k))
887 : else
888 6498196 : Tmlts(k) = -zSin(k) * depressT
889 : endif
890 :
891 : !-----------------------------------------------------------------
892 : ! Compute ice enthalpy
893 : ! If heat_capacity = F, zqin and zTin are used only for checking
894 : ! conservation.
895 : !-----------------------------------------------------------------
896 :
897 : !-----------------------------------------------------------------
898 : ! Compute ice temperatures from enthalpies using quadratic formula
899 : !-----------------------------------------------------------------
900 :
901 28705920 : if (ktherm == 2) then
902 22207724 : zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k))
903 : else
904 6498196 : zTin(k) = calculate_Tin_from_qin(zqin(k),Tmlts(k))
905 : endif
906 :
907 28705920 : if (l_brine) then
908 28402647 : Tmax = Tmlts(k)
909 : else ! fresh ice
910 303273 : Tmax = -zqin(k)*puny/(rhos*cp_ice*vicen)
911 : endif
912 :
913 : !-----------------------------------------------------------------
914 : ! Check for zTin > Tmax and zTin < Tmin
915 : !-----------------------------------------------------------------
916 28705920 : if (zTin(k) > Tmax) then
917 0 : tice_high = .true.
918 28705920 : elseif (zTin(k) < Tmin) then
919 0 : tice_low = .true.
920 : endif
921 :
922 : !-----------------------------------------------------------------
923 : ! If zTin is out of bounds, print diagnostics and exit.
924 : !-----------------------------------------------------------------
925 :
926 28705920 : if (tice_high .and. heat_capacity) then
927 :
928 0 : if (l_brine) then
929 0 : Tmax = Tmlts(k)
930 : else ! fresh ice
931 0 : Tmax = -zqin(k)*puny/(rhos*cp_ice*vicen)
932 : endif
933 :
934 0 : if (zTin(k) > Tmax) then
935 0 : write(warnstr,*) ' '
936 0 : call icepack_warnings_add(warnstr)
937 0 : write(warnstr,*) subname, 'Starting thermo, T > Tmax, layer', k
938 0 : call icepack_warnings_add(warnstr)
939 0 : write(warnstr,*) subname, 'k:', k
940 0 : call icepack_warnings_add(warnstr)
941 0 : write(warnstr,*) subname, 'zTin =',zTin(k),', Tmax=',Tmax
942 0 : call icepack_warnings_add(warnstr)
943 0 : write(warnstr,*) subname, 'zSin =',zSin(k)
944 0 : call icepack_warnings_add(warnstr)
945 0 : write(warnstr,*) subname, 'hin =',hin
946 0 : call icepack_warnings_add(warnstr)
947 0 : write(warnstr,*) subname, 'zqin =',zqin(k)
948 0 : call icepack_warnings_add(warnstr)
949 0 : write(warnstr,*) subname, 'qmlt=',enthalpy_of_melting(zSin(k))
950 0 : call icepack_warnings_add(warnstr)
951 0 : write(warnstr,*) subname, 'Tmlt=',Tmlts(k)
952 0 : call icepack_warnings_add(warnstr)
953 :
954 0 : if (ktherm == 2) then
955 0 : zqin(k) = enthalpy_of_melting(zSin(k)) - c1
956 0 : zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k))
957 0 : write(warnstr,*) subname, 'Corrected quantities'
958 0 : call icepack_warnings_add(warnstr)
959 0 : write(warnstr,*) subname, 'zqin=',zqin(k)
960 0 : call icepack_warnings_add(warnstr)
961 0 : write(warnstr,*) subname, 'zTin=',zTin(k)
962 0 : call icepack_warnings_add(warnstr)
963 : else
964 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
965 0 : call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, T > Tmax, layer" )
966 0 : return
967 : endif
968 : endif
969 : endif ! tice_high
970 :
971 28705920 : if (tice_low .and. heat_capacity) then
972 :
973 0 : if (zTin(k) < Tmin) then
974 0 : write(warnstr,*) ' '
975 0 : call icepack_warnings_add(warnstr)
976 0 : write(warnstr,*) subname, 'Starting thermo T < Tmin, layer', k
977 0 : call icepack_warnings_add(warnstr)
978 0 : write(warnstr,*) subname, 'zTin =', zTin(k)
979 0 : call icepack_warnings_add(warnstr)
980 0 : write(warnstr,*) subname, 'Tmin =', Tmin
981 0 : call icepack_warnings_add(warnstr)
982 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
983 0 : call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, T < Tmin, layer" )
984 0 : return
985 : endif
986 : endif ! tice_low
987 :
988 : !-----------------------------------------------------------------
989 : ! correct roundoff error
990 : !-----------------------------------------------------------------
991 :
992 28705920 : if (ktherm /= 2) then
993 :
994 6498196 : if (zTin(k) > c0) then
995 43773 : zTin(k) = c0
996 43773 : zqin(k) = -rhoi*Lfresh
997 : endif
998 :
999 : endif
1000 :
1001 : ! echmod: is this necessary?
1002 : ! if (ktherm == 1) then
1003 : ! if (zTin(k)>= -zSin(k)*depressT) then
1004 : ! zTin(k) = -zSin(k)*depressT - puny
1005 : ! zqin(k) = -rhoi*cp_ocn*zSin(k)*depressT
1006 : ! endif
1007 : ! endif
1008 :
1009 : !-----------------------------------------------------------------
1010 : ! initial energy per unit area of ice/snow, relative to 0 C
1011 : !-----------------------------------------------------------------
1012 :
1013 33066714 : einit = einit + hilyr*zqin(k)
1014 :
1015 : enddo ! nilyr
1016 :
1017 : end subroutine init_vertical_profile
1018 :
1019 : !=======================================================================
1020 : !
1021 : ! Compute growth and/or melting at the top and bottom surfaces.
1022 : ! Convert snow to ice if necessary.
1023 : !
1024 : ! authors William H. Lipscomb, LANL
1025 : ! C. M. Bitz, UW
1026 :
1027 4360794 : subroutine thickness_changes (nilyr, nslyr, &
1028 : dt, yday, &
1029 : efinal, &
1030 : hin, hilyr, &
1031 : hsn, hslyr, &
1032 4360794 : zqin, zqsn, &
1033 : fbot, Tbot, &
1034 : flatn, fsurfn, &
1035 : fcondtopn, fcondbotn, &
1036 : fsnow, hsn_new, &
1037 : fhocnn, evapn, &
1038 : evapsn, evapin, &
1039 : meltt, melts, &
1040 : meltb, &
1041 : congel, snoice, &
1042 : mlt_onset, frz_onset,&
1043 4360794 : zSin, sss, &
1044 : dsnow)
1045 :
1046 : integer (kind=int_kind), intent(in) :: &
1047 : nilyr , & ! number of ice layers
1048 : nslyr ! number of snow layers
1049 :
1050 : real (kind=dbl_kind), intent(in) :: &
1051 : dt , & ! time step
1052 : yday ! day of the year
1053 :
1054 : real (kind=dbl_kind), intent(in) :: &
1055 : fbot , & ! ice-ocean heat flux at bottom surface (W/m^2)
1056 : Tbot , & ! ice bottom surface temperature (deg C)
1057 : fsnow , & ! snowfall rate (kg m-2 s-1)
1058 : flatn , & ! surface downward latent heat (W m-2)
1059 : fsurfn , & ! net flux to top surface, excluding fcondtopn
1060 : fcondtopn ! downward cond flux at top surface (W m-2)
1061 :
1062 : real (kind=dbl_kind), intent(inout) :: &
1063 : fcondbotn ! downward cond flux at bottom surface (W m-2)
1064 :
1065 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1066 : zqin , & ! ice layer enthalpy (J m-3)
1067 : zqsn ! snow layer enthalpy (J m-3)
1068 :
1069 : real (kind=dbl_kind), intent(inout) :: &
1070 : hilyr , & ! ice layer thickness (m)
1071 : hslyr ! snow layer thickness (m)
1072 :
1073 : real (kind=dbl_kind), intent(inout) :: &
1074 : meltt , & ! top ice melt (m/step-->cm/day)
1075 : melts , & ! snow melt (m/step-->cm/day)
1076 : meltb , & ! basal ice melt (m/step-->cm/day)
1077 : congel , & ! basal ice growth (m/step-->cm/day)
1078 : snoice , & ! snow-ice formation (m/step-->cm/day)
1079 : dsnow , & ! snow formation (m/step-->cm/day)
1080 : ! iage , & ! ice age (s)
1081 : mlt_onset , & ! day of year that sfc melting begins
1082 : frz_onset ! day of year that freezing begins (congel or frazil)
1083 :
1084 : real (kind=dbl_kind), intent(inout) :: &
1085 : hin , & ! total ice thickness (m)
1086 : hsn ! total snow thickness (m)
1087 :
1088 : real (kind=dbl_kind), intent(out):: &
1089 : efinal ! final energy of melting (J m-2)
1090 :
1091 : real (kind=dbl_kind), intent(out):: &
1092 : fhocnn , & ! fbot, corrected for any surplus energy (W m-2)
1093 : evapn , & ! ice/snow mass sublimated/condensed (kg m-2 s-1)
1094 : evapsn , & ! ice/snow mass sublimated/condensed over snow (kg m-2 s-1)
1095 : evapin ! ice/snow mass sublimated/condensed over ice (kg m-2 s-1)
1096 :
1097 : real (kind=dbl_kind), intent(out):: &
1098 : hsn_new ! thickness of new snow (m)
1099 :
1100 : ! changes to zSin in this subroutine are not reloaded into the
1101 : ! trcrn array for ktherm /= 2, so we could remove ktherm=2 conditionals
1102 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1103 : zSin ! ice layer salinity (ppt)
1104 :
1105 : real (kind=dbl_kind), intent(in) :: &
1106 : sss ! ocean salinity (PSU)
1107 :
1108 : ! local variables
1109 :
1110 : integer (kind=int_kind) :: &
1111 : k ! vertical index
1112 :
1113 : real (kind=dbl_kind) :: &
1114 1486692 : esub , & ! energy for sublimation, > 0 (J m-2)
1115 1486692 : econ , & ! energy for condensation, < 0 (J m-2)
1116 1486692 : etop_mlt , & ! energy for top melting, > 0 (J m-2)
1117 1486692 : ebot_mlt , & ! energy for bottom melting, > 0 (J m-2)
1118 1486692 : ebot_gro , & ! energy for bottom growth, < 0 (J m-2)
1119 1486692 : emlt_atm , & ! total energy of brine, from atmosphere (J m-2)
1120 1486692 : emlt_ocn ! total energy of brine, to ocean (J m-2)
1121 :
1122 : real (kind=dbl_kind) :: &
1123 1486692 : qbotmax , & ! max enthalpy of ice growing at bottom
1124 1486692 : dhi , & ! change in ice thickness
1125 1486692 : dhs , & ! change in snow thickness
1126 1486692 : Ti , & ! ice temperature
1127 1486692 : Ts , & ! snow temperature
1128 1486692 : qbot , & ! enthalpy of ice growing at bottom surface (J m-3)
1129 1486692 : qsub , & ! energy/unit volume to sublimate ice/snow (J m-3)
1130 1486692 : hqtot , & ! sum of h*q for two layers
1131 1486692 : wk1 , & ! temporary variable
1132 1486692 : zqsnew , & ! enthalpy of new snow (J m-3)
1133 1486692 : hstot , & ! snow thickness including new snow (m)
1134 1486692 : Tmlts ! melting temperature
1135 :
1136 : real (kind=dbl_kind), dimension (nilyr+1) :: &
1137 18484716 : zi1 , & ! depth of ice layer boundaries (m)
1138 18484716 : zi2 ! adjusted depths, with equal hilyr (m)
1139 :
1140 : real (kind=dbl_kind), dimension (nslyr+1) :: &
1141 10530496 : zs1 , & ! depth of snow layer boundaries (m)
1142 10530496 : zs2 ! adjusted depths, with equal hslyr (m)
1143 :
1144 : real (kind=dbl_kind), dimension (nilyr) :: &
1145 15610614 : dzi ! ice layer thickness after growth/melting
1146 :
1147 : real (kind=dbl_kind), dimension (nslyr) :: &
1148 9043804 : dzs ! snow layer thickness after growth/melting
1149 :
1150 : real (kind=dbl_kind), dimension (nilyr) :: &
1151 16998024 : qm , & ! energy of melting (J m-3) = zqin in BL99 formulation
1152 17097306 : qmlt ! enthalpy of melted ice (J m-3) = zero in BL99 formulation
1153 :
1154 : real (kind=dbl_kind) :: &
1155 1486692 : qbotm , &
1156 1486692 : qbotp , &
1157 1486692 : qbot0
1158 :
1159 : character(len=*),parameter :: subname='(thickness_changes)'
1160 :
1161 : !-----------------------------------------------------------------
1162 : ! Initialize
1163 : !-----------------------------------------------------------------
1164 :
1165 4360794 : dhi = c0
1166 4360794 : dhs = c0
1167 4360794 : hsn_new = c0
1168 :
1169 33066714 : do k = 1, nilyr
1170 33066714 : dzi(k) = hilyr
1171 : enddo
1172 :
1173 9630500 : do k = 1, nslyr
1174 9630500 : dzs(k) = hslyr
1175 : enddo
1176 :
1177 33066714 : do k = 1, nilyr
1178 28705920 : if (ktherm == 2) then
1179 22207724 : qmlt(k) = enthalpy_of_melting(zSin(k))
1180 : else
1181 6498196 : qmlt(k) = c0
1182 : endif
1183 28705920 : qm(k) = zqin(k) - qmlt(k)
1184 28705920 : emlt_atm = c0
1185 33066714 : emlt_ocn = c0
1186 : enddo
1187 :
1188 : !-----------------------------------------------------------------
1189 : ! For l_brine = false (fresh ice), check for temperatures > 0.
1190 : ! Melt ice or snow as needed to bring temperatures back to 0.
1191 : ! For l_brine = true, this should not be necessary.
1192 : !-----------------------------------------------------------------
1193 :
1194 4360794 : if (.not. l_brine) then
1195 :
1196 606546 : do k = 1, nslyr
1197 303273 : Ts = (Lfresh + zqsn(k)/rhos) / cp_ice
1198 606546 : if (Ts > c0) then
1199 0 : dhs = cp_ice*Ts*dzs(k) / Lfresh
1200 0 : dzs(k) = dzs(k) - dhs
1201 0 : zqsn(k) = -rhos*Lfresh
1202 : endif
1203 : enddo
1204 :
1205 606546 : do k = 1, nilyr
1206 303273 : Ti = (Lfresh + zqin(k)/rhoi) / cp_ice
1207 606546 : if (Ti > c0) then
1208 0 : dhi = cp_ice*Ti*dzi(k) / Lfresh
1209 0 : dzi(k) = dzi(k) - dhi
1210 0 : zqin(k) = -rhoi*Lfresh
1211 : endif
1212 : enddo ! k
1213 :
1214 : endif ! .not. l_brine
1215 :
1216 : !-----------------------------------------------------------------
1217 : ! Compute energy available for sublimation/condensation, top melt,
1218 : ! and bottom growth/melt.
1219 : !-----------------------------------------------------------------
1220 :
1221 4360794 : wk1 = -flatn * dt
1222 4360794 : esub = max(wk1, c0) ! energy for sublimation, > 0
1223 4360794 : econ = min(wk1, c0) ! energy for condensation, < 0
1224 :
1225 4360794 : wk1 = (fsurfn - fcondtopn) * dt
1226 4360794 : etop_mlt = max(wk1, c0) ! etop_mlt > 0
1227 :
1228 4360794 : wk1 = (fcondbotn - fbot) * dt
1229 4360794 : ebot_mlt = max(wk1, c0) ! ebot_mlt > 0
1230 4360794 : ebot_gro = min(wk1, c0) ! ebot_gro < 0
1231 :
1232 : !--------------------------------------------------------------
1233 : ! Condensation (evapn > 0)
1234 : ! Note: evapn here has unit of kg/m^2. Divide by dt later.
1235 : ! This is the only case in which energy from the atmosphere
1236 : ! is used for changes in the brine energy (emlt_atm).
1237 : !--------------------------------------------------------------
1238 :
1239 4360794 : evapn = c0 ! initialize
1240 4360794 : evapsn = c0 ! initialize
1241 4360794 : evapin = c0 ! initialize
1242 :
1243 4360794 : if (hsn > puny) then ! add snow with enthalpy zqsn(1)
1244 3827465 : dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0
1245 3827465 : dzs(1) = dzs(1) + dhs
1246 3827465 : evapn = evapn + dhs*rhos
1247 3827465 : evapsn = evapsn + dhs*rhos
1248 : else ! add ice with enthalpy zqin(1)
1249 533329 : dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0
1250 533329 : dzi(1) = dzi(1) + dhi
1251 533329 : evapn = evapn + dhi*rhoi
1252 533329 : evapin = evapin + dhi*rhoi
1253 : ! enthalpy of melt water
1254 533329 : emlt_atm = emlt_atm - qmlt(1) * dhi
1255 : endif
1256 :
1257 : !--------------------------------------------------------------
1258 : ! Grow ice (bottom)
1259 : !--------------------------------------------------------------
1260 :
1261 4360794 : if (ktherm == 2) then
1262 :
1263 3172532 : qbotm = enthalpy_mush(Tbot, sss)
1264 3172532 : qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
1265 3172532 : qbot0 = qbotm - qbotp
1266 :
1267 3172532 : dhi = ebot_gro / qbotp ! dhi > 0
1268 :
1269 3172532 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
1270 3172532 : hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
1271 3172532 : emlt_ocn = emlt_ocn - qbot0 * dhi
1272 :
1273 : else
1274 :
1275 1188262 : Tmlts = -zSin(nilyr) * depressT
1276 :
1277 : ! enthalpy of new ice growing at bottom surface
1278 1188262 : if (heat_capacity) then
1279 884989 : if (l_brine) then
1280 884989 : qbotmax = -p5*rhoi*Lfresh ! max enthalpy of ice growing at bottom
1281 : qbot = -rhoi * (cp_ice * (Tmlts-Tbot) &
1282 : + Lfresh * (c1-Tmlts/Tbot) &
1283 884989 : - cp_ocn * Tmlts)
1284 884989 : qbot = min (qbot, qbotmax) ! in case Tbot is close to Tmlt
1285 : else
1286 0 : qbot = -rhoi * (-cp_ice * Tbot + Lfresh)
1287 : endif
1288 : else ! zero layer
1289 303273 : qbot = -rhoi * Lfresh
1290 : endif
1291 :
1292 1188262 : dhi = ebot_gro / qbot ! dhi > 0
1293 :
1294 1188262 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbot
1295 1188262 : hstot = c0
1296 : endif ! ktherm
1297 :
1298 4360794 : dzi(nilyr) = dzi(nilyr) + dhi
1299 4360794 : if (dzi(nilyr) > puny) then
1300 4360794 : zqin(nilyr) = hqtot / dzi(nilyr)
1301 4360794 : if (ktherm == 2) then
1302 3172532 : zSin(nilyr) = hstot / dzi(nilyr)
1303 3172532 : qmlt(nilyr) = enthalpy_of_melting(zSin(nilyr))
1304 : else
1305 1188262 : qmlt(nilyr) = c0
1306 : endif
1307 4360794 : qm(nilyr) = zqin(nilyr) - qmlt(nilyr)
1308 : endif
1309 :
1310 : ! update ice age due to freezing (new ice age = dt)
1311 : ! if (tr_iage) &
1312 : ! iage = (iage*hin + dt*dhi) / (hin + dhi)
1313 :
1314 : ! history diagnostics
1315 4360794 : congel = congel + dhi
1316 4360794 : if (dhi > puny .and. frz_onset < puny) &
1317 94 : frz_onset = yday
1318 :
1319 9630500 : do k = 1, nslyr
1320 :
1321 : !--------------------------------------------------------------
1322 : ! Remove internal snow melt
1323 : !--------------------------------------------------------------
1324 :
1325 5269706 : if (ktherm == 2 .and. zqsn(k) > -rhos * Lfresh) then
1326 :
1327 23270 : dhs = max(-dzs(k), &
1328 56686 : -((zqsn(k) + rhos*Lfresh) / (rhos*Lfresh)) * dzs(k))
1329 56686 : dzs(k) = dzs(k) + dhs
1330 56686 : zqsn(k) = -rhos * Lfresh
1331 56686 : melts = melts - dhs
1332 : ! delta E = zqsn(k) + rhos * Lfresh
1333 :
1334 : endif
1335 :
1336 : !--------------------------------------------------------------
1337 : ! Sublimation of snow (evapn < 0)
1338 : !--------------------------------------------------------------
1339 :
1340 5269706 : qsub = zqsn(k) - rhos*Lvap ! qsub < 0
1341 5269706 : dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0
1342 5269706 : dzs(k) = dzs(k) + dhs
1343 5269706 : esub = esub - dhs*qsub
1344 5269706 : esub = max(esub, c0) ! in case of roundoff error
1345 5269706 : evapn = evapn + dhs*rhos
1346 5269706 : evapsn = evapsn + dhs*rhos
1347 :
1348 : !--------------------------------------------------------------
1349 : ! Melt snow (top)
1350 : !--------------------------------------------------------------
1351 :
1352 5269706 : dhs = max(-dzs(k), etop_mlt/zqsn(k))
1353 5269706 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1354 5269706 : etop_mlt = etop_mlt - dhs*zqsn(k)
1355 5269706 : etop_mlt = max(etop_mlt, c0) ! in case of roundoff error
1356 :
1357 : ! history diagnostics
1358 5269706 : if (dhs < -puny .and. mlt_onset < puny) &
1359 22 : mlt_onset = yday
1360 9630500 : melts = melts - dhs
1361 :
1362 : enddo ! nslyr
1363 :
1364 33066714 : do k = 1, nilyr
1365 :
1366 : !--------------------------------------------------------------
1367 : ! Sublimation of ice (evapn < 0)
1368 : !--------------------------------------------------------------
1369 :
1370 28705920 : qsub = qm(k) - rhoi*Lvap ! qsub < 0
1371 28705920 : dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0
1372 28705920 : dzi(k) = dzi(k) + dhi
1373 28705920 : esub = esub - dhi*qsub
1374 28705920 : esub = max(esub, c0)
1375 28705920 : evapn = evapn + dhi*rhoi
1376 28705920 : evapin = evapin + dhi*rhoi
1377 28705920 : emlt_ocn = emlt_ocn - qmlt(k) * dhi
1378 :
1379 : !--------------------------------------------------------------
1380 : ! Melt ice (top)
1381 : !--------------------------------------------------------------
1382 :
1383 28705920 : if (qm(k) < c0) then
1384 28705771 : dhi = max(-dzi(k), etop_mlt/qm(k))
1385 : else
1386 149 : qm(k) = c0
1387 149 : dhi = -dzi(k)
1388 : endif
1389 28705920 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1390 :
1391 28705920 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1392 28705920 : etop_mlt = max(etop_mlt - dhi*qm(k), c0)
1393 :
1394 : ! history diagnostics
1395 28705920 : if (dhi < -puny .and. mlt_onset < puny) &
1396 92 : mlt_onset = yday
1397 33066714 : meltt = meltt - dhi
1398 :
1399 : enddo ! nilyr
1400 :
1401 33066714 : do k = nilyr, 1, -1
1402 :
1403 : !--------------------------------------------------------------
1404 : ! Melt ice (bottom)
1405 : !--------------------------------------------------------------
1406 :
1407 28705920 : if (qm(k) < c0) then
1408 28705771 : dhi = max(-dzi(k), ebot_mlt/qm(k))
1409 : else
1410 149 : qm(k) = c0
1411 149 : dhi = -dzi(k)
1412 : endif
1413 28705920 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1414 :
1415 28705920 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1416 28705920 : ebot_mlt = max(ebot_mlt - dhi*qm(k), c0)
1417 :
1418 : ! history diagnostics
1419 33066714 : meltb = meltb -dhi
1420 :
1421 : enddo ! nilyr
1422 :
1423 9630500 : do k = nslyr, 1, -1
1424 :
1425 : !--------------------------------------------------------------
1426 : ! Melt snow (only if all the ice has melted)
1427 : !--------------------------------------------------------------
1428 :
1429 5269706 : dhs = max(-dzs(k), ebot_mlt/zqsn(k))
1430 5269706 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1431 5269706 : ebot_mlt = ebot_mlt - dhs*zqsn(k)
1432 9630500 : ebot_mlt = max(ebot_mlt, c0)
1433 :
1434 : enddo ! nslyr
1435 :
1436 : !-----------------------------------------------------------------
1437 : ! Compute heat flux used by the ice (<=0).
1438 : ! fhocn is the available ocean heat that is left after use by ice
1439 : !-----------------------------------------------------------------
1440 :
1441 : fhocnn = fbot &
1442 4360794 : + (esub + etop_mlt + ebot_mlt)/dt
1443 :
1444 : !---!-----------------------------------------------------------------
1445 : !---! Add new snowfall at top surface.
1446 : !---!-----------------------------------------------------------------
1447 :
1448 : !----------------------------------------------------------------
1449 : ! NOTE: If heat flux diagnostics are to work, new snow should
1450 : ! have T = 0 (i.e. q = -rhos*Lfresh) and should not be
1451 : ! converted to rain.
1452 : !----------------------------------------------------------------
1453 :
1454 4360794 : if (fsnow > c0) then
1455 :
1456 3475124 : hsn_new = fsnow/rhos * dt
1457 3475124 : zqsnew = -rhos*Lfresh
1458 3475124 : hstot = dzs(1) + hsn_new
1459 :
1460 3475124 : if (hstot > c0) then
1461 1146473 : zqsn(1) = (dzs(1) * zqsn(1) &
1462 3475124 : + hsn_new * zqsnew) / hstot
1463 : ! avoid roundoff errors
1464 3475124 : zqsn(1) = min(zqsn(1), -rhos*Lfresh)
1465 :
1466 3475124 : dzs(1) = hstot
1467 : endif
1468 : endif
1469 :
1470 : !-----------------------------------------------------------------
1471 : ! Find the new ice and snow thicknesses.
1472 : !-----------------------------------------------------------------
1473 :
1474 4360794 : hin = c0
1475 4360794 : hsn = c0
1476 :
1477 33066714 : do k = 1, nilyr
1478 33066714 : hin = hin + dzi(k)
1479 : enddo ! k
1480 :
1481 9630500 : do k = 1, nslyr
1482 5269706 : hsn = hsn + dzs(k)
1483 9630500 : dsnow = dsnow + dzs(k) - hslyr
1484 : enddo ! k
1485 :
1486 : !-------------------------------------------------------------------
1487 : ! Convert snow to ice if snow lies below freeboard.
1488 : !-------------------------------------------------------------------
1489 :
1490 4360794 : if (ktherm /= 2) &
1491 : call freeboard (nslyr, &
1492 : snoice, &
1493 : hin, hsn, &
1494 0 : zqin, zqsn, &
1495 0 : dzi, dzs, &
1496 1188262 : dsnow)
1497 4360794 : if (icepack_warnings_aborted(subname)) return
1498 :
1499 : !---!-------------------------------------------------------------------
1500 : !---! Repartition the ice and snow into equal-thickness layers,
1501 : !---! conserving energy.
1502 : !---!-------------------------------------------------------------------
1503 :
1504 : !-----------------------------------------------------------------
1505 : ! Compute desired layer thicknesses.
1506 : !-----------------------------------------------------------------
1507 :
1508 4360794 : if (hin > c0) then
1509 4360787 : hilyr = hin / real(nilyr,kind=dbl_kind)
1510 : else
1511 7 : hin = c0
1512 7 : hilyr = c0
1513 : endif
1514 4360794 : if (hsn > c0) then
1515 3806634 : hslyr = hsn / real(nslyr,kind=dbl_kind)
1516 : else
1517 554160 : hsn = c0
1518 554160 : hslyr = c0
1519 : endif
1520 :
1521 : !-----------------------------------------------------------------
1522 : ! Compute depths zi1 of old layers (unequal thickness).
1523 : ! Compute depths zi2 of new layers (equal thickness).
1524 : !-----------------------------------------------------------------
1525 :
1526 4360794 : zi1(1) = c0
1527 4360794 : zi1(1+nilyr) = hin
1528 :
1529 4360794 : zi2(1) = c0
1530 4360794 : zi2(1+nilyr) = hin
1531 :
1532 4360794 : if (heat_capacity) then
1533 :
1534 28402647 : do k = 1, nilyr-1
1535 24345126 : zi1(k+1) = zi1(k) + dzi(k)
1536 28402647 : zi2(k+1) = zi2(k) + hilyr
1537 : enddo
1538 :
1539 : !-----------------------------------------------------------------
1540 : ! Conserving energy, compute the enthalpy of the new equal layers.
1541 : !-----------------------------------------------------------------
1542 :
1543 : call adjust_enthalpy (nilyr, &
1544 0 : zi1, zi2, &
1545 : hilyr, hin, &
1546 4057521 : zqin)
1547 4057521 : if (icepack_warnings_aborted(subname)) return
1548 :
1549 4057521 : if (ktherm == 2) &
1550 : call adjust_enthalpy (nilyr, &
1551 0 : zi1, zi2, &
1552 : hilyr, hin, &
1553 3172532 : zSin)
1554 4057521 : if (icepack_warnings_aborted(subname)) return
1555 :
1556 : else ! zero layer (nilyr=1)
1557 :
1558 303273 : zqin(1) = -rhoi * Lfresh
1559 303273 : zqsn(1) = -rhos * Lfresh
1560 :
1561 : endif
1562 :
1563 4360794 : if (nslyr > 1) then
1564 :
1565 : !-----------------------------------------------------------------
1566 : ! Compute depths zs1 of old layers (unequal thickness).
1567 : ! Compute depths zs2 of new layers (equal thickness).
1568 : !-----------------------------------------------------------------
1569 :
1570 227228 : zs1(1) = c0
1571 227228 : zs1(1+nslyr) = hsn
1572 :
1573 227228 : zs2(1) = c0
1574 227228 : zs2(1+nslyr) = hsn
1575 :
1576 1136140 : do k = 1, nslyr-1
1577 908912 : zs1(k+1) = zs1(k) + dzs(k)
1578 1136140 : zs2(k+1) = zs2(k) + hslyr
1579 : enddo
1580 :
1581 : !-----------------------------------------------------------------
1582 : ! Conserving energy, compute the enthalpy of the new equal layers.
1583 : !-----------------------------------------------------------------
1584 :
1585 : call adjust_enthalpy (nslyr, &
1586 0 : zs1, zs2, &
1587 : hslyr, hsn, &
1588 227228 : zqsn)
1589 227228 : if (icepack_warnings_aborted(subname)) return
1590 :
1591 : endif ! nslyr > 1
1592 :
1593 : !-----------------------------------------------------------------
1594 : ! Remove very thin snow layers (ktherm = 2)
1595 : !-----------------------------------------------------------------
1596 :
1597 4360794 : if (ktherm == 2) then
1598 6345064 : do k = 1, nslyr
1599 6345064 : if (hsn <= puny) then
1600 : fhocnn = fhocnn &
1601 76923 : + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
1602 76923 : zqsn(k) = -rhos*Lfresh
1603 76923 : hslyr = c0
1604 : endif
1605 : enddo
1606 : endif
1607 :
1608 : !-----------------------------------------------------------------
1609 : ! Compute final ice-snow energy, including the energy of
1610 : ! sublimated/condensed ice.
1611 : !-----------------------------------------------------------------
1612 :
1613 4360794 : efinal = -evapn*Lvap
1614 4360794 : evapn = evapn/dt
1615 4360794 : evapsn = evapsn/dt
1616 4360794 : evapin = evapin/dt
1617 :
1618 9630500 : do k = 1, nslyr
1619 9630500 : efinal = efinal + hslyr*zqsn(k)
1620 : enddo
1621 :
1622 33066714 : do k = 1, nilyr
1623 33066714 : efinal = efinal + hilyr*zqin(k)
1624 : enddo ! k
1625 :
1626 4360794 : if (ktherm < 2) then
1627 1188262 : emlt_atm = c0
1628 1188262 : emlt_ocn = c0
1629 : endif
1630 :
1631 : ! melt water is no longer zero enthalpy with ktherm=2
1632 4360794 : fhocnn = fhocnn + emlt_ocn/dt
1633 4360794 : efinal = efinal + emlt_atm ! for conservation check
1634 :
1635 : end subroutine thickness_changes
1636 :
1637 : !=======================================================================
1638 : !
1639 : ! If there is enough snow to lower the ice/snow interface below
1640 : ! sea level, convert enough snow to ice to bring the interface back
1641 : ! to sea level.
1642 : !
1643 : ! authors William H. Lipscomb, LANL
1644 : ! Elizabeth C. Hunke, LANL
1645 :
1646 1188262 : subroutine freeboard (nslyr, &
1647 : snoice, &
1648 : hin, hsn, &
1649 2376524 : zqin, zqsn, &
1650 2376524 : dzi, dzs, &
1651 : dsnow)
1652 :
1653 : integer (kind=int_kind), intent(in) :: &
1654 : nslyr ! number of snow layers
1655 :
1656 : ! real (kind=dbl_kind), intent(in) :: &
1657 : ! dt ! time step
1658 :
1659 : real (kind=dbl_kind), &
1660 : intent(inout) :: &
1661 : snoice , & ! snow-ice formation (m/step-->cm/day)
1662 : dsnow ! change in snow thickness after snow-ice formation (m)
1663 : ! iage ! ice age (s)
1664 :
1665 : real (kind=dbl_kind), &
1666 : intent(inout) :: &
1667 : hin , & ! ice thickness (m)
1668 : hsn ! snow thickness (m)
1669 :
1670 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1671 : zqsn ! snow layer enthalpy (J m-3)
1672 :
1673 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1674 : zqin , & ! ice layer enthalpy (J m-3)
1675 : dzi , & ! ice layer thicknesses (m)
1676 : dzs ! snow layer thicknesses (m)
1677 :
1678 : ! local variables
1679 :
1680 : integer (kind=int_kind) :: &
1681 : k ! vertical index
1682 :
1683 : real (kind=dbl_kind) :: &
1684 422506 : dhin , & ! change in ice thickness (m)
1685 422506 : dhsn , & ! change in snow thickness (m)
1686 422506 : hqs ! sum of h*q for snow (J m-2)
1687 :
1688 : real (kind=dbl_kind) :: &
1689 422506 : wk1 , & ! temporary variable
1690 422506 : dhs ! snow to remove from layer (m)
1691 :
1692 : character(len=*),parameter :: subname='(freeboard)'
1693 :
1694 : !-----------------------------------------------------------------
1695 : ! Determine whether snow lies below freeboard.
1696 : !-----------------------------------------------------------------
1697 :
1698 1188262 : dhin = c0
1699 1188262 : dhsn = c0
1700 1188262 : hqs = c0
1701 :
1702 1188262 : wk1 = hsn - hin*(rhow-rhoi)/rhos
1703 :
1704 1188262 : if (wk1 > puny .and. hsn > puny) then ! snow below freeboard
1705 17248 : dhsn = min(wk1*rhoi/rhow, hsn) ! snow to remove
1706 17248 : dhin = dhsn * rhos/rhoi ! ice to add
1707 : endif
1708 :
1709 : !-----------------------------------------------------------------
1710 : ! Adjust snow layer thickness.
1711 : ! Compute energy to transfer from snow to ice.
1712 : !-----------------------------------------------------------------
1713 :
1714 3285436 : do k = nslyr, 1, -1
1715 3285436 : if (dhin > puny) then
1716 40816 : dhs = min(dhsn, dzs(k)) ! snow to remove from layer
1717 40816 : hsn = hsn - dhs
1718 40816 : dsnow = dsnow -dhs !new snow addition term
1719 40816 : dzs(k) = dzs(k) - dhs
1720 40816 : dhsn = dhsn - dhs
1721 40816 : dhsn = max(dhsn,c0)
1722 40816 : hqs = hqs + dhs * zqsn(k)
1723 : endif ! dhin > puny
1724 : enddo
1725 :
1726 : !-----------------------------------------------------------------
1727 : ! Transfer volume and energy from snow to top ice layer.
1728 : !-----------------------------------------------------------------
1729 :
1730 1188262 : if (dhin > puny) then
1731 : ! update ice age due to freezing (new ice age = dt)
1732 : ! if (tr_iage) &
1733 : ! iage = (iage*hin+dt*dhin)/(hin+dhin)
1734 :
1735 17248 : wk1 = dzi(1) + dhin
1736 17248 : hin = hin + dhin
1737 17248 : zqin(1) = (dzi(1)*zqin(1) + hqs) / wk1
1738 17248 : dzi(1) = wk1
1739 :
1740 : ! history diagnostic
1741 17248 : snoice = snoice + dhin
1742 : endif ! dhin > puny
1743 :
1744 1188262 : end subroutine freeboard
1745 :
1746 : !=======================================================================
1747 : !
1748 : ! Conserving energy, compute the new enthalpy of equal-thickness ice
1749 : ! or snow layers.
1750 : !
1751 : ! authors William H. Lipscomb, LANL
1752 : ! C. M. Bitz, UW
1753 :
1754 7457281 : subroutine adjust_enthalpy (nlyr, &
1755 7457281 : z1, z2, &
1756 : hlyr, hn, &
1757 7457281 : qn)
1758 :
1759 : integer (kind=int_kind), intent(in) :: &
1760 : nlyr ! number of layers (nilyr or nslyr)
1761 :
1762 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1763 : z1 , & ! interface depth for old, unequal layers (m)
1764 : z2 ! interface depth for new, equal layers (m)
1765 :
1766 : real (kind=dbl_kind), intent(in) :: &
1767 : hlyr ! new layer thickness (m)
1768 :
1769 : real (kind=dbl_kind), intent(in) :: &
1770 : hn ! total thickness (m)
1771 :
1772 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1773 : qn ! layer quantity (enthalpy, salinity...)
1774 :
1775 : ! local variables
1776 :
1777 : integer (kind=int_kind) :: &
1778 : k, k1, k2 ! vertical indices
1779 :
1780 : real (kind=dbl_kind) :: &
1781 2524146 : hovlp ! overlap between old and new layers (m)
1782 :
1783 : real (kind=dbl_kind) :: &
1784 2524146 : rhlyr ! 1./hlyr
1785 :
1786 : real (kind=dbl_kind), dimension (nlyr) :: &
1787 29898330 : hq ! h * q for a layer
1788 :
1789 : character(len=*),parameter :: subname='(adjust_enthalpy)'
1790 :
1791 : !-----------------------------------------------------------------
1792 : ! Compute reciprocal layer thickness.
1793 : !-----------------------------------------------------------------
1794 :
1795 7457281 : rhlyr = c0
1796 7457281 : if (hn > puny) rhlyr = c1 / hlyr
1797 :
1798 : !-----------------------------------------------------------------
1799 : ! Compute h*q for new layers (k2) given overlap with old layers (k1)
1800 : !-----------------------------------------------------------------
1801 :
1802 59203792 : do k2 = 1, nlyr
1803 59203792 : hq(k2) = c0
1804 : enddo ! k
1805 7457281 : k1 = 1
1806 7457281 : k2 = 1
1807 102711834 : do while (k1 <= nlyr .and. k2 <= nlyr)
1808 64466692 : hovlp = min (z1(k1+1), z2(k2+1)) &
1809 159721245 : - max (z1(k1), z2(k2))
1810 95254553 : hovlp = max (hovlp, c0)
1811 95254553 : hq(k2) = hq(k2) + hovlp*qn(k1)
1812 102711834 : if (z1(k1+1) > z2(k2+1)) then
1813 43508042 : k2 = k2 + 1
1814 : else
1815 51746511 : k1 = k1 + 1
1816 : endif
1817 : enddo ! while
1818 :
1819 : !-----------------------------------------------------------------
1820 : ! Compute new enthalpies.
1821 : !-----------------------------------------------------------------
1822 :
1823 59203792 : do k = 1, nlyr
1824 59203792 : qn(k) = hq(k) * rhlyr
1825 : enddo ! k
1826 :
1827 7457281 : end subroutine adjust_enthalpy
1828 :
1829 : !=======================================================================
1830 : !
1831 : ! Check for energy conservation by comparing the change in energy
1832 : ! to the net energy input.
1833 : !
1834 : ! authors William H. Lipscomb, LANL
1835 : ! C. M. Bitz, UW
1836 : ! Adrian K. Turner, LANL
1837 :
1838 4360794 : subroutine conservation_check_vthermo(dt, &
1839 : fsurfn, flatn, &
1840 : fhocnn, fswint, &
1841 : fsnow, &
1842 : einit, einter, &
1843 : efinal, &
1844 : fcondtopn,fcondbotn, &
1845 : fadvocn, fbot )
1846 :
1847 : real (kind=dbl_kind), intent(in) :: &
1848 : dt ! time step
1849 :
1850 : real (kind=dbl_kind), intent(in) :: &
1851 : fsurfn , & ! net flux to top surface, excluding fcondtopn
1852 : flatn , & ! surface downward latent heat (W m-2)
1853 : fhocnn , & ! fbot, corrected for any surplus energy
1854 : fswint , & ! SW absorbed in ice interior, below surface (W m-2)
1855 : fsnow , & ! snowfall rate (kg m-2 s-1)
1856 : fcondtopn , &
1857 : fadvocn , &
1858 : fbot
1859 :
1860 : real (kind=dbl_kind), intent(in) :: &
1861 : einit , & ! initial energy of melting (J m-2)
1862 : einter , & ! intermediate energy of melting (J m-2)
1863 : efinal , & ! final energy of melting (J m-2)
1864 : fcondbotn
1865 :
1866 : ! local variables
1867 :
1868 : real (kind=dbl_kind) :: &
1869 1486692 : einp , & ! energy input during timestep (J m-2)
1870 1486692 : ferr ! energy conservation error (W m-2)
1871 :
1872 : character(len=*),parameter :: subname='(conservation_check_vthermo)'
1873 :
1874 : !----------------------------------------------------------------
1875 : ! If energy is not conserved, print diagnostics and exit.
1876 : !----------------------------------------------------------------
1877 :
1878 : !-----------------------------------------------------------------
1879 : ! Note that fsurf - flat = fsw + flw + fsens; i.e., the latent
1880 : ! heat is not included in the energy input, since (efinal - einit)
1881 : ! is the energy change in the system ice + vapor, and the latent
1882 : ! heat lost by the ice is equal to that gained by the vapor.
1883 : !-----------------------------------------------------------------
1884 :
1885 : einp = (fsurfn - flatn + fswint - fhocnn &
1886 4360794 : - fsnow*Lfresh - fadvocn) * dt
1887 4360794 : ferr = abs(efinal-einit-einp) / dt
1888 :
1889 4360794 : if (ferr > ferrmax) then
1890 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1891 0 : call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" )
1892 :
1893 0 : write(warnstr,*) subname, 'Thermo energy conservation error'
1894 0 : call icepack_warnings_add(warnstr)
1895 0 : write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr
1896 0 : call icepack_warnings_add(warnstr)
1897 0 : write(warnstr,*) subname, 'Energy error (J) =', ferr*dt
1898 0 : call icepack_warnings_add(warnstr)
1899 0 : write(warnstr,*) subname, 'Initial energy =', einit
1900 0 : call icepack_warnings_add(warnstr)
1901 0 : write(warnstr,*) subname, 'Final energy =', efinal
1902 0 : call icepack_warnings_add(warnstr)
1903 0 : write(warnstr,*) subname, 'efinal - einit =', efinal-einit
1904 0 : call icepack_warnings_add(warnstr)
1905 0 : write(warnstr,*) subname, 'fsurfn,flatn,fswint,fhocn, fsnow*Lfresh:'
1906 0 : call icepack_warnings_add(warnstr)
1907 0 : write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn, fsnow*Lfresh
1908 0 : call icepack_warnings_add(warnstr)
1909 0 : write(warnstr,*) subname, 'Input energy =', einp
1910 0 : call icepack_warnings_add(warnstr)
1911 0 : write(warnstr,*) subname, 'fbot,fcondbot:'
1912 0 : call icepack_warnings_add(warnstr)
1913 0 : write(warnstr,*) subname, fbot,fcondbotn
1914 0 : call icepack_warnings_add(warnstr)
1915 :
1916 : ! if (ktherm == 2) then
1917 0 : write(warnstr,*) subname, 'Intermediate energy =', einter
1918 0 : call icepack_warnings_add(warnstr)
1919 0 : write(warnstr,*) subname, 'efinal - einter =', &
1920 0 : efinal-einter
1921 0 : call icepack_warnings_add(warnstr)
1922 0 : write(warnstr,*) subname, 'einter - einit =', &
1923 0 : einter-einit
1924 0 : call icepack_warnings_add(warnstr)
1925 0 : write(warnstr,*) subname, 'Conduction Error =', (einter-einit) &
1926 0 : - (fcondtopn*dt - fcondbotn*dt + fswint*dt)
1927 0 : call icepack_warnings_add(warnstr)
1928 0 : write(warnstr,*) subname, 'Melt/Growth Error =', (einter-einit) &
1929 0 : + ferr*dt - (fcondtopn*dt - fcondbotn*dt + fswint*dt)
1930 0 : call icepack_warnings_add(warnstr)
1931 0 : write(warnstr,*) subname, 'Advection Error =', fadvocn*dt
1932 0 : call icepack_warnings_add(warnstr)
1933 : ! endif
1934 :
1935 : ! write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn
1936 : ! call icepack_warnings_add(warnstr)
1937 :
1938 0 : write(warnstr,*) subname, 'dt*(fsurfn, flatn, fswint, fhocn, fsnow*Lfresh, fadvocn):'
1939 0 : call icepack_warnings_add(warnstr)
1940 0 : write(warnstr,*) subname, fsurfn*dt, flatn*dt, &
1941 0 : fswint*dt, fhocnn*dt, &
1942 0 : fsnow*Lfresh*dt, fadvocn*dt
1943 0 : call icepack_warnings_add(warnstr)
1944 0 : return
1945 : endif
1946 :
1947 : end subroutine conservation_check_vthermo
1948 :
1949 : !=======================================================================
1950 : !
1951 : ! Given the vertical thermo state variables (hin, hsn),
1952 : ! compute the new ice state variables (vicen, vsnon).
1953 : ! Zero out state variables if ice has melted entirely.
1954 : !
1955 : ! authors William H. Lipscomb, LANL
1956 : ! C. M. Bitz, UW
1957 : ! Elizabeth C. Hunke, LANL
1958 :
1959 4360794 : subroutine update_state_vthermo(nilyr, nslyr, &
1960 : Tf, Tsf, &
1961 : hin, hsn, &
1962 4360794 : zqin, zSin, &
1963 4360794 : zqsn, &
1964 : aicen, vicen, &
1965 : vsnon)
1966 :
1967 : integer (kind=int_kind), intent(in) :: &
1968 : nilyr , & ! number of ice layers
1969 : nslyr ! number of snow layers
1970 :
1971 : real (kind=dbl_kind), intent(in) :: &
1972 : Tf ! freezing temperature (C)
1973 :
1974 : real (kind=dbl_kind), intent(inout) :: &
1975 : Tsf ! ice/snow surface temperature, Tsfcn
1976 :
1977 : real (kind=dbl_kind), intent(in) :: &
1978 : hin , & ! ice thickness (m)
1979 : hsn ! snow thickness (m)
1980 :
1981 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1982 : zqin , & ! ice layer enthalpy (J m-3)
1983 : zSin , & ! ice salinity (ppt)
1984 : zqsn ! snow layer enthalpy (J m-3)
1985 :
1986 : real (kind=dbl_kind), intent(inout) :: &
1987 : aicen , & ! concentration of ice
1988 : vicen , & ! volume per unit area of ice (m)
1989 : vsnon ! volume per unit area of snow (m)
1990 :
1991 : ! local variables
1992 :
1993 : integer (kind=int_kind) :: &
1994 : k ! ice layer index
1995 :
1996 : character(len=*),parameter :: subname='(update_state_vthermo)'
1997 :
1998 4360794 : if (hin <= c0) then
1999 7 : aicen = c0
2000 7 : vicen = c0
2001 7 : vsnon = c0
2002 7 : Tsf = Tf
2003 56 : do k = 1, nilyr
2004 56 : zqin(k) = c0
2005 : enddo
2006 7 : if (ktherm == 2) then
2007 56 : do k = 1, nilyr
2008 56 : zSin(k) = c0
2009 : enddo
2010 : endif
2011 14 : do k = 1, nslyr
2012 14 : zqsn(k) = c0
2013 : enddo
2014 : else
2015 : ! aicen is already up to date
2016 4360787 : vicen = aicen * hin
2017 4360787 : vsnon = aicen * hsn
2018 : endif
2019 :
2020 4360794 : end subroutine update_state_vthermo
2021 :
2022 : !=======================================================================
2023 : !autodocument_start icepack_step_therm1
2024 : ! Driver for thermodynamic changes not needed for coupling:
2025 : ! transport in thickness space, lateral growth and melting.
2026 : !
2027 : ! authors: William H. Lipscomb, LANL
2028 : ! Elizabeth C. Hunke, LANL
2029 :
2030 1370160 : subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
2031 1370160 : aicen_init , &
2032 1370160 : vicen_init , vsnon_init , &
2033 1370160 : aice , aicen , &
2034 1370160 : vice , vicen , &
2035 1370160 : vsno , vsnon , &
2036 : uvel , vvel , &
2037 2740320 : Tsfc , zqsn , &
2038 2740320 : zqin , zSin , &
2039 1370160 : alvl , vlvl , &
2040 1370160 : apnd , hpnd , &
2041 1370160 : ipnd , &
2042 1370160 : iage , FY , &
2043 1370160 : aerosno , aeroice , &
2044 1370160 : isosno , isoice , &
2045 : uatm , vatm , &
2046 : wind , zlvl , &
2047 : Qa , rhoa , &
2048 1370160 : Qa_iso , &
2049 : Tair , Tref , &
2050 : Qref , Uref , &
2051 1370160 : Qref_iso , &
2052 : Cdn_atm_ratio, &
2053 : Cdn_ocn , Cdn_ocn_skin, &
2054 : Cdn_ocn_floe, Cdn_ocn_keel, &
2055 : Cdn_atm , Cdn_atm_skin, &
2056 : Cdn_atm_floe, Cdn_atm_pond, &
2057 : Cdn_atm_rdg , hfreebd , &
2058 : hdraft , hridge , &
2059 : distrdg , hkeel , &
2060 : dkeel , lfloe , &
2061 : dfloe , &
2062 : strax , stray , &
2063 : strairxT , strairyT , &
2064 : potT , sst , &
2065 : sss , Tf , &
2066 : strocnxT , strocnyT , &
2067 : fbot , &
2068 : Tbot , Tsnice , &
2069 : frzmlt , rside , &
2070 : fside , &
2071 : fsnow , frain , &
2072 : fpond , &
2073 1370160 : fsurf , fsurfn , &
2074 1370160 : fcondtop , fcondtopn , &
2075 1370160 : fcondbot , fcondbotn , &
2076 1370160 : fswsfcn , fswintn , &
2077 1370160 : fswthrun , &
2078 1370160 : fswthrun_vdr, &
2079 1370160 : fswthrun_vdf, &
2080 1370160 : fswthrun_idr, &
2081 1370160 : fswthrun_idf, &
2082 : fswabs , &
2083 : flwout , &
2084 1370160 : Sswabsn , Iswabsn , &
2085 : flw , &
2086 1370160 : fsens , fsensn , &
2087 1370160 : flat , flatn , &
2088 : evap , &
2089 : evaps , evapi , &
2090 : fresh , fsalt , &
2091 : fhocn , &
2092 : fswthru , &
2093 : fswthru_vdr , &
2094 : fswthru_vdf , &
2095 : fswthru_idr , &
2096 : fswthru_idf , &
2097 1370160 : flatn_f , fsensn_f , &
2098 1370160 : fsurfn_f , fcondtopn_f , &
2099 1370160 : faero_atm , faero_ocn , &
2100 1370160 : fiso_atm , fiso_ocn , &
2101 1370160 : fiso_evap , &
2102 : HDO_ocn , H2_16O_ocn , &
2103 : H2_18O_ocn , &
2104 1370160 : dhsn , ffracn , &
2105 1370160 : meltt , melttn , &
2106 1370160 : meltb , meltbn , &
2107 1370160 : melts , meltsn , &
2108 1370160 : congel , congeln , &
2109 1370160 : snoice , snoicen , &
2110 1370160 : dsnown , &
2111 : lmask_n , lmask_s , &
2112 : mlt_onset , frz_onset , &
2113 : yday , prescribed_ice)
2114 :
2115 : integer (kind=int_kind), intent(in) :: &
2116 : ncat , & ! number of thickness categories
2117 : nilyr , & ! number of ice layers
2118 : nslyr ! number of snow layers
2119 :
2120 : real (kind=dbl_kind), intent(in) :: &
2121 : dt , & ! time step
2122 : uvel , & ! x-component of velocity (m/s)
2123 : vvel , & ! y-component of velocity (m/s)
2124 : strax , & ! wind stress components (N/m^2)
2125 : stray , & !
2126 : yday ! day of year
2127 :
2128 : logical (kind=log_kind), intent(in) :: &
2129 : lmask_n , & ! northern hemisphere mask
2130 : lmask_s ! southern hemisphere mask
2131 :
2132 : logical (kind=log_kind), intent(in), optional :: &
2133 : prescribed_ice ! if .true., use prescribed ice instead of computed
2134 :
2135 : real (kind=dbl_kind), intent(inout) :: &
2136 : aice , & ! sea ice concentration
2137 : vice , & ! volume per unit area of ice (m)
2138 : vsno , & ! volume per unit area of snow (m)
2139 : zlvl , & ! atm level height (m)
2140 : uatm , & ! wind velocity components (m/s)
2141 : vatm , &
2142 : wind , & ! wind speed (m/s)
2143 : potT , & ! air potential temperature (K)
2144 : Tair , & ! air temperature (K)
2145 : Qa , & ! specific humidity (kg/kg)
2146 : rhoa , & ! air density (kg/m^3)
2147 : frain , & ! rainfall rate (kg/m^2 s)
2148 : fsnow , & ! snowfall rate (kg/m^2 s)
2149 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
2150 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
2151 : fsalt , & ! salt flux to ocean (kg/m^2/s)
2152 : fhocn , & ! net heat flux to ocean (W/m^2)
2153 : fswthru , & ! shortwave penetrating to ocean (W/m^2)
2154 : fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2)
2155 : fcondtop , & ! top surface conductive flux (W/m^2)
2156 : fcondbot , & ! bottom surface conductive flux (W/m^2)
2157 : fsens , & ! sensible heat flux (W/m^2)
2158 : flat , & ! latent heat flux (W/m^2)
2159 : fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2)
2160 : flw , & ! incoming longwave radiation (W/m^2)
2161 : flwout , & ! outgoing longwave radiation (W/m^2)
2162 : evap , & ! evaporative water flux (kg/m^2/s)
2163 : evaps , & ! evaporative water flux over snow (kg/m^2/s)
2164 : evapi , & ! evaporative water flux over ice (kg/m^2/s)
2165 : congel , & ! basal ice growth (m/step-->cm/day)
2166 : snoice , & ! snow-ice formation (m/step-->cm/day)
2167 : Tref , & ! 2m atm reference temperature (K)
2168 : Qref , & ! 2m atm reference spec humidity (kg/kg)
2169 : Uref , & ! 10m atm reference wind speed (m/s)
2170 : Cdn_atm , & ! atm drag coefficient
2171 : Cdn_ocn , & ! ocn drag coefficient
2172 : hfreebd , & ! freeboard (m)
2173 : hdraft , & ! draft of ice + snow column (Stoessel1993)
2174 : hridge , & ! ridge height
2175 : distrdg , & ! distance between ridges
2176 : hkeel , & ! keel depth
2177 : dkeel , & ! distance between keels
2178 : lfloe , & ! floe length
2179 : dfloe , & ! distance between floes
2180 : Cdn_atm_skin, & ! neutral skin drag coefficient
2181 : Cdn_atm_floe, & ! neutral floe edge drag coefficient
2182 : Cdn_atm_pond, & ! neutral pond edge drag coefficient
2183 : Cdn_atm_rdg , & ! neutral ridge drag coefficient
2184 : Cdn_ocn_skin, & ! skin drag coefficient
2185 : Cdn_ocn_floe, & ! floe edge drag coefficient
2186 : Cdn_ocn_keel, & ! keel drag coefficient
2187 : Cdn_atm_ratio,& ! ratio drag atm / neutral drag atm
2188 : strairxT , & ! stress on ice by air, x-direction
2189 : strairyT , & ! stress on ice by air, y-direction
2190 : strocnxT , & ! ice-ocean stress, x-direction
2191 : strocnyT , & ! ice-ocean stress, y-direction
2192 : fbot , & ! ice-ocean heat flux at bottom surface (W/m^2)
2193 : frzmlt , & ! freezing/melting potential (W/m^2)
2194 : rside , & ! fraction of ice that melts laterally
2195 : fside , & ! lateral heat flux (W/m^2)
2196 : sst , & ! sea surface temperature (C)
2197 : Tf , & ! freezing temperature (C)
2198 : Tbot , & ! ice bottom surface temperature (deg C)
2199 : Tsnice , & ! snow ice interface temperature (deg C)
2200 : sss , & ! sea surface salinity (ppt)
2201 : meltt , & ! top ice melt (m/step-->cm/day)
2202 : melts , & ! snow melt (m/step-->cm/day)
2203 : meltb , & ! basal ice melt (m/step-->cm/day)
2204 : mlt_onset , & ! day of year that sfc melting begins
2205 : frz_onset ! day of year that freezing begins (congel or frazil)
2206 :
2207 : real (kind=dbl_kind), intent(inout), optional :: &
2208 : fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
2209 : fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2)
2210 : fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2)
2211 : fswthru_idf ! nir dif shortwave penetrating to ocean (W/m^2)
2212 :
2213 : real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
2214 : Qa_iso , & ! isotope specific humidity (kg/kg)
2215 : Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg)
2216 : fiso_atm , & ! isotope deposition rate (kg/m^2 s)
2217 : fiso_ocn , & ! isotope flux to ocean (kg/m^2/s)
2218 : fiso_evap ! isotope evaporation (kg/m^2/s)
2219 :
2220 : real (kind=dbl_kind), optional, intent(in) :: &
2221 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
2222 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
2223 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
2224 :
2225 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
2226 : aicen_init , & ! fractional area of ice
2227 : vicen_init , & ! volume per unit area of ice (m)
2228 : vsnon_init , & ! volume per unit area of snow (m)
2229 : aicen , & ! concentration of ice
2230 : vicen , & ! volume per unit area of ice (m)
2231 : vsnon , & ! volume per unit area of snow (m)
2232 : Tsfc , & ! ice/snow surface temperature, Tsfcn
2233 : alvl , & ! level ice area fraction
2234 : vlvl , & ! level ice volume fraction
2235 : apnd , & ! melt pond area fraction
2236 : hpnd , & ! melt pond depth (m)
2237 : ipnd , & ! melt pond refrozen lid thickness (m)
2238 : iage , & ! volume-weighted ice age
2239 : FY , & ! area-weighted first-year ice area
2240 : fsurfn , & ! net flux to top surface, excluding fcondtop
2241 : fcondtopn , & ! downward cond flux at top surface (W m-2)
2242 : fcondbotn , & ! downward cond flux at bottom surface (W m-2)
2243 : flatn , & ! latent heat flux (W m-2)
2244 : fsensn , & ! sensible heat flux (W m-2)
2245 : fsurfn_f , & ! net flux to top surface, excluding fcondtop
2246 : fcondtopn_f , & ! downward cond flux at top surface (W m-2)
2247 : flatn_f , & ! latent heat flux (W m-2)
2248 : fsensn_f , & ! sensible heat flux (W m-2)
2249 : fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
2250 : fswthrun , & ! SW through ice to ocean (W/m^2)
2251 : fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
2252 : faero_atm , & ! aerosol deposition rate (kg/m^2 s)
2253 : faero_ocn , & ! aerosol flux to ocean (kg/m^2/s)
2254 : dhsn , & ! depth difference for snow on sea ice and pond ice
2255 : ffracn , & ! fraction of fsurfn used to melt ipond
2256 : meltsn , & ! snow melt (m)
2257 : melttn , & ! top ice melt (m)
2258 : meltbn , & ! bottom ice melt (m)
2259 : congeln , & ! congelation ice growth (m)
2260 : snoicen , & ! snow-ice growth (m)
2261 : dsnown ! change in snow thickness (m/step-->cm/day)
2262 :
2263 : real (kind=dbl_kind), optional, dimension(:), intent(inout) :: &
2264 : fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
2265 : fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
2266 : fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
2267 : fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
2268 :
2269 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
2270 : zqsn , & ! snow layer enthalpy (J m-3)
2271 : zqin , & ! ice layer enthalpy (J m-3)
2272 : zSin , & ! internal ice layer salinities
2273 : Sswabsn , & ! SW radiation absorbed in snow layers (W m-2)
2274 : Iswabsn ! SW radiation absorbed in ice layers (W m-2)
2275 :
2276 : real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: &
2277 : aerosno , & ! snow aerosol tracer (kg/m^2)
2278 : aeroice ! ice aerosol tracer (kg/m^2)
2279 :
2280 : real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: &
2281 : isosno , & ! snow isotope tracer (kg/m^2)
2282 : isoice ! ice isotope tracer (kg/m^2)
2283 : !autodocument_end
2284 :
2285 : ! local variables
2286 :
2287 : integer (kind=int_kind) :: &
2288 : n ! category index
2289 :
2290 : real (kind=dbl_kind) :: &
2291 474288 : worka, workb ! temporary variables
2292 :
2293 : ! 2D coupler variables (computed for each category, then aggregated)
2294 : real (kind=dbl_kind) :: &
2295 474288 : fswabsn , & ! shortwave absorbed by ice (W/m^2)
2296 474288 : flwoutn , & ! upward LW at surface (W/m^2)
2297 474288 : evapn , & ! flux of vapor, atmos to ice (kg m-2 s-1)
2298 474288 : evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1)
2299 474288 : evapin , & ! flux of vapor, atmos to ice over ice (kg m-2 s-1)
2300 474288 : freshn , & ! flux of water, ice to ocean (kg/m^2/s)
2301 474288 : fsaltn , & ! flux of salt, ice to ocean (kg/m^2/s)
2302 474288 : fhocnn , & ! fbot corrected for leftover energy (W/m^2)
2303 474288 : strairxn , & ! air/ice zonal stress, (N/m^2)
2304 474288 : strairyn , & ! air/ice meridional stress, (N/m^2)
2305 474288 : Cdn_atm_ratio_n, & ! drag coefficient ratio
2306 474288 : Trefn , & ! air tmp reference level (K)
2307 474288 : Urefn , & ! air speed reference level (m/s)
2308 474288 : Qrefn , & ! air sp hum reference level (kg/kg)
2309 474288 : shcoef , & ! transfer coefficient for sensible heat
2310 474288 : lhcoef , & ! transfer coefficient for latent heat
2311 474288 : rfrac ! water fraction retained for melt ponds
2312 :
2313 : real (kind=dbl_kind), dimension(n_iso) :: &
2314 2266032 : Qrefn_iso , & ! isotope air sp hum reference level (kg/kg)
2315 2266032 : fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s)
2316 3214608 : fiso_evapn ! isotope evaporation (kg/m^2/s)
2317 :
2318 : real (kind=dbl_kind), allocatable, dimension(:,:) :: &
2319 1370160 : l_isosno , & ! local snow isotope tracer (kg/m^2)
2320 1370160 : l_isoice ! local ice isotope tracer (kg/m^2)
2321 :
2322 : real (kind=dbl_kind), allocatable, dimension(:) :: &
2323 1370160 : l_Qa_iso , & ! local isotope specific humidity (kg/kg)
2324 1370160 : l_Qref_iso , & ! local isotope 2m atm reference spec humidity (kg/kg)
2325 1370160 : l_fiso_atm , & ! local isotope deposition rate (kg/m^2 s)
2326 1370160 : l_fiso_ocn , & ! local isotope flux to ocean (kg/m^2/s)
2327 1370160 : l_fiso_evap ! local isotope evaporation (kg/m^2/s)
2328 :
2329 : real (kind=dbl_kind) :: &
2330 474288 : l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg)
2331 474288 : l_H2_16O_ocn, & ! local ocean concentration of H2_16O (kg/kg)
2332 474288 : l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg)
2333 :
2334 : real (kind=dbl_kind) :: &
2335 474288 : l_fswthru_vdr , & ! vis dir SW through ice to ocean (W/m^2)
2336 474288 : l_fswthru_vdf , & ! vis dif SW through ice to ocean (W/m^2)
2337 474288 : l_fswthru_idr , & ! nir dir SW through ice to ocean (W/m^2)
2338 474288 : l_fswthru_idf ! nir dif SW through ice to ocean (W/m^2)
2339 :
2340 : real (kind=dbl_kind), dimension(:), allocatable :: &
2341 1370160 : l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
2342 1370160 : l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
2343 1370160 : l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
2344 1370160 : l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
2345 :
2346 : real (kind=dbl_kind) :: &
2347 474288 : pond ! water retained in ponds (m)
2348 :
2349 : character(len=*),parameter :: subname='(icepack_step_therm1)'
2350 :
2351 : !-----------------------------------------------------------------
2352 : ! allocate local optional arguments
2353 : !-----------------------------------------------------------------
2354 :
2355 1370160 : if (present(isosno) ) then
2356 1370160 : allocate(l_isosno(size(isosno,dim=1),size(isosno,dim=2)))
2357 8757168 : l_isosno = isosno
2358 : else
2359 0 : allocate(l_isosno(1,1))
2360 0 : l_isosno = c0
2361 : endif
2362 :
2363 1370160 : if (present(isoice) ) then
2364 1370160 : allocate(l_isoice(size(isoice,dim=1),size(isoice,dim=2)))
2365 8757168 : l_isoice = isoice
2366 : else
2367 0 : allocate(l_isoice(1,1))
2368 0 : l_isoice = c0
2369 : endif
2370 :
2371 1370160 : if (present(Qa_iso) ) then
2372 1370160 : allocate(l_Qa_iso(size(Qa_iso)))
2373 5480640 : l_Qa_iso = Qa_iso
2374 : else
2375 0 : allocate(l_Qa_iso(1))
2376 0 : l_Qa_iso = c0
2377 : endif
2378 :
2379 1370160 : if (present(Qref_iso) ) then
2380 1370160 : allocate(l_Qref_iso(size(Qref_iso)))
2381 5480640 : l_Qref_iso = Qref_iso
2382 : else
2383 0 : allocate(l_Qref_iso(1))
2384 0 : l_Qref_iso = c0
2385 : endif
2386 :
2387 1370160 : if (present(fiso_atm) ) then
2388 1370160 : allocate(l_fiso_atm(size(fiso_atm)))
2389 5480640 : l_fiso_atm = fiso_atm
2390 : else
2391 0 : allocate(l_fiso_atm(1))
2392 0 : l_fiso_atm = c0
2393 : endif
2394 :
2395 1370160 : if (present(fiso_ocn) ) then
2396 1370160 : allocate(l_fiso_ocn(size(fiso_ocn)))
2397 5480640 : l_fiso_ocn = fiso_ocn
2398 : else
2399 0 : allocate(l_fiso_ocn(1))
2400 0 : l_fiso_ocn = c0
2401 : endif
2402 :
2403 1370160 : if (present(fiso_evap) ) then
2404 1370160 : allocate(l_fiso_evap(size(fiso_evap)))
2405 5480640 : l_fiso_evap = fiso_evap
2406 : else
2407 0 : allocate(l_fiso_evap(1))
2408 0 : l_fiso_evap = c0
2409 : endif
2410 :
2411 1370160 : l_HDO_ocn = c0
2412 1370160 : if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn
2413 :
2414 1370160 : l_H2_16O_ocn = c0
2415 1370160 : if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
2416 :
2417 1370160 : l_H2_18O_ocn = c0
2418 1370160 : if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
2419 :
2420 1370160 : l_fswthru_vdr = c0
2421 1370160 : if (present(fswthru_vdr) ) l_fswthru_vdr = fswthru_vdr
2422 :
2423 1370160 : l_fswthru_vdf = c0
2424 1370160 : if (present(fswthru_vdf) ) l_fswthru_vdf = fswthru_vdf
2425 :
2426 1370160 : l_fswthru_idr = c0
2427 1370160 : if (present(fswthru_idr) ) l_fswthru_idr = fswthru_idr
2428 :
2429 1370160 : l_fswthru_idf = c0
2430 1370160 : if (present(fswthru_idf) ) l_fswthru_idf = fswthru_idf
2431 :
2432 1370160 : allocate(l_fswthrun_vdr(ncat))
2433 1370160 : allocate(l_fswthrun_vdf(ncat))
2434 1370160 : allocate(l_fswthrun_idr(ncat))
2435 1370160 : allocate(l_fswthrun_idf(ncat))
2436 :
2437 7834848 : l_fswthrun_vdr = c0
2438 7834848 : if (present(fswthrun_vdr) ) l_fswthrun_vdr = fswthrun_vdr
2439 :
2440 7834848 : l_fswthrun_vdf = c0
2441 7834848 : if (present(fswthrun_vdf) ) l_fswthrun_vdf = fswthrun_vdf
2442 :
2443 7834848 : l_fswthrun_idr = c0
2444 7834848 : if (present(fswthrun_idr) ) l_fswthrun_idr = fswthrun_idr
2445 :
2446 7834848 : l_fswthrun_idf = c0
2447 7834848 : if (present(fswthrun_idf) ) l_fswthrun_idf = fswthrun_idf
2448 :
2449 : !-----------------------------------------------------------------
2450 : ! Adjust frzmlt to account for ice-ocean heat fluxes since last
2451 : ! call to coupler.
2452 : ! Compute lateral and bottom heat fluxes.
2453 : !-----------------------------------------------------------------
2454 :
2455 : call frzmlt_bottom_lateral (dt, ncat, &
2456 : nilyr, nslyr, &
2457 : aice, frzmlt, &
2458 0 : vicen, vsnon, &
2459 0 : zqin, zqsn, &
2460 : sst, Tf, &
2461 : ustar_min, &
2462 : fbot_xfer_type, &
2463 : strocnxT, strocnyT, &
2464 : Tbot, fbot, &
2465 : rside, Cdn_ocn, &
2466 1370160 : fside)
2467 :
2468 1370160 : if (icepack_warnings_aborted(subname)) return
2469 :
2470 : !-----------------------------------------------------------------
2471 : ! Update the neutral drag coefficients to account for form drag
2472 : ! Oceanic and atmospheric drag coefficients
2473 : !-----------------------------------------------------------------
2474 :
2475 1370160 : if (formdrag) then
2476 0 : call neutral_drag_coeffs (apnd , &
2477 0 : hpnd , ipnd , &
2478 0 : alvl , vlvl , &
2479 : aice , vice, &
2480 0 : vsno , aicen , &
2481 0 : vicen , &
2482 : Cdn_ocn , Cdn_ocn_skin, &
2483 : Cdn_ocn_floe, Cdn_ocn_keel, &
2484 : Cdn_atm , Cdn_atm_skin, &
2485 : Cdn_atm_floe, Cdn_atm_pond, &
2486 : Cdn_atm_rdg , hfreebd , &
2487 : hdraft , hridge , &
2488 : distrdg , hkeel , &
2489 : dkeel , lfloe , &
2490 96528 : dfloe , ncat)
2491 96528 : if (icepack_warnings_aborted(subname)) return
2492 : endif
2493 :
2494 7834848 : do n = 1, ncat
2495 :
2496 6464688 : meltsn (n) = c0
2497 6464688 : melttn (n) = c0
2498 6464688 : meltbn (n) = c0
2499 6464688 : congeln(n) = c0
2500 6464688 : snoicen(n) = c0
2501 6464688 : dsnown (n) = c0
2502 :
2503 6464688 : Trefn = c0
2504 6464688 : Qrefn = c0
2505 7387008 : Qrefn_iso(:) = c0
2506 7387008 : fiso_ocnn(:) = c0
2507 7387008 : fiso_evapn(:) = c0
2508 6464688 : Urefn = c0
2509 6464688 : lhcoef = c0
2510 6464688 : shcoef = c0
2511 6464688 : worka = c0
2512 6464688 : workb = c0
2513 :
2514 6464688 : if (aicen_init(n) > puny) then
2515 :
2516 4360794 : if (calc_Tsfc .or. calc_strair) then
2517 :
2518 : !-----------------------------------------------------------------
2519 : ! Atmosphere boundary layer calculation; compute coefficients
2520 : ! for sensible and latent heat fluxes.
2521 : !
2522 : ! NOTE: The wind stress is computed here for later use if
2523 : ! calc_strair = .true. Otherwise, the wind stress
2524 : ! components are set to the data values.
2525 : !-----------------------------------------------------------------
2526 :
2527 : call icepack_atm_boundary('ice', &
2528 1486692 : Tsfc(n), potT, &
2529 : uatm, vatm, &
2530 : wind, zlvl, &
2531 : Qa, rhoa, &
2532 : strairxn, strairyn, &
2533 : Trefn, Qrefn, &
2534 : worka, workb, &
2535 : lhcoef, shcoef, &
2536 : Cdn_atm, &
2537 : Cdn_atm_ratio_n, &
2538 : Qa_iso=l_Qa_iso, &
2539 0 : Qref_iso=Qrefn_iso, &
2540 : uvel=uvel, vvel=vvel, &
2541 4360794 : Uref=Urefn)
2542 4360794 : if (icepack_warnings_aborted(subname)) return
2543 :
2544 : endif ! calc_Tsfc or calc_strair
2545 :
2546 4360794 : if (.not.(calc_strair)) then
2547 : #ifndef CICE_IN_NEMO
2548 : ! Set to data values (on T points)
2549 0 : strairxn = strax
2550 0 : strairyn = stray
2551 : #else
2552 : ! NEMO wind stress is supplied on u grid, multipied
2553 : ! by ice concentration and set directly in evp, so
2554 : ! strairxT/yT = 0. Zero u-components here for safety.
2555 : strairxn = c0
2556 : strairyn = c0
2557 : #endif
2558 : endif
2559 :
2560 : !-----------------------------------------------------------------
2561 : ! Update ice age
2562 : ! This is further adjusted for freezing in the thermodynamics.
2563 : ! Melting does not alter the ice age.
2564 : !-----------------------------------------------------------------
2565 :
2566 4360794 : if (tr_iage) call increment_age (dt, iage(n))
2567 4360794 : if (icepack_warnings_aborted(subname)) return
2568 4360794 : if (tr_FY) call update_FYarea (dt, &
2569 : lmask_n, lmask_s, &
2570 229760 : yday, FY(n))
2571 4360794 : if (icepack_warnings_aborted(subname)) return
2572 :
2573 : !-----------------------------------------------------------------
2574 : ! Vertical thermodynamics: Heat conduction, growth and melting.
2575 : !-----------------------------------------------------------------
2576 :
2577 4360794 : if (.not.(calc_Tsfc)) then
2578 :
2579 : ! If not calculating surface temperature and fluxes, set
2580 : ! surface fluxes (flatn, fsurfn, and fcondtopn) to be used
2581 : ! in thickness_changes
2582 :
2583 : ! hadgem routine sets fluxes to default values in ice-only mode
2584 0 : call set_sfcflux(aicen (n), &
2585 0 : flatn_f (n), fsensn_f (n), &
2586 0 : fcondtopn_f(n), &
2587 0 : fsurfn_f (n), &
2588 0 : flatn (n), fsensn (n), &
2589 0 : fsurfn (n), &
2590 0 : fcondtopn (n))
2591 0 : if (icepack_warnings_aborted(subname)) return
2592 : endif
2593 :
2594 : call thermo_vertical(nilyr, nslyr, &
2595 1486692 : dt, aicen (n), &
2596 1486692 : vicen (n), vsnon (n), &
2597 1486692 : Tsfc (n), zSin (:,n), &
2598 0 : zqin (:,n), zqsn (:,n), &
2599 1486692 : apnd (n), hpnd (n), &
2600 : tr_pond_topo, &
2601 : flw, potT, &
2602 : Qa, rhoa, &
2603 : fsnow, fpond, &
2604 : fbot, Tbot, &
2605 : Tsnice, sss, &
2606 : lhcoef, shcoef, &
2607 1486692 : fswsfcn (n), fswintn (n), &
2608 0 : Sswabsn(:,n), Iswabsn(:,n), &
2609 1486692 : fsurfn (n), fcondtopn(n), &
2610 1486692 : fcondbotn(n), &
2611 1486692 : fsensn (n), flatn (n), &
2612 : flwoutn, evapn, &
2613 : evapsn, evapin, &
2614 : freshn, fsaltn, &
2615 : fhocnn, &
2616 1486692 : melttn (n), meltsn (n), &
2617 1486692 : meltbn (n), &
2618 1486692 : congeln (n), snoicen (n), &
2619 : mlt_onset, frz_onset, &
2620 1486692 : yday, dsnown (n), &
2621 7334178 : prescribed_ice)
2622 :
2623 4360794 : if (icepack_warnings_aborted(subname)) then
2624 0 : call icepack_warnings_add(subname//' ice: Vertical thermo error: ')
2625 0 : return
2626 : endif
2627 :
2628 : !-----------------------------------------------------------------
2629 : ! Total absorbed shortwave radiation
2630 : !-----------------------------------------------------------------
2631 :
2632 4360794 : fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)
2633 :
2634 : !-----------------------------------------------------------------
2635 : ! Aerosol update
2636 : !-----------------------------------------------------------------
2637 :
2638 4360794 : if (tr_aero) then
2639 : call update_aerosol (dt, &
2640 : nilyr, nslyr, n_aero, &
2641 0 : melttn (n), meltsn (n), &
2642 0 : meltbn (n), congeln (n), &
2643 0 : snoicen (n), fsnow, &
2644 0 : aerosno(:,:,n), aeroice(:,:,n), &
2645 0 : aicen_init (n), vicen_init (n), &
2646 0 : vsnon_init (n), &
2647 0 : vicen (n), vsnon (n), &
2648 0 : aicen (n), &
2649 229765 : faero_atm , faero_ocn)
2650 229765 : if (icepack_warnings_aborted(subname)) return
2651 : endif
2652 :
2653 4360794 : if (tr_iso) then
2654 : call update_isotope (dt = dt, &
2655 : nilyr = nilyr, nslyr = nslyr, &
2656 0 : meltt = melttn(n),melts = meltsn(n), &
2657 0 : meltb = meltbn(n),congel=congeln(n), &
2658 0 : snoice=snoicen(n),evap=evapn, &
2659 0 : fsnow=fsnow, Tsfc=Tsfc(n), &
2660 0 : Qref_iso=Qrefn_iso(:), &
2661 0 : isosno=l_isosno(:,n),isoice=l_isoice(:,n), &
2662 0 : aice_old=aicen_init(n),vice_old=vicen_init(n), &
2663 0 : vsno_old=vsnon_init(n), &
2664 0 : vicen=vicen(n),vsnon=vsnon(n), &
2665 0 : aicen=aicen(n), &
2666 : fiso_atm=l_fiso_atm(:), &
2667 0 : fiso_evapn=fiso_evapn(:), &
2668 0 : fiso_ocnn=fiso_ocnn(:), &
2669 : HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, &
2670 229760 : H2_18O_ocn=l_H2_18O_ocn)
2671 229760 : if (icepack_warnings_aborted(subname)) return
2672 : endif
2673 : endif ! aicen_init
2674 :
2675 : !-----------------------------------------------------------------
2676 : ! Melt ponds
2677 : ! If using tr_pond_cesm, the full calculation is performed here.
2678 : ! If using tr_pond_topo, the rest of the calculation is done after
2679 : ! the surface fluxes are merged, below.
2680 : !-----------------------------------------------------------------
2681 :
2682 : !call ice_timer_start(timer_ponds)
2683 6464688 : if (tr_pond) then
2684 :
2685 5402880 : if (tr_pond_cesm) then
2686 307440 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2687 : call compute_ponds_cesm(dt, hi_min, &
2688 : pndaspect, rfrac, &
2689 0 : melttn(n), meltsn(n), &
2690 : frain, &
2691 0 : aicen (n), vicen (n), &
2692 0 : Tsfc (n), &
2693 307440 : apnd (n), hpnd (n))
2694 307440 : if (icepack_warnings_aborted(subname)) return
2695 :
2696 5095440 : elseif (tr_pond_lvl) then
2697 4305360 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2698 : call compute_ponds_lvl(dt, nilyr, &
2699 : ktherm, &
2700 : hi_min, &
2701 : dpscale, frzpnd, &
2702 : pndaspect, rfrac, &
2703 1670640 : melttn(n), meltsn(n), &
2704 : frain, Tair, &
2705 1670640 : fsurfn(n), &
2706 1670640 : dhsn (n), ffracn(n), &
2707 1670640 : aicen (n), vicen (n), &
2708 1670640 : vsnon (n), &
2709 0 : zqin(:,n), zSin(:,n), &
2710 1670640 : Tsfc (n), alvl (n), &
2711 1670640 : apnd (n), hpnd (n), &
2712 5976000 : ipnd (n))
2713 4305360 : if (icepack_warnings_aborted(subname)) return
2714 :
2715 790080 : elseif (tr_pond_topo) then
2716 790080 : if (aicen_init(n) > puny) then
2717 :
2718 : ! collect liquid water in ponds
2719 : ! assume salt still runs off
2720 457021 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2721 80554 : pond = rfrac/rhofresh * (melttn(n)*rhoi &
2722 80554 : + meltsn(n)*rhos &
2723 457021 : + frain *dt)
2724 :
2725 : ! if pond does not exist, create new pond over full ice area
2726 : ! otherwise increase pond depth without changing pond area
2727 457021 : if (apnd(n) < puny) then
2728 272114 : hpnd(n) = c0
2729 272114 : apnd(n) = c1
2730 : endif
2731 457021 : hpnd(n) = (pond + hpnd(n)*apnd(n)) / apnd(n)
2732 457021 : fpond = fpond + pond * aicen(n) ! m
2733 : endif ! aicen_init
2734 : endif
2735 :
2736 : endif ! tr_pond
2737 : !call ice_timer_stop(timer_ponds)
2738 :
2739 : !-----------------------------------------------------------------
2740 : ! Increment area-weighted fluxes.
2741 : !-----------------------------------------------------------------
2742 :
2743 6464688 : if (aicen_init(n) > puny) &
2744 1486692 : call merge_fluxes (aicen=aicen_init(n), &
2745 : flw=flw, &
2746 : strairxn=strairxn, strairyn=strairyn,&
2747 : Cdn_atm_ratio_n=Cdn_atm_ratio_n, &
2748 1486692 : fsurfn=fsurfn(n), fcondtopn=fcondtopn(n),&
2749 1486692 : fcondbotn=fcondbotn(n), &
2750 1486692 : fsensn=fsensn(n), flatn=flatn(n), &
2751 : fswabsn=fswabsn, flwoutn=flwoutn, &
2752 : evapn=evapn, &
2753 : evapsn=evapsn, evapin=evapin, &
2754 : Trefn=Trefn, Qrefn=Qrefn, &
2755 : freshn=freshn, fsaltn=fsaltn, &
2756 : fhocnn=fhocnn, &
2757 1486692 : fswthrun=fswthrun(n), &
2758 0 : fswthrun_vdr=l_fswthrun_vdr(n), &
2759 0 : fswthrun_vdf=l_fswthrun_vdf(n), &
2760 0 : fswthrun_idr=l_fswthrun_idr(n), &
2761 0 : fswthrun_idf=l_fswthrun_idf(n), &
2762 : strairxT=strairxT, strairyT=strairyT,&
2763 : Cdn_atm_ratio=Cdn_atm_ratio, &
2764 : fsurf=fsurf, fcondtop=fcondtop,&
2765 : fcondbot=fcondbot, &
2766 : fsens=fsens, flat=flat, &
2767 : fswabs=fswabs, flwout=flwout, &
2768 : evap=evap, &
2769 : evaps=evaps, evapi=evapi, &
2770 : Tref=Tref, Qref=Qref, &
2771 : fresh=fresh, fsalt=fsalt, &
2772 : fhocn=fhocn, &
2773 : fswthru=fswthru, &
2774 : fswthru_vdr=l_fswthru_vdr, &
2775 : fswthru_vdf=l_fswthru_vdf, &
2776 : fswthru_idr=l_fswthru_idr, &
2777 : fswthru_idf=l_fswthru_idf, &
2778 1486692 : melttn=melttn (n), meltsn=meltsn(n), &
2779 1486692 : meltbn=meltbn (n), congeln=congeln(n),&
2780 1486692 : snoicen=snoicen(n), &
2781 : meltt=meltt, melts=melts, &
2782 : meltb=meltb, congel=congel, &
2783 : snoice=snoice, &
2784 : Uref=Uref, Urefn=Urefn, &
2785 : Qref_iso=l_Qref_iso, &
2786 0 : Qrefn_iso=Qrefn_iso, &
2787 : fiso_ocn=l_fiso_ocn, &
2788 0 : fiso_ocnn=fiso_ocnn, &
2789 : fiso_evap=l_fiso_evap, &
2790 5847486 : fiso_evapn=fiso_evapn)
2791 :
2792 7834848 : if (icepack_warnings_aborted(subname)) return
2793 :
2794 : enddo ! ncat
2795 :
2796 8757168 : if (present(isosno) ) isosno = l_isosno
2797 8757168 : if (present(isoice) ) isoice = l_isoice
2798 5480640 : if (present(Qa_iso) ) Qa_iso = l_Qa_iso
2799 5480640 : if (present(Qref_iso) ) Qref_iso = l_Qref_iso
2800 5480640 : if (present(fiso_atm) ) fiso_atm = l_fiso_atm
2801 5480640 : if (present(fiso_ocn) ) fiso_ocn = l_fiso_ocn
2802 5480640 : if (present(fiso_evap)) fiso_evap= l_fiso_evap
2803 7834848 : if (present(fswthrun_vdr)) fswthrun_vdr= l_fswthrun_vdr
2804 7834848 : if (present(fswthrun_vdf)) fswthrun_vdf= l_fswthrun_vdf
2805 7834848 : if (present(fswthrun_idr)) fswthrun_idr= l_fswthrun_idr
2806 7834848 : if (present(fswthrun_idf)) fswthrun_idf= l_fswthrun_idf
2807 1370160 : if (present(fswthru_vdr)) fswthru_vdr= l_fswthru_vdr
2808 1370160 : if (present(fswthru_vdf)) fswthru_vdf= l_fswthru_vdf
2809 1370160 : if (present(fswthru_idr)) fswthru_idr= l_fswthru_idr
2810 1370160 : if (present(fswthru_idf)) fswthru_idf= l_fswthru_idf
2811 1370160 : deallocate(l_isosno)
2812 1370160 : deallocate(l_isoice)
2813 1370160 : deallocate(l_Qa_iso)
2814 1370160 : deallocate(l_Qref_iso)
2815 1370160 : deallocate(l_fiso_atm)
2816 1370160 : deallocate(l_fiso_ocn)
2817 1370160 : deallocate(l_fiso_evap)
2818 1370160 : deallocate(l_fswthrun_vdr)
2819 1370160 : deallocate(l_fswthrun_vdf)
2820 1370160 : deallocate(l_fswthrun_idr)
2821 1370160 : deallocate(l_fswthrun_idf)
2822 :
2823 : !-----------------------------------------------------------------
2824 : ! Calculate ponds from the topographic scheme
2825 : !-----------------------------------------------------------------
2826 : !call ice_timer_start(timer_ponds)
2827 1370160 : if (tr_pond_topo) then
2828 : call compute_ponds_topo(dt, ncat, nilyr, &
2829 : ktherm, heat_capacity, &
2830 0 : aice, aicen, &
2831 0 : vice, vicen, &
2832 0 : vsno, vsnon, &
2833 : meltt, &
2834 : fsurf, fpond, &
2835 0 : Tsfc, Tf, &
2836 0 : zqin, zSin, &
2837 158016 : apnd, hpnd, ipnd )
2838 158016 : if (icepack_warnings_aborted(subname)) return
2839 : endif
2840 : !call ice_timer_stop(timer_ponds)
2841 :
2842 1370160 : end subroutine icepack_step_therm1
2843 :
2844 : !=======================================================================
2845 :
2846 : end module icepack_therm_vertical
2847 :
2848 : !=======================================================================
|
