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 4017728 : subroutine thermo_vertical (nilyr, nslyr, &
79 : dt, aicen, &
80 : vicen, vsnon, &
81 4017728 : Tsf, zSin, &
82 4017728 : 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 4017728 : 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 1361620 : dhi , & ! change in ice thickness
208 1361620 : dhs ! change in snow thickness
209 :
210 : ! 2D state variables (thickness, temperature)
211 :
212 : real (kind=dbl_kind) :: &
213 1361620 : hilyr , & ! ice layer thickness
214 1361620 : hslyr , & ! snow layer thickness
215 1361620 : hin , & ! ice thickness (m)
216 1361620 : hsn , & ! snow thickness (m)
217 1361620 : hsn_new , & ! thickness of new snow (m)
218 1361620 : worki , & ! local work array
219 1361620 : works ! local work array
220 :
221 : real (kind=dbl_kind), dimension (nilyr) :: &
222 18284700 : zTin ! internal ice layer temperatures
223 :
224 : real (kind=dbl_kind), dimension (nslyr) :: &
225 7063184 : zTsn ! internal snow layer temperatures
226 :
227 : ! other 2D flux and energy variables
228 :
229 : real (kind=dbl_kind) :: &
230 1361620 : einit , & ! initial energy of melting (J m-2)
231 1361620 : efinal , & ! final energy of melting (J m-2)
232 1361620 : einter ! intermediate energy
233 :
234 : real (kind=dbl_kind) :: &
235 1361620 : fadvocn ! advective heat flux to ocean
236 :
237 : character(len=*),parameter :: subname='(thermo_vertical)'
238 :
239 : !-----------------------------------------------------------------
240 : ! Initialize
241 : !-----------------------------------------------------------------
242 :
243 4017728 : flwoutn = c0
244 4017728 : evapn = c0
245 4017728 : evapsn = c0
246 4017728 : evapin = c0
247 4017728 : freshn = c0
248 4017728 : fsaltn = c0
249 4017728 : fhocnn = c0
250 4017728 : fadvocn = c0
251 :
252 4017728 : meltt = c0
253 4017728 : meltb = c0
254 4017728 : melts = c0
255 4017728 : congel = c0
256 4017728 : snoice = c0
257 4017728 : dsnow = c0
258 8944368 : zTsn(:) = c0
259 30322186 : zTin(:) = c0
260 :
261 4017728 : if (calc_Tsfc) then
262 4017728 : fsensn = c0
263 4017728 : flatn = c0
264 4017728 : fsurfn = c0
265 4017728 : 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 4017728 : einit )
281 4017728 : if (icepack_warnings_aborted(subname)) return
282 :
283 : ! Save initial ice and snow thickness (for fresh and fsalt)
284 4017728 : worki = hin
285 4017728 : works = hsn
286 :
287 : !-----------------------------------------------------------------
288 : ! Compute new surface temperature and internal ice and snow
289 : ! temperatures.
290 : !-----------------------------------------------------------------
291 :
292 4017728 : if (heat_capacity) then ! usual case
293 :
294 3714455 : 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 2829466 : fadvocn, snoice)
314 2829466 : 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 4017728 : einter = c0
370 8944368 : do k = 1, nslyr
371 8944368 : einter = einter + hslyr * zqsn(k)
372 : enddo ! k
373 30322186 : do k = 1, nilyr
374 30322186 : einter = einter + hilyr * zqin(k)
375 : enddo ! k
376 :
377 4017728 : Tsnice = c0
378 4017728 : if ((hslyr+hilyr) > puny) then
379 4017728 : if (hslyr > puny) then
380 3490220 : Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr)
381 : else
382 527508 : Tsnice = Tsf
383 : endif
384 : endif
385 :
386 4017728 : 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 4017728 : dsnow)
412 4017728 : 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 4017728 : fadvocn, fbot )
426 4017728 : 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 4017728 : dhi = hin - worki
448 4017728 : dhs = hsn - works - hsn_new
449 :
450 4017728 : freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt
451 4017728 : fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
452 4017728 : fhocnn = fhocnn + fadvocn ! for ktherm=2
453 :
454 4017728 : 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 4017728 : vicen, vsnon)
470 4017728 : 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 1273632 : subroutine frzmlt_bottom_lateral (dt, ncat, &
488 : nilyr, nslyr, &
489 : aice, frzmlt, &
490 2547264 : vicen, vsnon, &
491 2547264 : qicen, qsnon, &
492 : sst, Tf, &
493 : ustar_min, &
494 439248 : 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 439248 : etot , & ! total energy in column
543 439248 : qavg ! average enthalpy in column (approximate)
544 :
545 : real (kind=dbl_kind) :: &
546 439248 : deltaT , & ! SST - Tbot >= 0
547 439248 : ustar , & ! skin friction velocity for fbot (m/s)
548 439248 : wlat , & ! lateral melt rate (m/s)
549 439248 : xtmp ! temporary variable
550 :
551 : ! Parameters for bottom melting
552 :
553 : real (kind=dbl_kind) :: &
554 439248 : 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 1273632 : rside = c0
573 1273632 : fside = c0
574 1273632 : Tbot = Tf
575 1273632 : fbot = c0
576 1273632 : wlat = c0
577 :
578 1273632 : 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 477501 : deltaT = max((sst-Tbot),c0)
586 :
587 : ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
588 477501 : ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow)
589 477501 : ustar = max (ustar,ustar_min)
590 :
591 477501 : 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 477501 : cpchr = -cp_ocn*rhow*0.006_dbl_kind
598 : endif
599 :
600 477501 : fbot = cpchr * deltaT * ustar ! < 0
601 477501 : 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 477501 : wlat = m1 * deltaT**m2 ! Maykut & Perovich
613 477501 : rside = wlat*dt*pi/(floeshape*floediam) ! Steele
614 477501 : rside = max(c0,min(rside,c1))
615 :
616 : !-----------------------------------------------------------------
617 : ! Compute heat flux associated with this value of rside.
618 : !-----------------------------------------------------------------
619 :
620 2862878 : do n = 1, ncat
621 :
622 2385377 : etot = c0
623 2385377 : qavg = c0
624 :
625 : ! melting energy/unit area in each column, etot < 0
626 :
627 5735954 : do k = 1, nslyr
628 3350577 : etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
629 5735954 : qavg = qavg + qsnon(k,n)
630 : enddo
631 :
632 17796634 : do k = 1, nilyr
633 15411257 : etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind)
634 17796634 : qavg = qavg + qicen(k,n)
635 : enddo ! nilyr
636 :
637 : ! lateral heat flux, fside < 0
638 2862878 : if (tr_fsd) then ! floe size distribution
639 58070 : 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 477501 : xtmp = frzmlt/(fbot + fside + puny)
651 477501 : xtmp = min(xtmp, c1)
652 477501 : fbot = fbot * xtmp
653 477501 : rside = rside * xtmp
654 477501 : fside = fside * xtmp
655 :
656 : endif
657 :
658 1273632 : 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 4017728 : subroutine init_vertical_profile(nilyr, nslyr, &
670 : aicen, vicen, &
671 : vsnon, &
672 : hin, hilyr, &
673 : hsn, hslyr, &
674 4017728 : zqin, zTin, &
675 4017728 : zqsn, zTsn, &
676 4017728 : 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 14266972 : Tmlts ! melting temperature
712 :
713 : integer (kind=int_kind) :: &
714 : k ! ice layer index
715 :
716 : real (kind=dbl_kind) :: &
717 1361620 : rnslyr, & ! real(nslyr)
718 1361620 : 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 4017728 : rnslyr = real(nslyr,kind=dbl_kind)
733 :
734 4017728 : tsno_high = .false.
735 4017728 : tice_high = .false.
736 4017728 : tsno_low = .false.
737 4017728 : tice_low = .false.
738 :
739 4017728 : einit = c0
740 :
741 : !-----------------------------------------------------------------
742 : ! Surface temperature, ice and snow thickness
743 : ! Initialize internal energy
744 : !-----------------------------------------------------------------
745 :
746 4017728 : hin = vicen / aicen
747 4017728 : hsn = vsnon / aicen
748 4017728 : hilyr = hin / real(nilyr,kind=dbl_kind)
749 4017728 : 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 8944368 : 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 4926640 : if (hslyr > hs_min/rnslyr .and. heat_capacity) then
767 : ! zqsn < 0
768 959838 : Tmax = -zqsn(k)*puny*rnslyr / &
769 2895687 : (rhos*cp_ice*vsnon)
770 : else
771 2030953 : zqsn (k) = -rhos * Lfresh
772 2030953 : Tmax = puny
773 : endif
774 :
775 : !-----------------------------------------------------------------
776 : ! Compute snow temperatures from enthalpies.
777 : ! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0.
778 : !-----------------------------------------------------------------
779 4926640 : zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice
780 :
781 : !-----------------------------------------------------------------
782 : ! Check for zTsn > Tmax (allowing for roundoff error) and zTsn < Tmin.
783 : !-----------------------------------------------------------------
784 8944368 : if (zTsn(k) > Tmax) then
785 0 : tsno_high = .true.
786 4926640 : 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 4017728 : 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 4017728 : 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 8944368 : do k = 1, nslyr
852 :
853 4926640 : if (zTsn(k) > c0) then ! correct roundoff error
854 8265 : zTsn(k) = c0
855 8265 : zqsn(k) = -rhos*Lfresh
856 : endif
857 :
858 : !-----------------------------------------------------------------
859 : ! initial energy per unit area of ice/snow, relative to 0 C
860 : !-----------------------------------------------------------------
861 8944368 : einit = einit + hslyr*zqsn(k)
862 :
863 : enddo ! nslyr
864 :
865 30322186 : do k = 1, nilyr
866 :
867 : !---------------------------------------------------------------------
868 : ! Use initial salinity profile for thin ice
869 : !---------------------------------------------------------------------
870 :
871 26304458 : 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 26304458 : if (ktherm == 2) then
886 19806262 : 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 26304458 : if (ktherm == 2) then
902 19806262 : 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 26304458 : if (l_brine) then
908 26001185 : 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 26304458 : if (zTin(k) > Tmax) then
917 0 : tice_high = .true.
918 26304458 : 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 26304458 : 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 26304458 : 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 26304458 : 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 30322186 : 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 4017728 : subroutine thickness_changes (nilyr, nslyr, &
1028 : dt, yday, &
1029 : efinal, &
1030 : hin, hilyr, &
1031 : hsn, hslyr, &
1032 4017728 : 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 4017728 : 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 1361620 : esub , & ! energy for sublimation, > 0 (J m-2)
1115 1361620 : econ , & ! energy for condensation, < 0 (J m-2)
1116 1361620 : etop_mlt , & ! energy for top melting, > 0 (J m-2)
1117 1361620 : ebot_mlt , & ! energy for bottom melting, > 0 (J m-2)
1118 1361620 : ebot_gro , & ! energy for bottom growth, < 0 (J m-2)
1119 1361620 : emlt_atm , & ! total energy of brine, from atmosphere (J m-2)
1120 1361620 : emlt_ocn ! total energy of brine, to ocean (J m-2)
1121 :
1122 : real (kind=dbl_kind) :: &
1123 1361620 : qbotmax , & ! max enthalpy of ice growing at bottom
1124 1361620 : dhi , & ! change in ice thickness
1125 1361620 : dhs , & ! change in snow thickness
1126 1361620 : Ti , & ! ice temperature
1127 1361620 : Ts , & ! snow temperature
1128 1361620 : qbot , & ! enthalpy of ice growing at bottom surface (J m-3)
1129 1361620 : qsub , & ! energy/unit volume to sublimate ice/snow (J m-3)
1130 1361620 : hqtot , & ! sum of h*q for two layers
1131 1361620 : wk1 , & ! temporary variable
1132 1361620 : zqsnew , & ! enthalpy of new snow (J m-3)
1133 1361620 : hstot , & ! snow thickness including new snow (m)
1134 1361620 : Tmlts ! melting temperature
1135 :
1136 : real (kind=dbl_kind), dimension (nilyr+1) :: &
1137 16923080 : zi1 , & ! depth of ice layer boundaries (m)
1138 16923080 : zi2 ! adjusted depths, with equal hilyr (m)
1139 :
1140 : real (kind=dbl_kind), dimension (nslyr+1) :: &
1141 9719292 : zs1 , & ! depth of snow layer boundaries (m)
1142 9719292 : zs2 ! adjusted depths, with equal hslyr (m)
1143 :
1144 : real (kind=dbl_kind), dimension (nilyr) :: &
1145 14266972 : dzi ! ice layer thickness after growth/melting
1146 :
1147 : real (kind=dbl_kind), dimension (nslyr) :: &
1148 8357672 : dzs ! snow layer thickness after growth/melting
1149 :
1150 : real (kind=dbl_kind), dimension (nilyr) :: &
1151 15561460 : qm , & ! energy of melting (J m-3) = zqin in BL99 formulation
1152 15628592 : qmlt ! enthalpy of melted ice (J m-3) = zero in BL99 formulation
1153 :
1154 : real (kind=dbl_kind) :: &
1155 1361620 : qbotm , &
1156 1361620 : qbotp , &
1157 1361620 : qbot0
1158 :
1159 : character(len=*),parameter :: subname='(thickness_changes)'
1160 :
1161 : !-----------------------------------------------------------------
1162 : ! Initialize
1163 : !-----------------------------------------------------------------
1164 :
1165 4017728 : dhi = c0
1166 4017728 : dhs = c0
1167 4017728 : hsn_new = c0
1168 :
1169 30322186 : do k = 1, nilyr
1170 30322186 : dzi(k) = hilyr
1171 : enddo
1172 :
1173 8944368 : do k = 1, nslyr
1174 8944368 : dzs(k) = hslyr
1175 : enddo
1176 :
1177 30322186 : do k = 1, nilyr
1178 26304458 : if (ktherm == 2) then
1179 19806262 : qmlt(k) = enthalpy_of_melting(zSin(k))
1180 : else
1181 6498196 : qmlt(k) = c0
1182 : endif
1183 26304458 : qm(k) = zqin(k) - qmlt(k)
1184 26304458 : emlt_atm = c0
1185 30322186 : 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 4017728 : 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 4017728 : wk1 = -flatn * dt
1222 4017728 : esub = max(wk1, c0) ! energy for sublimation, > 0
1223 4017728 : econ = min(wk1, c0) ! energy for condensation, < 0
1224 :
1225 4017728 : wk1 = (fsurfn - fcondtopn) * dt
1226 4017728 : etop_mlt = max(wk1, c0) ! etop_mlt > 0
1227 :
1228 4017728 : wk1 = (fcondbotn - fbot) * dt
1229 4017728 : ebot_mlt = max(wk1, c0) ! ebot_mlt > 0
1230 4017728 : 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 4017728 : evapn = c0 ! initialize
1240 4017728 : evapsn = c0 ! initialize
1241 4017728 : evapin = c0 ! initialize
1242 :
1243 4017728 : if (hsn > puny) then ! add snow with enthalpy zqsn(1)
1244 3490220 : dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0
1245 3490220 : dzs(1) = dzs(1) + dhs
1246 3490220 : evapn = evapn + dhs*rhos
1247 3490220 : evapsn = evapsn + dhs*rhos
1248 : else ! add ice with enthalpy zqin(1)
1249 527508 : dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0
1250 527508 : dzi(1) = dzi(1) + dhi
1251 527508 : evapn = evapn + dhi*rhoi
1252 527508 : evapin = evapin + dhi*rhoi
1253 : ! enthalpy of melt water
1254 527508 : emlt_atm = emlt_atm - qmlt(1) * dhi
1255 : endif
1256 :
1257 : !--------------------------------------------------------------
1258 : ! Grow ice (bottom)
1259 : !--------------------------------------------------------------
1260 :
1261 4017728 : if (ktherm == 2) then
1262 :
1263 2829466 : qbotm = enthalpy_mush(Tbot, sss)
1264 2829466 : qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
1265 2829466 : qbot0 = qbotm - qbotp
1266 :
1267 2829466 : dhi = ebot_gro / qbotp ! dhi > 0
1268 :
1269 2829466 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
1270 2829466 : hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
1271 2829466 : 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 4017728 : dzi(nilyr) = dzi(nilyr) + dhi
1299 4017728 : if (dzi(nilyr) > puny) then
1300 4017728 : zqin(nilyr) = hqtot / dzi(nilyr)
1301 4017728 : if (ktherm == 2) then
1302 2829466 : zSin(nilyr) = hstot / dzi(nilyr)
1303 2829466 : qmlt(nilyr) = enthalpy_of_melting(zSin(nilyr))
1304 : else
1305 1188262 : qmlt(nilyr) = c0
1306 : endif
1307 4017728 : 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 4017728 : congel = congel + dhi
1316 4017728 : if (dhi > puny .and. frz_onset < puny) &
1317 87 : frz_onset = yday
1318 :
1319 8944368 : do k = 1, nslyr
1320 :
1321 : !--------------------------------------------------------------
1322 : ! Remove internal snow melt
1323 : !--------------------------------------------------------------
1324 :
1325 4926640 : if (ktherm == 2 .and. zqsn(k) > -rhos * Lfresh) then
1326 :
1327 20124 : dhs = max(-dzs(k), &
1328 50394 : -((zqsn(k) + rhos*Lfresh) / (rhos*Lfresh)) * dzs(k))
1329 50394 : dzs(k) = dzs(k) + dhs
1330 50394 : zqsn(k) = -rhos * Lfresh
1331 50394 : melts = melts - dhs
1332 : ! delta E = zqsn(k) + rhos * Lfresh
1333 :
1334 : endif
1335 :
1336 : !--------------------------------------------------------------
1337 : ! Sublimation of snow (evapn < 0)
1338 : !--------------------------------------------------------------
1339 :
1340 4926640 : qsub = zqsn(k) - rhos*Lvap ! qsub < 0
1341 4926640 : dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0
1342 4926640 : dzs(k) = dzs(k) + dhs
1343 4926640 : esub = esub - dhs*qsub
1344 4926640 : esub = max(esub, c0) ! in case of roundoff error
1345 4926640 : evapn = evapn + dhs*rhos
1346 4926640 : evapsn = evapsn + dhs*rhos
1347 :
1348 : !--------------------------------------------------------------
1349 : ! Melt snow (top)
1350 : !--------------------------------------------------------------
1351 :
1352 4926640 : dhs = max(-dzs(k), etop_mlt/zqsn(k))
1353 4926640 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1354 4926640 : etop_mlt = etop_mlt - dhs*zqsn(k)
1355 4926640 : etop_mlt = max(etop_mlt, c0) ! in case of roundoff error
1356 :
1357 : ! history diagnostics
1358 4926640 : if (dhs < -puny .and. mlt_onset < puny) &
1359 22 : mlt_onset = yday
1360 8944368 : melts = melts - dhs
1361 :
1362 : enddo ! nslyr
1363 :
1364 30322186 : do k = 1, nilyr
1365 :
1366 : !--------------------------------------------------------------
1367 : ! Sublimation of ice (evapn < 0)
1368 : !--------------------------------------------------------------
1369 :
1370 26304458 : qsub = qm(k) - rhoi*Lvap ! qsub < 0
1371 26304458 : dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0
1372 26304458 : dzi(k) = dzi(k) + dhi
1373 26304458 : esub = esub - dhi*qsub
1374 26304458 : esub = max(esub, c0)
1375 26304458 : evapn = evapn + dhi*rhoi
1376 26304458 : evapin = evapin + dhi*rhoi
1377 26304458 : emlt_ocn = emlt_ocn - qmlt(k) * dhi
1378 :
1379 : !--------------------------------------------------------------
1380 : ! Melt ice (top)
1381 : !--------------------------------------------------------------
1382 :
1383 26304458 : if (qm(k) < c0) then
1384 26304332 : dhi = max(-dzi(k), etop_mlt/qm(k))
1385 : else
1386 126 : qm(k) = c0
1387 126 : dhi = -dzi(k)
1388 : endif
1389 26304458 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1390 :
1391 26304458 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1392 26304458 : etop_mlt = max(etop_mlt - dhi*qm(k), c0)
1393 :
1394 : ! history diagnostics
1395 26304458 : if (dhi < -puny .and. mlt_onset < puny) &
1396 83 : mlt_onset = yday
1397 30322186 : meltt = meltt - dhi
1398 :
1399 : enddo ! nilyr
1400 :
1401 30322186 : do k = nilyr, 1, -1
1402 :
1403 : !--------------------------------------------------------------
1404 : ! Melt ice (bottom)
1405 : !--------------------------------------------------------------
1406 :
1407 26304458 : if (qm(k) < c0) then
1408 26304332 : dhi = max(-dzi(k), ebot_mlt/qm(k))
1409 : else
1410 126 : qm(k) = c0
1411 126 : dhi = -dzi(k)
1412 : endif
1413 26304458 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1414 :
1415 26304458 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1416 26304458 : ebot_mlt = max(ebot_mlt - dhi*qm(k), c0)
1417 :
1418 : ! history diagnostics
1419 30322186 : meltb = meltb -dhi
1420 :
1421 : enddo ! nilyr
1422 :
1423 8944368 : do k = nslyr, 1, -1
1424 :
1425 : !--------------------------------------------------------------
1426 : ! Melt snow (only if all the ice has melted)
1427 : !--------------------------------------------------------------
1428 :
1429 4926640 : dhs = max(-dzs(k), ebot_mlt/zqsn(k))
1430 4926640 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1431 4926640 : ebot_mlt = ebot_mlt - dhs*zqsn(k)
1432 8944368 : 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 4017728 : + (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 4017728 : if (fsnow > c0) then
1455 :
1456 3153688 : hsn_new = fsnow/rhos * dt
1457 3153688 : zqsnew = -rhos*Lfresh
1458 3153688 : hstot = dzs(1) + hsn_new
1459 :
1460 3153688 : if (hstot > c0) then
1461 1032216 : zqsn(1) = (dzs(1) * zqsn(1) &
1462 3153688 : + hsn_new * zqsnew) / hstot
1463 : ! avoid roundoff errors
1464 3153688 : zqsn(1) = min(zqsn(1), -rhos*Lfresh)
1465 :
1466 3153688 : dzs(1) = hstot
1467 : endif
1468 : endif
1469 :
1470 : !-----------------------------------------------------------------
1471 : ! Find the new ice and snow thicknesses.
1472 : !-----------------------------------------------------------------
1473 :
1474 4017728 : hin = c0
1475 4017728 : hsn = c0
1476 :
1477 30322186 : do k = 1, nilyr
1478 30322186 : hin = hin + dzi(k)
1479 : enddo ! k
1480 :
1481 8944368 : do k = 1, nslyr
1482 4926640 : hsn = hsn + dzs(k)
1483 8944368 : dsnow = dsnow + dzs(k) - hslyr
1484 : enddo ! k
1485 :
1486 : !-------------------------------------------------------------------
1487 : ! Convert snow to ice if snow lies below freeboard.
1488 : !-------------------------------------------------------------------
1489 :
1490 4017728 : if (ktherm /= 2) &
1491 : call freeboard (nslyr, &
1492 : snoice, &
1493 : hin, hsn, &
1494 0 : zqin, zqsn, &
1495 0 : dzi, dzs, &
1496 1188262 : dsnow)
1497 4017728 : 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 4017728 : if (hin > c0) then
1509 4017721 : hilyr = hin / real(nilyr,kind=dbl_kind)
1510 : else
1511 7 : hin = c0
1512 7 : hilyr = c0
1513 : endif
1514 4017728 : if (hsn > c0) then
1515 3471896 : hslyr = hsn / real(nslyr,kind=dbl_kind)
1516 : else
1517 545832 : hsn = c0
1518 545832 : 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 4017728 : zi1(1) = c0
1527 4017728 : zi1(1+nilyr) = hin
1528 :
1529 4017728 : zi2(1) = c0
1530 4017728 : zi2(1+nilyr) = hin
1531 :
1532 4017728 : if (heat_capacity) then
1533 :
1534 26001185 : do k = 1, nilyr-1
1535 22286730 : zi1(k+1) = zi1(k) + dzi(k)
1536 26001185 : 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 3714455 : zqin)
1547 3714455 : if (icepack_warnings_aborted(subname)) return
1548 :
1549 3714455 : if (ktherm == 2) &
1550 : call adjust_enthalpy (nilyr, &
1551 0 : zi1, zi2, &
1552 : hilyr, hin, &
1553 2829466 : zSin)
1554 3714455 : 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 4017728 : 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 4017728 : if (ktherm == 2) then
1598 5658932 : do k = 1, nslyr
1599 5658932 : if (hsn <= puny) then
1600 : fhocnn = fhocnn &
1601 68595 : + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
1602 68595 : zqsn(k) = -rhos*Lfresh
1603 68595 : 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 4017728 : efinal = -evapn*Lvap
1614 4017728 : evapn = evapn/dt
1615 4017728 : evapsn = evapsn/dt
1616 4017728 : evapin = evapin/dt
1617 :
1618 8944368 : do k = 1, nslyr
1619 8944368 : efinal = efinal + hslyr*zqsn(k)
1620 : enddo
1621 :
1622 30322186 : do k = 1, nilyr
1623 30322186 : efinal = efinal + hilyr*zqin(k)
1624 : enddo ! k
1625 :
1626 4017728 : 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 4017728 : fhocnn = fhocnn + emlt_ocn/dt
1633 4017728 : 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 6771149 : subroutine adjust_enthalpy (nlyr, &
1755 6771149 : z1, z2, &
1756 : hlyr, hn, &
1757 6771149 : 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 2274002 : hovlp ! overlap between old and new layers (m)
1782 :
1783 : real (kind=dbl_kind) :: &
1784 2274002 : rhlyr ! 1./hlyr
1785 :
1786 : real (kind=dbl_kind), dimension (nlyr) :: &
1787 27025202 : hq ! h * q for a layer
1788 :
1789 : character(len=*),parameter :: subname='(adjust_enthalpy)'
1790 :
1791 : !-----------------------------------------------------------------
1792 : ! Compute reciprocal layer thickness.
1793 : !-----------------------------------------------------------------
1794 :
1795 6771149 : rhlyr = c0
1796 6771149 : 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 53714736 : do k2 = 1, nlyr
1803 53714736 : hq(k2) = c0
1804 : enddo ! k
1805 6771149 : k1 = 1
1806 6771149 : k2 = 1
1807 93105986 : do while (k1 <= nlyr .and. k2 <= nlyr)
1808 57962948 : hovlp = min (z1(k1+1), z2(k2+1)) &
1809 144297785 : - max (z1(k1), z2(k2))
1810 86334837 : hovlp = max (hovlp, c0)
1811 86334837 : hq(k2) = hq(k2) + hovlp*qn(k1)
1812 93105986 : if (z1(k1+1) > z2(k2+1)) then
1813 39391250 : k2 = k2 + 1
1814 : else
1815 46943587 : k1 = k1 + 1
1816 : endif
1817 : enddo ! while
1818 :
1819 : !-----------------------------------------------------------------
1820 : ! Compute new enthalpies.
1821 : !-----------------------------------------------------------------
1822 :
1823 53714736 : do k = 1, nlyr
1824 53714736 : qn(k) = hq(k) * rhlyr
1825 : enddo ! k
1826 :
1827 6771149 : 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 4017728 : 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 1361620 : einp , & ! energy input during timestep (J m-2)
1870 1361620 : 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 4017728 : - fsnow*Lfresh - fadvocn) * dt
1887 4017728 : ferr = abs(efinal-einit-einp) / dt
1888 :
1889 4017728 : 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 4017728 : subroutine update_state_vthermo(nilyr, nslyr, &
1960 : Tf, Tsf, &
1961 : hin, hsn, &
1962 4017728 : zqin, zSin, &
1963 4017728 : 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 4017728 : 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 4017721 : vicen = aicen * hin
2017 4017721 : vsnon = aicen * hsn
2018 : endif
2019 :
2020 4017728 : 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 1273632 : subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
2031 1273632 : aicen_init , &
2032 1273632 : vicen_init , vsnon_init , &
2033 1273632 : aice , aicen , &
2034 1273632 : vice , vicen , &
2035 1273632 : vsno , vsnon , &
2036 : uvel , vvel , &
2037 2547264 : Tsfc , zqsn , &
2038 2547264 : zqin , zSin , &
2039 1273632 : alvl , vlvl , &
2040 1273632 : apnd , hpnd , &
2041 1273632 : ipnd , &
2042 1273632 : iage , FY , &
2043 1273632 : aerosno , aeroice , &
2044 1273632 : isosno , isoice , &
2045 : uatm , vatm , &
2046 : wind , zlvl , &
2047 : Qa , rhoa , &
2048 1273632 : Qa_iso , &
2049 : Tair , Tref , &
2050 : Qref , Uref , &
2051 1273632 : 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 1273632 : fsurf , fsurfn , &
2074 1273632 : fcondtop , fcondtopn , &
2075 1273632 : fcondbot , fcondbotn , &
2076 1273632 : fswsfcn , fswintn , &
2077 1273632 : fswthrun , fswabs , &
2078 : flwout , &
2079 1273632 : Sswabsn , Iswabsn , &
2080 : flw , &
2081 1273632 : fsens , fsensn , &
2082 1273632 : flat , flatn , &
2083 : evap , &
2084 : evaps , evapi , &
2085 : fresh , fsalt , &
2086 : fhocn , fswthru , &
2087 1273632 : flatn_f , fsensn_f , &
2088 1273632 : fsurfn_f , fcondtopn_f , &
2089 1273632 : faero_atm , faero_ocn , &
2090 1273632 : fiso_atm , fiso_ocn , &
2091 1273632 : fiso_evap , &
2092 : HDO_ocn , H2_16O_ocn , &
2093 : H2_18O_ocn , &
2094 1273632 : dhsn , ffracn , &
2095 1273632 : meltt , melttn , &
2096 1273632 : meltb , meltbn , &
2097 1273632 : melts , meltsn , &
2098 1273632 : congel , congeln , &
2099 1273632 : snoice , snoicen , &
2100 1273632 : dsnown , &
2101 : lmask_n , lmask_s , &
2102 : mlt_onset , frz_onset , &
2103 : yday , prescribed_ice)
2104 :
2105 : integer (kind=int_kind), intent(in) :: &
2106 : ncat , & ! number of thickness categories
2107 : nilyr , & ! number of ice layers
2108 : nslyr ! number of snow layers
2109 :
2110 : real (kind=dbl_kind), intent(in) :: &
2111 : dt , & ! time step
2112 : uvel , & ! x-component of velocity (m/s)
2113 : vvel , & ! y-component of velocity (m/s)
2114 : strax , & ! wind stress components (N/m^2)
2115 : stray , & !
2116 : yday ! day of year
2117 :
2118 : logical (kind=log_kind), intent(in) :: &
2119 : lmask_n , & ! northern hemisphere mask
2120 : lmask_s ! southern hemisphere mask
2121 :
2122 : logical (kind=log_kind), intent(in), optional :: &
2123 : prescribed_ice ! if .true., use prescribed ice instead of computed
2124 :
2125 : real (kind=dbl_kind), intent(inout) :: &
2126 : aice , & ! sea ice concentration
2127 : vice , & ! volume per unit area of ice (m)
2128 : vsno , & ! volume per unit area of snow (m)
2129 : zlvl , & ! atm level height (m)
2130 : uatm , & ! wind velocity components (m/s)
2131 : vatm , &
2132 : wind , & ! wind speed (m/s)
2133 : potT , & ! air potential temperature (K)
2134 : Tair , & ! air temperature (K)
2135 : Qa , & ! specific humidity (kg/kg)
2136 : rhoa , & ! air density (kg/m^3)
2137 : frain , & ! rainfall rate (kg/m^2 s)
2138 : fsnow , & ! snowfall rate (kg/m^2 s)
2139 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
2140 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
2141 : fsalt , & ! salt flux to ocean (kg/m^2/s)
2142 : fhocn , & ! net heat flux to ocean (W/m^2)
2143 : fswthru , & ! shortwave penetrating to ocean (W/m^2)
2144 : fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2)
2145 : fcondtop , & ! top surface conductive flux (W/m^2)
2146 : fcondbot , & ! bottom surface conductive flux (W/m^2)
2147 : fsens , & ! sensible heat flux (W/m^2)
2148 : flat , & ! latent heat flux (W/m^2)
2149 : fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2)
2150 : flw , & ! incoming longwave radiation (W/m^2)
2151 : flwout , & ! outgoing longwave radiation (W/m^2)
2152 : evap , & ! evaporative water flux (kg/m^2/s)
2153 : evaps , & ! evaporative water flux over snow (kg/m^2/s)
2154 : evapi , & ! evaporative water flux over ice (kg/m^2/s)
2155 : congel , & ! basal ice growth (m/step-->cm/day)
2156 : snoice , & ! snow-ice formation (m/step-->cm/day)
2157 : Tref , & ! 2m atm reference temperature (K)
2158 : Qref , & ! 2m atm reference spec humidity (kg/kg)
2159 : Uref , & ! 10m atm reference wind speed (m/s)
2160 : Cdn_atm , & ! atm drag coefficient
2161 : Cdn_ocn , & ! ocn drag coefficient
2162 : hfreebd , & ! freeboard (m)
2163 : hdraft , & ! draft of ice + snow column (Stoessel1993)
2164 : hridge , & ! ridge height
2165 : distrdg , & ! distance between ridges
2166 : hkeel , & ! keel depth
2167 : dkeel , & ! distance between keels
2168 : lfloe , & ! floe length
2169 : dfloe , & ! distance between floes
2170 : Cdn_atm_skin, & ! neutral skin drag coefficient
2171 : Cdn_atm_floe, & ! neutral floe edge drag coefficient
2172 : Cdn_atm_pond, & ! neutral pond edge drag coefficient
2173 : Cdn_atm_rdg , & ! neutral ridge drag coefficient
2174 : Cdn_ocn_skin, & ! skin drag coefficient
2175 : Cdn_ocn_floe, & ! floe edge drag coefficient
2176 : Cdn_ocn_keel, & ! keel drag coefficient
2177 : Cdn_atm_ratio,& ! ratio drag atm / neutral drag atm
2178 : strairxT , & ! stress on ice by air, x-direction
2179 : strairyT , & ! stress on ice by air, y-direction
2180 : strocnxT , & ! ice-ocean stress, x-direction
2181 : strocnyT , & ! ice-ocean stress, y-direction
2182 : fbot , & ! ice-ocean heat flux at bottom surface (W/m^2)
2183 : frzmlt , & ! freezing/melting potential (W/m^2)
2184 : rside , & ! fraction of ice that melts laterally
2185 : fside , & ! lateral heat flux (W/m^2)
2186 : sst , & ! sea surface temperature (C)
2187 : Tf , & ! freezing temperature (C)
2188 : Tbot , & ! ice bottom surface temperature (deg C)
2189 : Tsnice , & ! snow ice interface temperature (deg C)
2190 : sss , & ! sea surface salinity (ppt)
2191 : meltt , & ! top ice melt (m/step-->cm/day)
2192 : melts , & ! snow melt (m/step-->cm/day)
2193 : meltb , & ! basal ice melt (m/step-->cm/day)
2194 : mlt_onset , & ! day of year that sfc melting begins
2195 : frz_onset ! day of year that freezing begins (congel or frazil)
2196 :
2197 : real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
2198 : Qa_iso , & ! isotope specific humidity (kg/kg)
2199 : Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg)
2200 : fiso_atm , & ! isotope deposition rate (kg/m^2 s)
2201 : fiso_ocn , & ! isotope flux to ocean (kg/m^2/s)
2202 : fiso_evap ! isotope evaporation (kg/m^2/s)
2203 :
2204 : real (kind=dbl_kind), optional, intent(in) :: &
2205 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
2206 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
2207 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
2208 :
2209 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
2210 : aicen_init , & ! fractional area of ice
2211 : vicen_init , & ! volume per unit area of ice (m)
2212 : vsnon_init , & ! volume per unit area of snow (m)
2213 : aicen , & ! concentration of ice
2214 : vicen , & ! volume per unit area of ice (m)
2215 : vsnon , & ! volume per unit area of snow (m)
2216 : Tsfc , & ! ice/snow surface temperature, Tsfcn
2217 : alvl , & ! level ice area fraction
2218 : vlvl , & ! level ice volume fraction
2219 : apnd , & ! melt pond area fraction
2220 : hpnd , & ! melt pond depth (m)
2221 : ipnd , & ! melt pond refrozen lid thickness (m)
2222 : iage , & ! volume-weighted ice age
2223 : FY , & ! area-weighted first-year ice area
2224 : fsurfn , & ! net flux to top surface, excluding fcondtop
2225 : fcondtopn , & ! downward cond flux at top surface (W m-2)
2226 : fcondbotn , & ! downward cond flux at bottom surface (W m-2)
2227 : flatn , & ! latent heat flux (W m-2)
2228 : fsensn , & ! sensible heat flux (W m-2)
2229 : fsurfn_f , & ! net flux to top surface, excluding fcondtop
2230 : fcondtopn_f , & ! downward cond flux at top surface (W m-2)
2231 : flatn_f , & ! latent heat flux (W m-2)
2232 : fsensn_f , & ! sensible heat flux (W m-2)
2233 : fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
2234 : fswthrun , & ! SW through ice to ocean (W/m^2)
2235 : fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
2236 : faero_atm , & ! aerosol deposition rate (kg/m^2 s)
2237 : faero_ocn , & ! aerosol flux to ocean (kg/m^2/s)
2238 : dhsn , & ! depth difference for snow on sea ice and pond ice
2239 : ffracn , & ! fraction of fsurfn used to melt ipond
2240 : meltsn , & ! snow melt (m)
2241 : melttn , & ! top ice melt (m)
2242 : meltbn , & ! bottom ice melt (m)
2243 : congeln , & ! congelation ice growth (m)
2244 : snoicen , & ! snow-ice growth (m)
2245 : dsnown ! change in snow thickness (m/step-->cm/day)
2246 :
2247 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
2248 : zqsn , & ! snow layer enthalpy (J m-3)
2249 : zqin , & ! ice layer enthalpy (J m-3)
2250 : zSin , & ! internal ice layer salinities
2251 : Sswabsn , & ! SW radiation absorbed in snow layers (W m-2)
2252 : Iswabsn ! SW radiation absorbed in ice layers (W m-2)
2253 :
2254 : real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: &
2255 : aerosno , & ! snow aerosol tracer (kg/m^2)
2256 : aeroice ! ice aerosol tracer (kg/m^2)
2257 :
2258 : real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: &
2259 : isosno , & ! snow isotope tracer (kg/m^2)
2260 : isoice ! ice isotope tracer (kg/m^2)
2261 : !autodocument_end
2262 :
2263 : ! local variables
2264 :
2265 : integer (kind=int_kind) :: &
2266 : n ! category index
2267 :
2268 : real (kind=dbl_kind) :: &
2269 439248 : worka, workb ! temporary variables
2270 :
2271 : ! 2D coupler variables (computed for each category, then aggregated)
2272 : real (kind=dbl_kind) :: &
2273 439248 : fswabsn , & ! shortwave absorbed by ice (W/m^2)
2274 439248 : flwoutn , & ! upward LW at surface (W/m^2)
2275 439248 : evapn , & ! flux of vapor, atmos to ice (kg m-2 s-1)
2276 439248 : evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1)
2277 439248 : evapin , & ! flux of vapor, atmos to ice over ice (kg m-2 s-1)
2278 439248 : freshn , & ! flux of water, ice to ocean (kg/m^2/s)
2279 439248 : fsaltn , & ! flux of salt, ice to ocean (kg/m^2/s)
2280 439248 : fhocnn , & ! fbot corrected for leftover energy (W/m^2)
2281 439248 : strairxn , & ! air/ice zonal stress, (N/m^2)
2282 439248 : strairyn , & ! air/ice meridional stress, (N/m^2)
2283 439248 : Cdn_atm_ratio_n, & ! drag coefficient ratio
2284 439248 : Trefn , & ! air tmp reference level (K)
2285 439248 : Urefn , & ! air speed reference level (m/s)
2286 439248 : Qrefn , & ! air sp hum reference level (kg/kg)
2287 439248 : shcoef , & ! transfer coefficient for sensible heat
2288 439248 : lhcoef , & ! transfer coefficient for latent heat
2289 439248 : rfrac ! water fraction retained for melt ponds
2290 :
2291 : real (kind=dbl_kind), dimension(n_iso) :: &
2292 2108016 : Qrefn_iso , & ! isotope air sp hum reference level (kg/kg)
2293 2108016 : fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s)
2294 2986512 : fiso_evapn ! isotope evaporation (kg/m^2/s)
2295 :
2296 : real (kind=dbl_kind), allocatable, dimension(:,:) :: &
2297 1273632 : l_isosno , & ! local snow isotope tracer (kg/m^2)
2298 1273632 : l_isoice ! local ice isotope tracer (kg/m^2)
2299 :
2300 : real (kind=dbl_kind), allocatable, dimension(:) :: &
2301 1273632 : l_Qa_iso , & ! local isotope specific humidity (kg/kg)
2302 1273632 : l_Qref_iso , & ! local isotope 2m atm reference spec humidity (kg/kg)
2303 1273632 : l_fiso_atm , & ! local isotope deposition rate (kg/m^2 s)
2304 1273632 : l_fiso_ocn , & ! local isotope flux to ocean (kg/m^2/s)
2305 1273632 : l_fiso_evap ! local isotope evaporation (kg/m^2/s)
2306 :
2307 : real (kind=dbl_kind) :: &
2308 439248 : l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg)
2309 439248 : l_H2_16O_ocn, & ! local ocean concentration of H2_16O (kg/kg)
2310 439248 : l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg)
2311 :
2312 : real (kind=dbl_kind) :: &
2313 439248 : pond ! water retained in ponds (m)
2314 :
2315 : character(len=*),parameter :: subname='(icepack_step_therm1)'
2316 :
2317 : !-----------------------------------------------------------------
2318 : ! allocate local optional arguments
2319 : !-----------------------------------------------------------------
2320 :
2321 1273632 : if (present(isosno) ) then
2322 1273632 : allocate(l_isosno(size(isosno,dim=1),size(isosno,dim=2)))
2323 8178000 : l_isosno = isosno
2324 : else
2325 0 : allocate(l_isosno(1,1))
2326 0 : l_isosno = c0
2327 : endif
2328 :
2329 1273632 : if (present(isoice) ) then
2330 1273632 : allocate(l_isoice(size(isoice,dim=1),size(isoice,dim=2)))
2331 8178000 : l_isoice = isoice
2332 : else
2333 0 : allocate(l_isoice(1,1))
2334 0 : l_isoice = c0
2335 : endif
2336 :
2337 1273632 : if (present(Qa_iso) ) then
2338 1273632 : allocate(l_Qa_iso(size(Qa_iso)))
2339 5094528 : l_Qa_iso = Qa_iso
2340 : else
2341 0 : allocate(l_Qa_iso(1))
2342 0 : l_Qa_iso = c0
2343 : endif
2344 :
2345 1273632 : if (present(Qref_iso) ) then
2346 1273632 : allocate(l_Qref_iso(size(Qref_iso)))
2347 5094528 : l_Qref_iso = Qref_iso
2348 : else
2349 0 : allocate(l_Qref_iso(1))
2350 0 : l_Qref_iso = c0
2351 : endif
2352 :
2353 1273632 : if (present(fiso_atm) ) then
2354 1273632 : allocate(l_fiso_atm(size(fiso_atm)))
2355 5094528 : l_fiso_atm = fiso_atm
2356 : else
2357 0 : allocate(l_fiso_atm(1))
2358 0 : l_fiso_atm = c0
2359 : endif
2360 :
2361 1273632 : if (present(fiso_ocn) ) then
2362 1273632 : allocate(l_fiso_ocn(size(fiso_ocn)))
2363 5094528 : l_fiso_ocn = fiso_ocn
2364 : else
2365 0 : allocate(l_fiso_ocn(1))
2366 0 : l_fiso_ocn = c0
2367 : endif
2368 :
2369 1273632 : if (present(fiso_evap) ) then
2370 1273632 : allocate(l_fiso_evap(size(fiso_evap)))
2371 5094528 : l_fiso_evap = fiso_evap
2372 : else
2373 0 : allocate(l_fiso_evap(1))
2374 0 : l_fiso_evap = c0
2375 : endif
2376 :
2377 1273632 : l_HDO_ocn = c0
2378 1273632 : if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn
2379 :
2380 1273632 : l_H2_16O_ocn = c0
2381 1273632 : if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
2382 :
2383 1273632 : l_H2_18O_ocn = c0
2384 1273632 : if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
2385 :
2386 : !-----------------------------------------------------------------
2387 : ! Adjust frzmlt to account for ice-ocean heat fluxes since last
2388 : ! call to coupler.
2389 : ! Compute lateral and bottom heat fluxes.
2390 : !-----------------------------------------------------------------
2391 :
2392 : call frzmlt_bottom_lateral (dt, ncat, &
2393 : nilyr, nslyr, &
2394 : aice, frzmlt, &
2395 0 : vicen, vsnon, &
2396 0 : zqin, zqsn, &
2397 : sst, Tf, &
2398 : ustar_min, &
2399 : fbot_xfer_type, &
2400 : strocnxT, strocnyT, &
2401 : Tbot, fbot, &
2402 : rside, Cdn_ocn, &
2403 1273632 : fside)
2404 :
2405 1273632 : if (icepack_warnings_aborted(subname)) return
2406 :
2407 : !-----------------------------------------------------------------
2408 : ! Update the neutral drag coefficients to account for form drag
2409 : ! Oceanic and atmospheric drag coefficients
2410 : !-----------------------------------------------------------------
2411 :
2412 1273632 : if (formdrag) then
2413 0 : call neutral_drag_coeffs (apnd , &
2414 0 : hpnd , ipnd , &
2415 0 : alvl , vlvl , &
2416 : aice , vice, &
2417 0 : vsno , aicen , &
2418 0 : vicen , &
2419 : Cdn_ocn , Cdn_ocn_skin, &
2420 : Cdn_ocn_floe, Cdn_ocn_keel, &
2421 : Cdn_atm , Cdn_atm_skin, &
2422 : Cdn_atm_floe, Cdn_atm_pond, &
2423 : Cdn_atm_rdg , hfreebd , &
2424 : hdraft , hridge , &
2425 : distrdg , hkeel , &
2426 : dkeel , lfloe , &
2427 96528 : dfloe , ncat)
2428 96528 : if (icepack_warnings_aborted(subname)) return
2429 : endif
2430 :
2431 7255680 : do n = 1, ncat
2432 :
2433 5982048 : meltsn (n) = c0
2434 5982048 : melttn (n) = c0
2435 5982048 : meltbn (n) = c0
2436 5982048 : congeln(n) = c0
2437 5982048 : snoicen(n) = c0
2438 5982048 : dsnown (n) = c0
2439 :
2440 5982048 : Trefn = c0
2441 5982048 : Qrefn = c0
2442 6904368 : Qrefn_iso(:) = c0
2443 6904368 : fiso_ocnn(:) = c0
2444 6904368 : fiso_evapn(:) = c0
2445 5982048 : Urefn = c0
2446 5982048 : lhcoef = c0
2447 5982048 : shcoef = c0
2448 5982048 : worka = c0
2449 5982048 : workb = c0
2450 :
2451 5982048 : if (aicen_init(n) > puny) then
2452 :
2453 4017728 : if (calc_Tsfc .or. calc_strair) then
2454 :
2455 : !-----------------------------------------------------------------
2456 : ! Atmosphere boundary layer calculation; compute coefficients
2457 : ! for sensible and latent heat fluxes.
2458 : !
2459 : ! NOTE: The wind stress is computed here for later use if
2460 : ! calc_strair = .true. Otherwise, the wind stress
2461 : ! components are set to the data values.
2462 : !-----------------------------------------------------------------
2463 :
2464 : call icepack_atm_boundary('ice', &
2465 1361620 : Tsfc(n), potT, &
2466 : uatm, vatm, &
2467 : wind, zlvl, &
2468 : Qa, rhoa, &
2469 : strairxn, strairyn, &
2470 : Trefn, Qrefn, &
2471 : worka, workb, &
2472 : lhcoef, shcoef, &
2473 : Cdn_atm, &
2474 : Cdn_atm_ratio_n, &
2475 : Qa_iso=l_Qa_iso, &
2476 0 : Qref_iso=Qrefn_iso, &
2477 : uvel=uvel, vvel=vvel, &
2478 4017728 : Uref=Urefn)
2479 4017728 : if (icepack_warnings_aborted(subname)) return
2480 :
2481 : endif ! calc_Tsfc or calc_strair
2482 :
2483 4017728 : if (.not.(calc_strair)) then
2484 : #ifndef CICE_IN_NEMO
2485 : ! Set to data values (on T points)
2486 0 : strairxn = strax
2487 0 : strairyn = stray
2488 : #else
2489 : ! NEMO wind stress is supplied on u grid, multipied
2490 : ! by ice concentration and set directly in evp, so
2491 : ! strairxT/yT = 0. Zero u-components here for safety.
2492 : strairxn = c0
2493 : strairyn = c0
2494 : #endif
2495 : endif
2496 :
2497 : !-----------------------------------------------------------------
2498 : ! Update ice age
2499 : ! This is further adjusted for freezing in the thermodynamics.
2500 : ! Melting does not alter the ice age.
2501 : !-----------------------------------------------------------------
2502 :
2503 4017728 : if (tr_iage) call increment_age (dt, iage(n))
2504 4017728 : if (icepack_warnings_aborted(subname)) return
2505 4017728 : if (tr_FY) call update_FYarea (dt, &
2506 : lmask_n, lmask_s, &
2507 229760 : yday, FY(n))
2508 4017728 : if (icepack_warnings_aborted(subname)) return
2509 :
2510 : !-----------------------------------------------------------------
2511 : ! Vertical thermodynamics: Heat conduction, growth and melting.
2512 : !-----------------------------------------------------------------
2513 :
2514 4017728 : if (.not.(calc_Tsfc)) then
2515 :
2516 : ! If not calculating surface temperature and fluxes, set
2517 : ! surface fluxes (flatn, fsurfn, and fcondtopn) to be used
2518 : ! in thickness_changes
2519 :
2520 : ! hadgem routine sets fluxes to default values in ice-only mode
2521 0 : call set_sfcflux(aicen (n), &
2522 0 : flatn_f (n), fsensn_f (n), &
2523 0 : fcondtopn_f(n), &
2524 0 : fsurfn_f (n), &
2525 0 : flatn (n), fsensn (n), &
2526 0 : fsurfn (n), &
2527 0 : fcondtopn (n))
2528 0 : if (icepack_warnings_aborted(subname)) return
2529 : endif
2530 :
2531 : call thermo_vertical(nilyr, nslyr, &
2532 1361620 : dt, aicen (n), &
2533 1361620 : vicen (n), vsnon (n), &
2534 1361620 : Tsfc (n), zSin (:,n), &
2535 0 : zqin (:,n), zqsn (:,n), &
2536 1361620 : apnd (n), hpnd (n), &
2537 : tr_pond_topo, &
2538 : flw, potT, &
2539 : Qa, rhoa, &
2540 : fsnow, fpond, &
2541 : fbot, Tbot, &
2542 : Tsnice, sss, &
2543 : lhcoef, shcoef, &
2544 1361620 : fswsfcn (n), fswintn (n), &
2545 0 : Sswabsn(:,n), Iswabsn(:,n), &
2546 1361620 : fsurfn (n), fcondtopn(n), &
2547 1361620 : fcondbotn(n), &
2548 1361620 : fsensn (n), flatn (n), &
2549 : flwoutn, evapn, &
2550 : evapsn, evapin, &
2551 : freshn, fsaltn, &
2552 : fhocnn, &
2553 1361620 : melttn (n), meltsn (n), &
2554 1361620 : meltbn (n), &
2555 1361620 : congeln (n), snoicen (n), &
2556 : mlt_onset, frz_onset, &
2557 1361620 : yday, dsnown (n), &
2558 6740968 : prescribed_ice)
2559 :
2560 4017728 : if (icepack_warnings_aborted(subname)) then
2561 0 : call icepack_warnings_add(subname//' ice: Vertical thermo error: ')
2562 0 : return
2563 : endif
2564 :
2565 : !-----------------------------------------------------------------
2566 : ! Total absorbed shortwave radiation
2567 : !-----------------------------------------------------------------
2568 :
2569 4017728 : fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)
2570 :
2571 : !-----------------------------------------------------------------
2572 : ! Aerosol update
2573 : !-----------------------------------------------------------------
2574 :
2575 4017728 : if (tr_aero) then
2576 : call update_aerosol (dt, &
2577 : nilyr, nslyr, n_aero, &
2578 0 : melttn (n), meltsn (n), &
2579 0 : meltbn (n), congeln (n), &
2580 0 : snoicen (n), fsnow, &
2581 0 : aerosno(:,:,n), aeroice(:,:,n), &
2582 0 : aicen_init (n), vicen_init (n), &
2583 0 : vsnon_init (n), &
2584 0 : vicen (n), vsnon (n), &
2585 0 : aicen (n), &
2586 229765 : faero_atm , faero_ocn)
2587 229765 : if (icepack_warnings_aborted(subname)) return
2588 : endif
2589 :
2590 4017728 : if (tr_iso) then
2591 : call update_isotope (dt = dt, &
2592 : nilyr = nilyr, nslyr = nslyr, &
2593 0 : meltt = melttn(n),melts = meltsn(n), &
2594 0 : meltb = meltbn(n),congel=congeln(n), &
2595 0 : snoice=snoicen(n),evap=evapn, &
2596 0 : fsnow=fsnow, Tsfc=Tsfc(n), &
2597 0 : Qref_iso=Qrefn_iso(:), &
2598 0 : isosno=l_isosno(:,n),isoice=l_isoice(:,n), &
2599 0 : aice_old=aicen_init(n),vice_old=vicen_init(n), &
2600 0 : vsno_old=vsnon_init(n), &
2601 0 : vicen=vicen(n),vsnon=vsnon(n), &
2602 0 : aicen=aicen(n), &
2603 : fiso_atm=l_fiso_atm(:), &
2604 0 : fiso_evapn=fiso_evapn(:), &
2605 0 : fiso_ocnn=fiso_ocnn(:), &
2606 : HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, &
2607 229760 : H2_18O_ocn=l_H2_18O_ocn)
2608 229760 : if (icepack_warnings_aborted(subname)) return
2609 : endif
2610 : endif ! aicen_init
2611 :
2612 : !-----------------------------------------------------------------
2613 : ! Melt ponds
2614 : ! If using tr_pond_cesm, the full calculation is performed here.
2615 : ! If using tr_pond_topo, the rest of the calculation is done after
2616 : ! the surface fluxes are merged, below.
2617 : !-----------------------------------------------------------------
2618 :
2619 : !call ice_timer_start(timer_ponds)
2620 5982048 : if (tr_pond) then
2621 :
2622 4920240 : if (tr_pond_cesm) then
2623 307440 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2624 : call compute_ponds_cesm(dt, hi_min, &
2625 : pndaspect, rfrac, &
2626 0 : melttn(n), meltsn(n), &
2627 : frain, &
2628 0 : aicen (n), vicen (n), &
2629 0 : Tsfc (n), &
2630 307440 : apnd (n), hpnd (n))
2631 307440 : if (icepack_warnings_aborted(subname)) return
2632 :
2633 4612800 : elseif (tr_pond_lvl) then
2634 3822720 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2635 : call compute_ponds_lvl(dt, nilyr, &
2636 : ktherm, &
2637 : hi_min, &
2638 : dpscale, frzpnd, &
2639 : pndaspect, rfrac, &
2640 1495440 : melttn(n), meltsn(n), &
2641 : frain, Tair, &
2642 1495440 : fsurfn(n), &
2643 1495440 : dhsn (n), ffracn(n), &
2644 1495440 : aicen (n), vicen (n), &
2645 1495440 : vsnon (n), &
2646 0 : zqin(:,n), zSin(:,n), &
2647 1495440 : Tsfc (n), alvl (n), &
2648 1495440 : apnd (n), hpnd (n), &
2649 5318160 : ipnd (n))
2650 3822720 : if (icepack_warnings_aborted(subname)) return
2651 :
2652 790080 : elseif (tr_pond_topo) then
2653 790080 : if (aicen_init(n) > puny) then
2654 :
2655 : ! collect liquid water in ponds
2656 : ! assume salt still runs off
2657 457021 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2658 80554 : pond = rfrac/rhofresh * (melttn(n)*rhoi &
2659 80554 : + meltsn(n)*rhos &
2660 457021 : + frain *dt)
2661 :
2662 : ! if pond does not exist, create new pond over full ice area
2663 : ! otherwise increase pond depth without changing pond area
2664 457021 : if (apnd(n) < puny) then
2665 272114 : hpnd(n) = c0
2666 272114 : apnd(n) = c1
2667 : endif
2668 457021 : hpnd(n) = (pond + hpnd(n)*apnd(n)) / apnd(n)
2669 457021 : fpond = fpond + pond * aicen(n) ! m
2670 : endif ! aicen_init
2671 : endif
2672 :
2673 : endif ! tr_pond
2674 : !call ice_timer_stop(timer_ponds)
2675 :
2676 : !-----------------------------------------------------------------
2677 : ! Increment area-weighted fluxes.
2678 : !-----------------------------------------------------------------
2679 :
2680 5982048 : if (aicen_init(n) > puny) &
2681 1361620 : call merge_fluxes (aicen_init(n), &
2682 : flw, &
2683 : strairxn, strairyn, &
2684 : Cdn_atm_ratio_n, &
2685 1361620 : fsurfn(n), fcondtopn(n), &
2686 1361620 : fcondbotn(n), &
2687 1361620 : fsensn(n), flatn(n), &
2688 : fswabsn, flwoutn, &
2689 : evapn, &
2690 : evapsn, evapin, &
2691 : Trefn, Qrefn, &
2692 : freshn, fsaltn, &
2693 1361620 : fhocnn, fswthrun(n), &
2694 : strairxT, strairyT, &
2695 : Cdn_atm_ratio, &
2696 : fsurf, fcondtop, &
2697 : fcondbot, &
2698 : fsens, flat, &
2699 : fswabs, flwout, &
2700 : evap, &
2701 : evaps, evapi, &
2702 : Tref, Qref, &
2703 : fresh, fsalt, &
2704 : fhocn, fswthru, &
2705 1361620 : melttn (n), meltsn(n), &
2706 1361620 : meltbn (n), congeln(n), &
2707 1361620 : snoicen(n), &
2708 : meltt, melts, &
2709 : meltb, congel, &
2710 : snoice, &
2711 : Uref=Uref, Urefn=Urefn, &
2712 : Qref_iso=l_Qref_iso, &
2713 0 : Qrefn_iso=Qrefn_iso, &
2714 : fiso_ocn=l_fiso_ocn, &
2715 0 : fiso_ocnn=fiso_ocnn, &
2716 : fiso_evap=l_fiso_evap, &
2717 4017728 : fiso_evapn=fiso_evapn)
2718 :
2719 7255680 : if (icepack_warnings_aborted(subname)) return
2720 :
2721 : enddo ! ncat
2722 :
2723 8178000 : if (present(isosno) ) isosno = l_isosno
2724 8178000 : if (present(isoice) ) isoice = l_isoice
2725 5094528 : if (present(Qa_iso) ) Qa_iso = l_Qa_iso
2726 5094528 : if (present(Qref_iso) ) Qref_iso = l_Qref_iso
2727 5094528 : if (present(fiso_atm) ) fiso_atm = l_fiso_atm
2728 5094528 : if (present(fiso_ocn) ) fiso_ocn = l_fiso_ocn
2729 5094528 : if (present(fiso_evap)) fiso_evap= l_fiso_evap
2730 1273632 : deallocate(l_isosno)
2731 1273632 : deallocate(l_isoice)
2732 1273632 : deallocate(l_Qa_iso)
2733 1273632 : deallocate(l_Qref_iso)
2734 1273632 : deallocate(l_fiso_atm)
2735 1273632 : deallocate(l_fiso_ocn)
2736 1273632 : deallocate(l_fiso_evap)
2737 :
2738 : !-----------------------------------------------------------------
2739 : ! Calculate ponds from the topographic scheme
2740 : !-----------------------------------------------------------------
2741 : !call ice_timer_start(timer_ponds)
2742 1273632 : if (tr_pond_topo) then
2743 : call compute_ponds_topo(dt, ncat, nilyr, &
2744 : ktherm, heat_capacity, &
2745 0 : aice, aicen, &
2746 0 : vice, vicen, &
2747 0 : vsno, vsnon, &
2748 : meltt, &
2749 : fsurf, fpond, &
2750 0 : Tsfc, Tf, &
2751 0 : zqin, zSin, &
2752 158016 : apnd, hpnd, ipnd )
2753 158016 : if (icepack_warnings_aborted(subname)) return
2754 : endif
2755 : !call ice_timer_stop(timer_ponds)
2756 :
2757 1273632 : end subroutine icepack_step_therm1
2758 :
2759 : !=======================================================================
2760 :
2761 : end module icepack_therm_vertical
2762 :
2763 : !=======================================================================
|
