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 4420832 : subroutine thermo_vertical (nilyr, nslyr, &
79 : dt, aicen, &
80 : vicen, vsnon, &
81 4420832 : Tsf, zSin, &
82 4420832 : 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 4420832 : 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 1617332 : dhi , & ! change in ice thickness
208 1617332 : dhs ! change in snow thickness
209 :
210 : ! 2D state variables (thickness, temperature)
211 :
212 : real (kind=dbl_kind) :: &
213 1617332 : hilyr , & ! ice layer thickness
214 1617332 : hslyr , & ! snow layer thickness
215 1617332 : hin , & ! ice thickness (m)
216 1617332 : hsn , & ! snow thickness (m)
217 1617332 : hsn_new , & ! thickness of new snow (m)
218 1617332 : worki , & ! local work array
219 1617332 : works ! local work array
220 :
221 : real (kind=dbl_kind), dimension (nilyr) :: &
222 21136604 : zTin ! internal ice layer temperatures
223 :
224 : real (kind=dbl_kind), dimension (nslyr) :: &
225 7977712 : zTsn ! internal snow layer temperatures
226 :
227 : ! other 2D flux and energy variables
228 :
229 : real (kind=dbl_kind) :: &
230 1617332 : einit , & ! initial energy of melting (J m-2)
231 1617332 : efinal , & ! final energy of melting (J m-2)
232 1617332 : einter ! intermediate energy
233 :
234 : real (kind=dbl_kind) :: &
235 1617332 : fadvocn ! advective heat flux to ocean
236 :
237 : character(len=*),parameter :: subname='(thermo_vertical)'
238 :
239 : !-----------------------------------------------------------------
240 : ! Initialize
241 : !-----------------------------------------------------------------
242 :
243 4420832 : flwoutn = c0
244 4420832 : evapn = c0
245 4420832 : evapsn = c0
246 4420832 : evapin = c0
247 4420832 : freshn = c0
248 4420832 : fsaltn = c0
249 4420832 : fhocnn = c0
250 4420832 : fadvocn = c0
251 :
252 4420832 : meltt = c0
253 4420832 : meltb = c0
254 4420832 : melts = c0
255 4420832 : congel = c0
256 4420832 : snoice = c0
257 4420832 : dsnow = c0
258 9750576 : zTsn(:) = c0
259 33547018 : zTin(:) = c0
260 :
261 4420832 : if (calc_Tsfc) then
262 4420832 : fsensn = c0
263 4420832 : flatn = c0
264 4420832 : fsurfn = c0
265 4420832 : 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 4420832 : einit )
281 4420832 : if (icepack_warnings_aborted(subname)) return
282 :
283 : ! Save initial ice and snow thickness (for fresh and fsalt)
284 4420832 : worki = hin
285 4420832 : works = hsn
286 :
287 : !-----------------------------------------------------------------
288 : ! Compute new surface temperature and internal ice and snow
289 : ! temperatures.
290 : !-----------------------------------------------------------------
291 :
292 4420832 : if (heat_capacity) then ! usual case
293 :
294 4117559 : 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 3232570 : fadvocn, snoice)
314 3232570 : 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 4420832 : einter = c0
370 9750576 : do k = 1, nslyr
371 9750576 : einter = einter + hslyr * zqsn(k)
372 : enddo ! k
373 33547018 : do k = 1, nilyr
374 33547018 : einter = einter + hilyr * zqin(k)
375 : enddo ! k
376 :
377 4420832 : Tsnice = c0
378 4420832 : if ((hslyr+hilyr) > puny) then
379 4420832 : if (hslyr > puny) then
380 3884832 : Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr)
381 : else
382 536000 : Tsnice = Tsf
383 : endif
384 : endif
385 :
386 4420832 : 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 4420832 : dsnow)
412 4420832 : 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 4420832 : fadvocn, fbot )
426 4420832 : 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 4420832 : dhi = hin - worki
448 4420832 : dhs = hsn - works - hsn_new
449 :
450 4420832 : freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt
451 4420832 : fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
452 4420832 : fhocnn = fhocnn + fadvocn ! for ktherm=2
453 :
454 4420832 : 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 4420832 : vicen, vsnon)
470 4420832 : 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 1384704 : subroutine frzmlt_bottom_lateral (dt, ncat, &
488 : nilyr, nslyr, &
489 : aice, frzmlt, &
490 2769408 : vicen, vsnon, &
491 2769408 : qicen, qsnon, &
492 : sst, Tf, &
493 : ustar_min, &
494 509328 : 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 509328 : etot , & ! total energy in column
543 509328 : qavg ! average enthalpy in column (approximate)
544 :
545 : real (kind=dbl_kind) :: &
546 509328 : deltaT , & ! SST - Tbot >= 0
547 509328 : ustar , & ! skin friction velocity for fbot (m/s)
548 509328 : wlat , & ! lateral melt rate (m/s)
549 509328 : xtmp ! temporary variable
550 :
551 : ! Parameters for bottom melting
552 :
553 : real (kind=dbl_kind) :: &
554 509328 : 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 1384704 : rside = c0
573 1384704 : fside = c0
574 1384704 : Tbot = Tf
575 1384704 : fbot = c0
576 1384704 : wlat = c0
577 :
578 1384704 : 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 510676 : deltaT = max((sst-Tbot),c0)
586 :
587 : ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
588 510676 : ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow)
589 510676 : ustar = max (ustar,ustar_min)
590 :
591 510676 : 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 510676 : cpchr = -cp_ocn*rhow*0.006_dbl_kind
598 : endif
599 :
600 510676 : fbot = cpchr * deltaT * ustar ! < 0
601 510676 : 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 510676 : wlat = m1 * deltaT**m2 ! Maykut & Perovich
613 510676 : rside = wlat*dt*pi/(floeshape*floediam) ! Steele
614 510676 : rside = max(c0,min(rside,c1))
615 :
616 : !-----------------------------------------------------------------
617 : ! Compute heat flux associated with this value of rside.
618 : !-----------------------------------------------------------------
619 :
620 3061928 : do n = 1, ncat
621 :
622 2551252 : etot = c0
623 2551252 : qavg = c0
624 :
625 : ! melting energy/unit area in each column, etot < 0
626 :
627 6067704 : do k = 1, nslyr
628 3516452 : etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
629 6067704 : qavg = qavg + qsnon(k,n)
630 : enddo
631 :
632 19123634 : do k = 1, nilyr
633 16572382 : etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind)
634 19123634 : qavg = qavg + qicen(k,n)
635 : enddo ! nilyr
636 :
637 : ! lateral heat flux, fside < 0
638 3061928 : if (tr_fsd) then ! floe size distribution
639 159235 : fside = fside + wlat*qavg
640 : else ! default floe size
641 2392017 : 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 510676 : xtmp = frzmlt/(fbot + fside + puny)
651 510676 : xtmp = min(xtmp, c1)
652 510676 : fbot = fbot * xtmp
653 510676 : rside = rside * xtmp
654 510676 : fside = fside * xtmp
655 :
656 : endif
657 :
658 1384704 : 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 4420832 : subroutine init_vertical_profile(nilyr, nslyr, &
670 : aicen, vicen, &
671 : vsnon, &
672 : hin, hilyr, &
673 : hsn, hslyr, &
674 4420832 : zqin, zTin, &
675 4420832 : zqsn, zTsn, &
676 4420832 : 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 16715772 : Tmlts ! melting temperature
712 :
713 : integer (kind=int_kind) :: &
714 : k ! ice layer index
715 :
716 : real (kind=dbl_kind) :: &
717 1617332 : rnslyr, & ! real(nslyr)
718 1617332 : 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 4420832 : rnslyr = real(nslyr,kind=dbl_kind)
733 :
734 4420832 : tsno_high = .false.
735 4420832 : tice_high = .false.
736 4420832 : tsno_low = .false.
737 4420832 : tice_low = .false.
738 :
739 4420832 : einit = c0
740 :
741 : !-----------------------------------------------------------------
742 : ! Surface temperature, ice and snow thickness
743 : ! Initialize internal energy
744 : !-----------------------------------------------------------------
745 :
746 4420832 : hin = vicen / aicen
747 4420832 : hsn = vsnon / aicen
748 4420832 : hilyr = hin / real(nilyr,kind=dbl_kind)
749 4420832 : 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 9750576 : 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 5329744 : if (hslyr > hs_min/rnslyr .and. heat_capacity) then
767 : ! zqsn < 0
768 1156213 : Tmax = -zqsn(k)*puny*rnslyr / &
769 3205826 : (rhos*cp_ice*vsnon)
770 : else
771 2123918 : zqsn (k) = -rhos * Lfresh
772 2123918 : Tmax = puny
773 : endif
774 :
775 : !-----------------------------------------------------------------
776 : ! Compute snow temperatures from enthalpies.
777 : ! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0.
778 : !-----------------------------------------------------------------
779 5329744 : zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice
780 :
781 : !-----------------------------------------------------------------
782 : ! Check for zTsn > Tmax (allowing for roundoff error) and zTsn < Tmin.
783 : !-----------------------------------------------------------------
784 9750576 : if (zTsn(k) > Tmax) then
785 0 : tsno_high = .true.
786 5329744 : 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 4420832 : 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 4420832 : 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 9750576 : do k = 1, nslyr
852 :
853 5329744 : if (zTsn(k) > c0) then ! correct roundoff error
854 9878 : zTsn(k) = c0
855 9878 : zqsn(k) = -rhos*Lfresh
856 : endif
857 :
858 : !-----------------------------------------------------------------
859 : ! initial energy per unit area of ice/snow, relative to 0 C
860 : !-----------------------------------------------------------------
861 9750576 : einit = einit + hslyr*zqsn(k)
862 :
863 : enddo ! nslyr
864 :
865 33547018 : do k = 1, nilyr
866 :
867 : !---------------------------------------------------------------------
868 : ! Use initial salinity profile for thin ice
869 : !---------------------------------------------------------------------
870 :
871 29126186 : 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 29126186 : if (ktherm == 2) then
886 22627990 : 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 29126186 : if (ktherm == 2) then
902 22627990 : 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 29126186 : if (l_brine) then
908 28822913 : 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 29126186 : if (zTin(k) > Tmax) then
917 0 : tice_high = .true.
918 29126186 : 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 29126186 : 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 29126186 : 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 29126186 : 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 33547018 : 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 4420832 : subroutine thickness_changes (nilyr, nslyr, &
1028 : dt, yday, &
1029 : efinal, &
1030 : hin, hilyr, &
1031 : hsn, hslyr, &
1032 4420832 : 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 4420832 : 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 1617332 : esub , & ! energy for sublimation, > 0 (J m-2)
1115 1617332 : econ , & ! energy for condensation, < 0 (J m-2)
1116 1617332 : etop_mlt , & ! energy for top melting, > 0 (J m-2)
1117 1617332 : ebot_mlt , & ! energy for bottom melting, > 0 (J m-2)
1118 1617332 : ebot_gro , & ! energy for bottom growth, < 0 (J m-2)
1119 1617332 : emlt_atm , & ! total energy of brine, from atmosphere (J m-2)
1120 1617332 : emlt_ocn ! total energy of brine, to ocean (J m-2)
1121 :
1122 : real (kind=dbl_kind) :: &
1123 1617332 : qbotmax , & ! max enthalpy of ice growing at bottom
1124 1617332 : dhi , & ! change in ice thickness
1125 1617332 : dhs , & ! change in snow thickness
1126 1617332 : Ti , & ! ice temperature
1127 1617332 : Ts , & ! snow temperature
1128 1617332 : qbot , & ! enthalpy of ice growing at bottom surface (J m-3)
1129 1617332 : qsub , & ! energy/unit volume to sublimate ice/snow (J m-3)
1130 1617332 : hqtot , & ! sum of h*q for two layers
1131 1617332 : wk1 , & ! temporary variable
1132 1617332 : zqsnew , & ! enthalpy of new snow (J m-3)
1133 1617332 : hstot , & ! snow thickness including new snow (m)
1134 1617332 : Tmlts ! melting temperature
1135 :
1136 : real (kind=dbl_kind), dimension (nilyr+1) :: &
1137 19519272 : zi1 , & ! depth of ice layer boundaries (m)
1138 19519272 : zi2 ! adjusted depths, with equal hilyr (m)
1139 :
1140 : real (kind=dbl_kind), dimension (nslyr+1) :: &
1141 10781212 : zs1 , & ! depth of snow layer boundaries (m)
1142 10781212 : zs2 ! adjusted depths, with equal hslyr (m)
1143 :
1144 : real (kind=dbl_kind), dimension (nilyr) :: &
1145 16715772 : dzi ! ice layer thickness after growth/melting
1146 :
1147 : real (kind=dbl_kind), dimension (nslyr) :: &
1148 9163880 : dzs ! snow layer thickness after growth/melting
1149 :
1150 : real (kind=dbl_kind), dimension (nilyr) :: &
1151 17901940 : qm , & ! energy of melting (J m-3) = zqin in BL99 formulation
1152 18333104 : qmlt ! enthalpy of melted ice (J m-3) = zero in BL99 formulation
1153 :
1154 : real (kind=dbl_kind) :: &
1155 1617332 : qbotm , &
1156 1617332 : qbotp , &
1157 1617332 : qbot0
1158 :
1159 : character(len=*),parameter :: subname='(thickness_changes)'
1160 :
1161 : !-----------------------------------------------------------------
1162 : ! Initialize
1163 : !-----------------------------------------------------------------
1164 :
1165 4420832 : dhi = c0
1166 4420832 : dhs = c0
1167 4420832 : hsn_new = c0
1168 :
1169 33547018 : do k = 1, nilyr
1170 33547018 : dzi(k) = hilyr
1171 : enddo
1172 :
1173 9750576 : do k = 1, nslyr
1174 9750576 : dzs(k) = hslyr
1175 : enddo
1176 :
1177 33547018 : do k = 1, nilyr
1178 29126186 : if (ktherm == 2) then
1179 22627990 : qmlt(k) = enthalpy_of_melting(zSin(k))
1180 : else
1181 6498196 : qmlt(k) = c0
1182 : endif
1183 29126186 : qm(k) = zqin(k) - qmlt(k)
1184 29126186 : emlt_atm = c0
1185 33547018 : 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 4420832 : 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 4420832 : wk1 = -flatn * dt
1222 4420832 : esub = max(wk1, c0) ! energy for sublimation, > 0
1223 4420832 : econ = min(wk1, c0) ! energy for condensation, < 0
1224 :
1225 4420832 : wk1 = (fsurfn - fcondtopn) * dt
1226 4420832 : etop_mlt = max(wk1, c0) ! etop_mlt > 0
1227 :
1228 4420832 : wk1 = (fcondbotn - fbot) * dt
1229 4420832 : ebot_mlt = max(wk1, c0) ! ebot_mlt > 0
1230 4420832 : 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 4420832 : evapn = c0 ! initialize
1240 4420832 : evapsn = c0 ! initialize
1241 4420832 : evapin = c0 ! initialize
1242 :
1243 4420832 : if (hsn > puny) then ! add snow with enthalpy zqsn(1)
1244 3884832 : dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0
1245 3884832 : dzs(1) = dzs(1) + dhs
1246 3884832 : evapn = evapn + dhs*rhos
1247 3884832 : evapsn = evapsn + dhs*rhos
1248 : else ! add ice with enthalpy zqin(1)
1249 536000 : dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0
1250 536000 : dzi(1) = dzi(1) + dhi
1251 536000 : evapn = evapn + dhi*rhoi
1252 536000 : evapin = evapin + dhi*rhoi
1253 : ! enthalpy of melt water
1254 536000 : emlt_atm = emlt_atm - qmlt(1) * dhi
1255 : endif
1256 :
1257 : !--------------------------------------------------------------
1258 : ! Grow ice (bottom)
1259 : !--------------------------------------------------------------
1260 :
1261 4420832 : if (ktherm == 2) then
1262 :
1263 3232570 : qbotm = enthalpy_mush(Tbot, sss)
1264 3232570 : qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
1265 3232570 : qbot0 = qbotm - qbotp
1266 :
1267 3232570 : dhi = ebot_gro / qbotp ! dhi > 0
1268 :
1269 3232570 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
1270 3232570 : hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
1271 3232570 : 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 4420832 : dzi(nilyr) = dzi(nilyr) + dhi
1299 4420832 : if (dzi(nilyr) > puny) then
1300 4420832 : zqin(nilyr) = hqtot / dzi(nilyr)
1301 4420832 : if (ktherm == 2) then
1302 3232570 : zSin(nilyr) = hstot / dzi(nilyr)
1303 3232570 : qmlt(nilyr) = enthalpy_of_melting(zSin(nilyr))
1304 : else
1305 1188262 : qmlt(nilyr) = c0
1306 : endif
1307 4420832 : 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 4420832 : congel = congel + dhi
1316 4420832 : if (dhi > puny .and. frz_onset < puny) &
1317 93 : frz_onset = yday
1318 :
1319 9750576 : do k = 1, nslyr
1320 :
1321 : !--------------------------------------------------------------
1322 : ! Remove internal snow melt
1323 : !--------------------------------------------------------------
1324 :
1325 5329744 : if (ktherm == 2 .and. zqsn(k) > -rhos * Lfresh) then
1326 :
1327 26715 : dhs = max(-dzs(k), &
1328 60131 : -((zqsn(k) + rhos*Lfresh) / (rhos*Lfresh)) * dzs(k))
1329 60131 : dzs(k) = dzs(k) + dhs
1330 60131 : zqsn(k) = -rhos * Lfresh
1331 60131 : melts = melts - dhs
1332 : ! delta E = zqsn(k) + rhos * Lfresh
1333 :
1334 : endif
1335 :
1336 : !--------------------------------------------------------------
1337 : ! Sublimation of snow (evapn < 0)
1338 : !--------------------------------------------------------------
1339 :
1340 5329744 : qsub = zqsn(k) - rhos*Lvap ! qsub < 0
1341 5329744 : dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0
1342 5329744 : dzs(k) = dzs(k) + dhs
1343 5329744 : esub = esub - dhs*qsub
1344 5329744 : esub = max(esub, c0) ! in case of roundoff error
1345 5329744 : evapn = evapn + dhs*rhos
1346 5329744 : evapsn = evapsn + dhs*rhos
1347 :
1348 : !--------------------------------------------------------------
1349 : ! Melt snow (top)
1350 : !--------------------------------------------------------------
1351 :
1352 5329744 : dhs = max(-dzs(k), etop_mlt/zqsn(k))
1353 5329744 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1354 5329744 : etop_mlt = etop_mlt - dhs*zqsn(k)
1355 5329744 : etop_mlt = max(etop_mlt, c0) ! in case of roundoff error
1356 :
1357 : ! history diagnostics
1358 5329744 : if (dhs < -puny .and. mlt_onset < puny) &
1359 22 : mlt_onset = yday
1360 9750576 : melts = melts - dhs
1361 :
1362 : enddo ! nslyr
1363 :
1364 33547018 : do k = 1, nilyr
1365 :
1366 : !--------------------------------------------------------------
1367 : ! Sublimation of ice (evapn < 0)
1368 : !--------------------------------------------------------------
1369 :
1370 29126186 : qsub = qm(k) - rhoi*Lvap ! qsub < 0
1371 29126186 : dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0
1372 29126186 : dzi(k) = dzi(k) + dhi
1373 29126186 : esub = esub - dhi*qsub
1374 29126186 : esub = max(esub, c0)
1375 29126186 : evapn = evapn + dhi*rhoi
1376 29126186 : evapin = evapin + dhi*rhoi
1377 29126186 : emlt_ocn = emlt_ocn - qmlt(k) * dhi
1378 :
1379 : !--------------------------------------------------------------
1380 : ! Melt ice (top)
1381 : !--------------------------------------------------------------
1382 :
1383 29126186 : if (qm(k) < c0) then
1384 29126038 : dhi = max(-dzi(k), etop_mlt/qm(k))
1385 : else
1386 148 : qm(k) = c0
1387 148 : dhi = -dzi(k)
1388 : endif
1389 29126186 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1390 :
1391 29126186 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1392 29126186 : etop_mlt = max(etop_mlt - dhi*qm(k), c0)
1393 :
1394 : ! history diagnostics
1395 29126186 : if (dhi < -puny .and. mlt_onset < puny) &
1396 92 : mlt_onset = yday
1397 33547018 : meltt = meltt - dhi
1398 :
1399 : enddo ! nilyr
1400 :
1401 33547018 : do k = nilyr, 1, -1
1402 :
1403 : !--------------------------------------------------------------
1404 : ! Melt ice (bottom)
1405 : !--------------------------------------------------------------
1406 :
1407 29126186 : if (qm(k) < c0) then
1408 29126038 : dhi = max(-dzi(k), ebot_mlt/qm(k))
1409 : else
1410 148 : qm(k) = c0
1411 148 : dhi = -dzi(k)
1412 : endif
1413 29126186 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1414 :
1415 29126186 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1416 29126186 : ebot_mlt = max(ebot_mlt - dhi*qm(k), c0)
1417 :
1418 : ! history diagnostics
1419 33547018 : meltb = meltb -dhi
1420 :
1421 : enddo ! nilyr
1422 :
1423 9750576 : do k = nslyr, 1, -1
1424 :
1425 : !--------------------------------------------------------------
1426 : ! Melt snow (only if all the ice has melted)
1427 : !--------------------------------------------------------------
1428 :
1429 5329744 : dhs = max(-dzs(k), ebot_mlt/zqsn(k))
1430 5329744 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1431 5329744 : ebot_mlt = ebot_mlt - dhs*zqsn(k)
1432 5329744 : ebot_mlt = max(ebot_mlt, c0)
1433 :
1434 : ! Add this to the snow melt (J. Zhu)
1435 9750576 : melts = melts - dhs
1436 :
1437 : enddo ! nslyr
1438 :
1439 : !-----------------------------------------------------------------
1440 : ! Compute heat flux used by the ice (<=0).
1441 : ! fhocn is the available ocean heat that is left after use by ice
1442 : !-----------------------------------------------------------------
1443 :
1444 : fhocnn = fbot &
1445 4420832 : + (esub + etop_mlt + ebot_mlt)/dt
1446 :
1447 : !---!-----------------------------------------------------------------
1448 : !---! Add new snowfall at top surface.
1449 : !---!-----------------------------------------------------------------
1450 :
1451 : !----------------------------------------------------------------
1452 : ! NOTE: If heat flux diagnostics are to work, new snow should
1453 : ! have T = 0 (i.e. q = -rhos*Lfresh) and should not be
1454 : ! converted to rain.
1455 : !----------------------------------------------------------------
1456 :
1457 4420832 : if (fsnow > c0) then
1458 :
1459 3524581 : hsn_new = fsnow/rhos * dt
1460 3524581 : zqsnew = -rhos*Lfresh
1461 3524581 : hstot = dzs(1) + hsn_new
1462 :
1463 3524581 : if (hstot > c0) then
1464 1266532 : zqsn(1) = (dzs(1) * zqsn(1) &
1465 3524581 : + hsn_new * zqsnew) / hstot
1466 : ! avoid roundoff errors
1467 3524581 : zqsn(1) = min(zqsn(1), -rhos*Lfresh)
1468 :
1469 3524581 : dzs(1) = hstot
1470 : endif
1471 : endif
1472 :
1473 : !-----------------------------------------------------------------
1474 : ! Find the new ice and snow thicknesses.
1475 : !-----------------------------------------------------------------
1476 :
1477 4420832 : hin = c0
1478 4420832 : hsn = c0
1479 :
1480 33547018 : do k = 1, nilyr
1481 33547018 : hin = hin + dzi(k)
1482 : enddo ! k
1483 :
1484 9750576 : do k = 1, nslyr
1485 5329744 : hsn = hsn + dzs(k)
1486 9750576 : dsnow = dsnow + dzs(k) - hslyr
1487 : enddo ! k
1488 :
1489 : !-------------------------------------------------------------------
1490 : ! Convert snow to ice if snow lies below freeboard.
1491 : !-------------------------------------------------------------------
1492 :
1493 4420832 : if (ktherm /= 2) &
1494 : call freeboard (nslyr, &
1495 : snoice, &
1496 : hin, hsn, &
1497 0 : zqin, zqsn, &
1498 0 : dzi, dzs, &
1499 1188262 : dsnow)
1500 4420832 : if (icepack_warnings_aborted(subname)) return
1501 :
1502 : !---!-------------------------------------------------------------------
1503 : !---! Repartition the ice and snow into equal-thickness layers,
1504 : !---! conserving energy.
1505 : !---!-------------------------------------------------------------------
1506 :
1507 : !-----------------------------------------------------------------
1508 : ! Compute desired layer thicknesses.
1509 : !-----------------------------------------------------------------
1510 :
1511 4420832 : if (hin > c0) then
1512 4420825 : hilyr = hin / real(nilyr,kind=dbl_kind)
1513 : else
1514 7 : hin = c0
1515 7 : hilyr = c0
1516 : endif
1517 4420832 : if (hsn > c0) then
1518 3862640 : hslyr = hsn / real(nslyr,kind=dbl_kind)
1519 : else
1520 558192 : hsn = c0
1521 558192 : hslyr = c0
1522 : endif
1523 :
1524 : !-----------------------------------------------------------------
1525 : ! Compute depths zi1 of old layers (unequal thickness).
1526 : ! Compute depths zi2 of new layers (equal thickness).
1527 : !-----------------------------------------------------------------
1528 :
1529 4420832 : zi1(1) = c0
1530 4420832 : zi1(1+nilyr) = hin
1531 :
1532 4420832 : zi2(1) = c0
1533 4420832 : zi2(1+nilyr) = hin
1534 :
1535 4420832 : if (heat_capacity) then
1536 :
1537 28822913 : do k = 1, nilyr-1
1538 24705354 : zi1(k+1) = zi1(k) + dzi(k)
1539 28822913 : zi2(k+1) = zi2(k) + hilyr
1540 : enddo
1541 :
1542 : !-----------------------------------------------------------------
1543 : ! Conserving energy, compute the enthalpy of the new equal layers.
1544 : !-----------------------------------------------------------------
1545 :
1546 : call adjust_enthalpy (nilyr, &
1547 0 : zi1, zi2, &
1548 : hilyr, hin, &
1549 4117559 : zqin)
1550 4117559 : if (icepack_warnings_aborted(subname)) return
1551 :
1552 4117559 : if (ktherm == 2) &
1553 : call adjust_enthalpy (nilyr, &
1554 0 : zi1, zi2, &
1555 : hilyr, hin, &
1556 3232570 : zSin)
1557 4117559 : if (icepack_warnings_aborted(subname)) return
1558 :
1559 : else ! zero layer (nilyr=1)
1560 :
1561 303273 : zqin(1) = -rhoi * Lfresh
1562 303273 : zqsn(1) = -rhos * Lfresh
1563 :
1564 : endif
1565 :
1566 4420832 : if (nslyr > 1) then
1567 :
1568 : !-----------------------------------------------------------------
1569 : ! Compute depths zs1 of old layers (unequal thickness).
1570 : ! Compute depths zs2 of new layers (equal thickness).
1571 : !-----------------------------------------------------------------
1572 :
1573 227228 : zs1(1) = c0
1574 227228 : zs1(1+nslyr) = hsn
1575 :
1576 227228 : zs2(1) = c0
1577 227228 : zs2(1+nslyr) = hsn
1578 :
1579 1136140 : do k = 1, nslyr-1
1580 908912 : zs1(k+1) = zs1(k) + dzs(k)
1581 1136140 : zs2(k+1) = zs2(k) + hslyr
1582 : enddo
1583 :
1584 : !-----------------------------------------------------------------
1585 : ! Conserving energy, compute the enthalpy of the new equal layers.
1586 : !-----------------------------------------------------------------
1587 :
1588 : call adjust_enthalpy (nslyr, &
1589 0 : zs1, zs2, &
1590 : hslyr, hsn, &
1591 227228 : zqsn)
1592 227228 : if (icepack_warnings_aborted(subname)) return
1593 :
1594 : endif ! nslyr > 1
1595 :
1596 : !-----------------------------------------------------------------
1597 : ! Remove very thin snow layers (ktherm = 2)
1598 : !-----------------------------------------------------------------
1599 :
1600 4420832 : if (ktherm == 2) then
1601 6465140 : do k = 1, nslyr
1602 6465140 : if (hsn <= puny) then
1603 : fhocnn = fhocnn &
1604 80955 : + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
1605 80955 : zqsn(k) = -rhos*Lfresh
1606 80955 : hslyr = c0
1607 : endif
1608 : enddo
1609 : endif
1610 :
1611 : !-----------------------------------------------------------------
1612 : ! Compute final ice-snow energy, including the energy of
1613 : ! sublimated/condensed ice.
1614 : !-----------------------------------------------------------------
1615 :
1616 4420832 : efinal = -evapn*Lvap
1617 4420832 : evapn = evapn/dt
1618 4420832 : evapsn = evapsn/dt
1619 4420832 : evapin = evapin/dt
1620 :
1621 9750576 : do k = 1, nslyr
1622 9750576 : efinal = efinal + hslyr*zqsn(k)
1623 : enddo
1624 :
1625 33547018 : do k = 1, nilyr
1626 33547018 : efinal = efinal + hilyr*zqin(k)
1627 : enddo ! k
1628 :
1629 4420832 : if (ktherm < 2) then
1630 1188262 : emlt_atm = c0
1631 1188262 : emlt_ocn = c0
1632 : endif
1633 :
1634 : ! melt water is no longer zero enthalpy with ktherm=2
1635 4420832 : fhocnn = fhocnn + emlt_ocn/dt
1636 4420832 : efinal = efinal + emlt_atm ! for conservation check
1637 :
1638 : end subroutine thickness_changes
1639 :
1640 : !=======================================================================
1641 : !
1642 : ! If there is enough snow to lower the ice/snow interface below
1643 : ! sea level, convert enough snow to ice to bring the interface back
1644 : ! to sea level.
1645 : !
1646 : ! authors William H. Lipscomb, LANL
1647 : ! Elizabeth C. Hunke, LANL
1648 :
1649 1188262 : subroutine freeboard (nslyr, &
1650 : snoice, &
1651 : hin, hsn, &
1652 2376524 : zqin, zqsn, &
1653 2376524 : dzi, dzs, &
1654 : dsnow)
1655 :
1656 : integer (kind=int_kind), intent(in) :: &
1657 : nslyr ! number of snow layers
1658 :
1659 : ! real (kind=dbl_kind), intent(in) :: &
1660 : ! dt ! time step
1661 :
1662 : real (kind=dbl_kind), &
1663 : intent(inout) :: &
1664 : snoice , & ! snow-ice formation (m/step-->cm/day)
1665 : dsnow ! change in snow thickness after snow-ice formation (m)
1666 : ! iage ! ice age (s)
1667 :
1668 : real (kind=dbl_kind), &
1669 : intent(inout) :: &
1670 : hin , & ! ice thickness (m)
1671 : hsn ! snow thickness (m)
1672 :
1673 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1674 : zqsn ! snow layer enthalpy (J m-3)
1675 :
1676 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1677 : zqin , & ! ice layer enthalpy (J m-3)
1678 : dzi , & ! ice layer thicknesses (m)
1679 : dzs ! snow layer thicknesses (m)
1680 :
1681 : ! local variables
1682 :
1683 : integer (kind=int_kind) :: &
1684 : k ! vertical index
1685 :
1686 : real (kind=dbl_kind) :: &
1687 422506 : dhin , & ! change in ice thickness (m)
1688 422506 : dhsn , & ! change in snow thickness (m)
1689 422506 : hqs ! sum of h*q for snow (J m-2)
1690 :
1691 : real (kind=dbl_kind) :: &
1692 422506 : wk1 , & ! temporary variable
1693 422506 : dhs ! snow to remove from layer (m)
1694 :
1695 : character(len=*),parameter :: subname='(freeboard)'
1696 :
1697 : !-----------------------------------------------------------------
1698 : ! Determine whether snow lies below freeboard.
1699 : !-----------------------------------------------------------------
1700 :
1701 1188262 : dhin = c0
1702 1188262 : dhsn = c0
1703 1188262 : hqs = c0
1704 :
1705 1188262 : wk1 = hsn - hin*(rhow-rhoi)/rhos
1706 :
1707 1188262 : if (wk1 > puny .and. hsn > puny) then ! snow below freeboard
1708 17248 : dhsn = min(wk1*rhoi/rhow, hsn) ! snow to remove
1709 17248 : dhin = dhsn * rhos/rhoi ! ice to add
1710 : endif
1711 :
1712 : !-----------------------------------------------------------------
1713 : ! Adjust snow layer thickness.
1714 : ! Compute energy to transfer from snow to ice.
1715 : !-----------------------------------------------------------------
1716 :
1717 3285436 : do k = nslyr, 1, -1
1718 3285436 : if (dhin > puny) then
1719 40816 : dhs = min(dhsn, dzs(k)) ! snow to remove from layer
1720 40816 : hsn = hsn - dhs
1721 40816 : dsnow = dsnow -dhs !new snow addition term
1722 40816 : dzs(k) = dzs(k) - dhs
1723 40816 : dhsn = dhsn - dhs
1724 40816 : dhsn = max(dhsn,c0)
1725 40816 : hqs = hqs + dhs * zqsn(k)
1726 : endif ! dhin > puny
1727 : enddo
1728 :
1729 : !-----------------------------------------------------------------
1730 : ! Transfer volume and energy from snow to top ice layer.
1731 : !-----------------------------------------------------------------
1732 :
1733 1188262 : if (dhin > puny) then
1734 : ! update ice age due to freezing (new ice age = dt)
1735 : ! if (tr_iage) &
1736 : ! iage = (iage*hin+dt*dhin)/(hin+dhin)
1737 :
1738 17248 : wk1 = dzi(1) + dhin
1739 17248 : hin = hin + dhin
1740 17248 : zqin(1) = (dzi(1)*zqin(1) + hqs) / wk1
1741 17248 : dzi(1) = wk1
1742 :
1743 : ! history diagnostic
1744 17248 : snoice = snoice + dhin
1745 : endif ! dhin > puny
1746 :
1747 1188262 : end subroutine freeboard
1748 :
1749 : !=======================================================================
1750 : !
1751 : ! Conserving energy, compute the new enthalpy of equal-thickness ice
1752 : ! or snow layers.
1753 : !
1754 : ! authors William H. Lipscomb, LANL
1755 : ! C. M. Bitz, UW
1756 :
1757 7577357 : subroutine adjust_enthalpy (nlyr, &
1758 7577357 : z1, z2, &
1759 : hlyr, hn, &
1760 7577357 : qn)
1761 :
1762 : integer (kind=int_kind), intent(in) :: &
1763 : nlyr ! number of layers (nilyr or nslyr)
1764 :
1765 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1766 : z1 , & ! interface depth for old, unequal layers (m)
1767 : z2 ! interface depth for new, equal layers (m)
1768 :
1769 : real (kind=dbl_kind), intent(in) :: &
1770 : hlyr ! new layer thickness (m)
1771 :
1772 : real (kind=dbl_kind), intent(in) :: &
1773 : hn ! total thickness (m)
1774 :
1775 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1776 : qn ! layer quantity (enthalpy, salinity...)
1777 :
1778 : ! local variables
1779 :
1780 : integer (kind=int_kind) :: &
1781 : k, k1, k2 ! vertical indices
1782 :
1783 : real (kind=dbl_kind) :: &
1784 2785426 : hovlp ! overlap between old and new layers (m)
1785 :
1786 : real (kind=dbl_kind) :: &
1787 2785426 : rhlyr ! 1./hlyr
1788 :
1789 : real (kind=dbl_kind), dimension (nlyr) :: &
1790 31706162 : hq ! h * q for a layer
1791 :
1792 : character(len=*),parameter :: subname='(adjust_enthalpy)'
1793 :
1794 : !-----------------------------------------------------------------
1795 : ! Compute reciprocal layer thickness.
1796 : !-----------------------------------------------------------------
1797 :
1798 7577357 : rhlyr = c0
1799 7577357 : if (hn > puny) rhlyr = c1 / hlyr
1800 :
1801 : !-----------------------------------------------------------------
1802 : ! Compute h*q for new layers (k2) given overlap with old layers (k1)
1803 : !-----------------------------------------------------------------
1804 :
1805 60164400 : do k2 = 1, nlyr
1806 60164400 : hq(k2) = c0
1807 : enddo ! k
1808 7577357 : k1 = 1
1809 7577357 : k2 = 1
1810 104392898 : do while (k1 <= nlyr .and. k2 <= nlyr)
1811 71259972 : hovlp = min (z1(k1+1), z2(k2+1)) &
1812 168075513 : - max (z1(k1), z2(k2))
1813 96815541 : hovlp = max (hovlp, c0)
1814 96815541 : hq(k2) = hq(k2) + hovlp*qn(k1)
1815 104392898 : if (z1(k1+1) > z2(k2+1)) then
1816 44228498 : k2 = k2 + 1
1817 : else
1818 52587043 : k1 = k1 + 1
1819 : endif
1820 : enddo ! while
1821 :
1822 : !-----------------------------------------------------------------
1823 : ! Compute new enthalpies.
1824 : !-----------------------------------------------------------------
1825 :
1826 60164400 : do k = 1, nlyr
1827 60164400 : qn(k) = hq(k) * rhlyr
1828 : enddo ! k
1829 :
1830 7577357 : end subroutine adjust_enthalpy
1831 :
1832 : !=======================================================================
1833 : !
1834 : ! Check for energy conservation by comparing the change in energy
1835 : ! to the net energy input.
1836 : !
1837 : ! authors William H. Lipscomb, LANL
1838 : ! C. M. Bitz, UW
1839 : ! Adrian K. Turner, LANL
1840 :
1841 4420832 : subroutine conservation_check_vthermo(dt, &
1842 : fsurfn, flatn, &
1843 : fhocnn, fswint, &
1844 : fsnow, &
1845 : einit, einter, &
1846 : efinal, &
1847 : fcondtopn,fcondbotn, &
1848 : fadvocn, fbot )
1849 :
1850 : real (kind=dbl_kind), intent(in) :: &
1851 : dt ! time step
1852 :
1853 : real (kind=dbl_kind), intent(in) :: &
1854 : fsurfn , & ! net flux to top surface, excluding fcondtopn
1855 : flatn , & ! surface downward latent heat (W m-2)
1856 : fhocnn , & ! fbot, corrected for any surplus energy
1857 : fswint , & ! SW absorbed in ice interior, below surface (W m-2)
1858 : fsnow , & ! snowfall rate (kg m-2 s-1)
1859 : fcondtopn , &
1860 : fadvocn , &
1861 : fbot
1862 :
1863 : real (kind=dbl_kind), intent(in) :: &
1864 : einit , & ! initial energy of melting (J m-2)
1865 : einter , & ! intermediate energy of melting (J m-2)
1866 : efinal , & ! final energy of melting (J m-2)
1867 : fcondbotn
1868 :
1869 : ! local variables
1870 :
1871 : real (kind=dbl_kind) :: &
1872 1617332 : einp , & ! energy input during timestep (J m-2)
1873 1617332 : ferr ! energy conservation error (W m-2)
1874 :
1875 : character(len=*),parameter :: subname='(conservation_check_vthermo)'
1876 :
1877 : !----------------------------------------------------------------
1878 : ! If energy is not conserved, print diagnostics and exit.
1879 : !----------------------------------------------------------------
1880 :
1881 : !-----------------------------------------------------------------
1882 : ! Note that fsurf - flat = fsw + flw + fsens; i.e., the latent
1883 : ! heat is not included in the energy input, since (efinal - einit)
1884 : ! is the energy change in the system ice + vapor, and the latent
1885 : ! heat lost by the ice is equal to that gained by the vapor.
1886 : !-----------------------------------------------------------------
1887 :
1888 : einp = (fsurfn - flatn + fswint - fhocnn &
1889 4420832 : - fsnow*Lfresh - fadvocn) * dt
1890 4420832 : ferr = abs(efinal-einit-einp) / dt
1891 :
1892 4420832 : if (ferr > ferrmax) then
1893 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1894 0 : call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" )
1895 :
1896 0 : write(warnstr,*) subname, 'Thermo energy conservation error'
1897 0 : call icepack_warnings_add(warnstr)
1898 0 : write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr
1899 0 : call icepack_warnings_add(warnstr)
1900 0 : write(warnstr,*) subname, 'Energy error (J) =', ferr*dt
1901 0 : call icepack_warnings_add(warnstr)
1902 0 : write(warnstr,*) subname, 'Initial energy =', einit
1903 0 : call icepack_warnings_add(warnstr)
1904 0 : write(warnstr,*) subname, 'Final energy =', efinal
1905 0 : call icepack_warnings_add(warnstr)
1906 0 : write(warnstr,*) subname, 'efinal - einit =', efinal-einit
1907 0 : call icepack_warnings_add(warnstr)
1908 0 : write(warnstr,*) subname, 'fsurfn,flatn,fswint,fhocn, fsnow*Lfresh:'
1909 0 : call icepack_warnings_add(warnstr)
1910 0 : write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn, fsnow*Lfresh
1911 0 : call icepack_warnings_add(warnstr)
1912 0 : write(warnstr,*) subname, 'Input energy =', einp
1913 0 : call icepack_warnings_add(warnstr)
1914 0 : write(warnstr,*) subname, 'fbot,fcondbot:'
1915 0 : call icepack_warnings_add(warnstr)
1916 0 : write(warnstr,*) subname, fbot,fcondbotn
1917 0 : call icepack_warnings_add(warnstr)
1918 :
1919 : ! if (ktherm == 2) then
1920 0 : write(warnstr,*) subname, 'Intermediate energy =', einter
1921 0 : call icepack_warnings_add(warnstr)
1922 0 : write(warnstr,*) subname, 'efinal - einter =', &
1923 0 : efinal-einter
1924 0 : call icepack_warnings_add(warnstr)
1925 0 : write(warnstr,*) subname, 'einter - einit =', &
1926 0 : einter-einit
1927 0 : call icepack_warnings_add(warnstr)
1928 0 : write(warnstr,*) subname, 'Conduction Error =', (einter-einit) &
1929 0 : - (fcondtopn*dt - fcondbotn*dt + fswint*dt)
1930 0 : call icepack_warnings_add(warnstr)
1931 0 : write(warnstr,*) subname, 'Melt/Growth Error =', (einter-einit) &
1932 0 : + ferr*dt - (fcondtopn*dt - fcondbotn*dt + fswint*dt)
1933 0 : call icepack_warnings_add(warnstr)
1934 0 : write(warnstr,*) subname, 'Advection Error =', fadvocn*dt
1935 0 : call icepack_warnings_add(warnstr)
1936 : ! endif
1937 :
1938 : ! write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn
1939 : ! call icepack_warnings_add(warnstr)
1940 :
1941 0 : write(warnstr,*) subname, 'dt*(fsurfn, flatn, fswint, fhocn, fsnow*Lfresh, fadvocn):'
1942 0 : call icepack_warnings_add(warnstr)
1943 0 : write(warnstr,*) subname, fsurfn*dt, flatn*dt, &
1944 0 : fswint*dt, fhocnn*dt, &
1945 0 : fsnow*Lfresh*dt, fadvocn*dt
1946 0 : call icepack_warnings_add(warnstr)
1947 0 : return
1948 : endif
1949 :
1950 : end subroutine conservation_check_vthermo
1951 :
1952 : !=======================================================================
1953 : !
1954 : ! Given the vertical thermo state variables (hin, hsn),
1955 : ! compute the new ice state variables (vicen, vsnon).
1956 : ! Zero out state variables if ice has melted entirely.
1957 : !
1958 : ! authors William H. Lipscomb, LANL
1959 : ! C. M. Bitz, UW
1960 : ! Elizabeth C. Hunke, LANL
1961 :
1962 4420832 : subroutine update_state_vthermo(nilyr, nslyr, &
1963 : Tf, Tsf, &
1964 : hin, hsn, &
1965 4420832 : zqin, zSin, &
1966 4420832 : zqsn, &
1967 : aicen, vicen, &
1968 : vsnon)
1969 :
1970 : integer (kind=int_kind), intent(in) :: &
1971 : nilyr , & ! number of ice layers
1972 : nslyr ! number of snow layers
1973 :
1974 : real (kind=dbl_kind), intent(in) :: &
1975 : Tf ! freezing temperature (C)
1976 :
1977 : real (kind=dbl_kind), intent(inout) :: &
1978 : Tsf ! ice/snow surface temperature, Tsfcn
1979 :
1980 : real (kind=dbl_kind), intent(in) :: &
1981 : hin , & ! ice thickness (m)
1982 : hsn ! snow thickness (m)
1983 :
1984 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1985 : zqin , & ! ice layer enthalpy (J m-3)
1986 : zSin , & ! ice salinity (ppt)
1987 : zqsn ! snow layer enthalpy (J m-3)
1988 :
1989 : real (kind=dbl_kind), intent(inout) :: &
1990 : aicen , & ! concentration of ice
1991 : vicen , & ! volume per unit area of ice (m)
1992 : vsnon ! volume per unit area of snow (m)
1993 :
1994 : ! local variables
1995 :
1996 : integer (kind=int_kind) :: &
1997 : k ! ice layer index
1998 :
1999 : character(len=*),parameter :: subname='(update_state_vthermo)'
2000 :
2001 4420832 : if (hin <= c0) then
2002 7 : aicen = c0
2003 7 : vicen = c0
2004 7 : vsnon = c0
2005 7 : Tsf = Tf
2006 56 : do k = 1, nilyr
2007 56 : zqin(k) = c0
2008 : enddo
2009 7 : if (ktherm == 2) then
2010 56 : do k = 1, nilyr
2011 56 : zSin(k) = c0
2012 : enddo
2013 : endif
2014 14 : do k = 1, nslyr
2015 14 : zqsn(k) = c0
2016 : enddo
2017 : else
2018 : ! aicen is already up to date
2019 4420825 : vicen = aicen * hin
2020 4420825 : vsnon = aicen * hsn
2021 : endif
2022 :
2023 4420832 : end subroutine update_state_vthermo
2024 :
2025 : !=======================================================================
2026 : !autodocument_start icepack_step_therm1
2027 : ! Driver for thermodynamic changes not needed for coupling:
2028 : ! transport in thickness space, lateral growth and melting.
2029 : !
2030 : ! authors: William H. Lipscomb, LANL
2031 : ! Elizabeth C. Hunke, LANL
2032 :
2033 1384704 : subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
2034 1384704 : aicen_init , &
2035 1384704 : vicen_init , vsnon_init , &
2036 1384704 : aice , aicen , &
2037 1384704 : vice , vicen , &
2038 1384704 : vsno , vsnon , &
2039 : uvel , vvel , &
2040 2769408 : Tsfc , zqsn , &
2041 2769408 : zqin , zSin , &
2042 1384704 : alvl , vlvl , &
2043 1384704 : apnd , hpnd , &
2044 1384704 : ipnd , &
2045 1384704 : iage , FY , &
2046 1384704 : aerosno , aeroice , &
2047 1384704 : isosno , isoice , &
2048 : uatm , vatm , &
2049 : wind , zlvl , &
2050 : Qa , rhoa , &
2051 1384704 : Qa_iso , &
2052 : Tair , Tref , &
2053 : Qref , Uref , &
2054 1384704 : Qref_iso , &
2055 : Cdn_atm_ratio, &
2056 : Cdn_ocn , Cdn_ocn_skin, &
2057 : Cdn_ocn_floe, Cdn_ocn_keel, &
2058 : Cdn_atm , Cdn_atm_skin, &
2059 : Cdn_atm_floe, Cdn_atm_pond, &
2060 : Cdn_atm_rdg , hfreebd , &
2061 : hdraft , hridge , &
2062 : distrdg , hkeel , &
2063 : dkeel , lfloe , &
2064 : dfloe , &
2065 : strax , stray , &
2066 : strairxT , strairyT , &
2067 : potT , sst , &
2068 : sss , Tf , &
2069 : strocnxT , strocnyT , &
2070 : fbot , &
2071 : Tbot , Tsnice , &
2072 : frzmlt , rside , &
2073 : fside , &
2074 : fsnow , frain , &
2075 : fpond , &
2076 1384704 : fsurf , fsurfn , &
2077 1384704 : fcondtop , fcondtopn , &
2078 1384704 : fcondbot , fcondbotn , &
2079 1384704 : fswsfcn , fswintn , &
2080 1384704 : fswthrun , &
2081 1384704 : fswthrun_vdr, &
2082 1384704 : fswthrun_vdf, &
2083 1384704 : fswthrun_idr, &
2084 1384704 : fswthrun_idf, &
2085 : fswabs , &
2086 : flwout , &
2087 1384704 : Sswabsn , Iswabsn , &
2088 : flw , &
2089 1384704 : fsens , fsensn , &
2090 1384704 : flat , flatn , &
2091 : evap , &
2092 : evaps , evapi , &
2093 : fresh , fsalt , &
2094 : fhocn , &
2095 : fswthru , &
2096 : fswthru_vdr , &
2097 : fswthru_vdf , &
2098 : fswthru_idr , &
2099 : fswthru_idf , &
2100 1384704 : flatn_f , fsensn_f , &
2101 1384704 : fsurfn_f , fcondtopn_f , &
2102 1384704 : faero_atm , faero_ocn , &
2103 1384704 : fiso_atm , fiso_ocn , &
2104 1384704 : fiso_evap , &
2105 : HDO_ocn , H2_16O_ocn , &
2106 : H2_18O_ocn , &
2107 1384704 : dhsn , ffracn , &
2108 1384704 : meltt , melttn , &
2109 1384704 : meltb , meltbn , &
2110 1384704 : melts , meltsn , &
2111 1384704 : congel , congeln , &
2112 1384704 : snoice , snoicen , &
2113 1384704 : dsnown , &
2114 : lmask_n , lmask_s , &
2115 : mlt_onset , frz_onset , &
2116 : yday , prescribed_ice)
2117 :
2118 : integer (kind=int_kind), intent(in) :: &
2119 : ncat , & ! number of thickness categories
2120 : nilyr , & ! number of ice layers
2121 : nslyr ! number of snow layers
2122 :
2123 : real (kind=dbl_kind), intent(in) :: &
2124 : dt , & ! time step
2125 : uvel , & ! x-component of velocity (m/s)
2126 : vvel , & ! y-component of velocity (m/s)
2127 : strax , & ! wind stress components (N/m^2)
2128 : stray , & !
2129 : yday ! day of year
2130 :
2131 : logical (kind=log_kind), intent(in) :: &
2132 : lmask_n , & ! northern hemisphere mask
2133 : lmask_s ! southern hemisphere mask
2134 :
2135 : logical (kind=log_kind), intent(in), optional :: &
2136 : prescribed_ice ! if .true., use prescribed ice instead of computed
2137 :
2138 : real (kind=dbl_kind), intent(inout) :: &
2139 : aice , & ! sea ice concentration
2140 : vice , & ! volume per unit area of ice (m)
2141 : vsno , & ! volume per unit area of snow (m)
2142 : zlvl , & ! atm level height (m)
2143 : uatm , & ! wind velocity components (m/s)
2144 : vatm , &
2145 : wind , & ! wind speed (m/s)
2146 : potT , & ! air potential temperature (K)
2147 : Tair , & ! air temperature (K)
2148 : Qa , & ! specific humidity (kg/kg)
2149 : rhoa , & ! air density (kg/m^3)
2150 : frain , & ! rainfall rate (kg/m^2 s)
2151 : fsnow , & ! snowfall rate (kg/m^2 s)
2152 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
2153 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
2154 : fsalt , & ! salt flux to ocean (kg/m^2/s)
2155 : fhocn , & ! net heat flux to ocean (W/m^2)
2156 : fswthru , & ! shortwave penetrating to ocean (W/m^2)
2157 : fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2)
2158 : fcondtop , & ! top surface conductive flux (W/m^2)
2159 : fcondbot , & ! bottom surface conductive flux (W/m^2)
2160 : fsens , & ! sensible heat flux (W/m^2)
2161 : flat , & ! latent heat flux (W/m^2)
2162 : fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2)
2163 : flw , & ! incoming longwave radiation (W/m^2)
2164 : flwout , & ! outgoing longwave radiation (W/m^2)
2165 : evap , & ! evaporative water flux (kg/m^2/s)
2166 : evaps , & ! evaporative water flux over snow (kg/m^2/s)
2167 : evapi , & ! evaporative water flux over ice (kg/m^2/s)
2168 : congel , & ! basal ice growth (m/step-->cm/day)
2169 : snoice , & ! snow-ice formation (m/step-->cm/day)
2170 : Tref , & ! 2m atm reference temperature (K)
2171 : Qref , & ! 2m atm reference spec humidity (kg/kg)
2172 : Uref , & ! 10m atm reference wind speed (m/s)
2173 : Cdn_atm , & ! atm drag coefficient
2174 : Cdn_ocn , & ! ocn drag coefficient
2175 : hfreebd , & ! freeboard (m)
2176 : hdraft , & ! draft of ice + snow column (Stoessel1993)
2177 : hridge , & ! ridge height
2178 : distrdg , & ! distance between ridges
2179 : hkeel , & ! keel depth
2180 : dkeel , & ! distance between keels
2181 : lfloe , & ! floe length
2182 : dfloe , & ! distance between floes
2183 : Cdn_atm_skin, & ! neutral skin drag coefficient
2184 : Cdn_atm_floe, & ! neutral floe edge drag coefficient
2185 : Cdn_atm_pond, & ! neutral pond edge drag coefficient
2186 : Cdn_atm_rdg , & ! neutral ridge drag coefficient
2187 : Cdn_ocn_skin, & ! skin drag coefficient
2188 : Cdn_ocn_floe, & ! floe edge drag coefficient
2189 : Cdn_ocn_keel, & ! keel drag coefficient
2190 : Cdn_atm_ratio,& ! ratio drag atm / neutral drag atm
2191 : strairxT , & ! stress on ice by air, x-direction
2192 : strairyT , & ! stress on ice by air, y-direction
2193 : strocnxT , & ! ice-ocean stress, x-direction
2194 : strocnyT , & ! ice-ocean stress, y-direction
2195 : fbot , & ! ice-ocean heat flux at bottom surface (W/m^2)
2196 : frzmlt , & ! freezing/melting potential (W/m^2)
2197 : rside , & ! fraction of ice that melts laterally
2198 : fside , & ! lateral heat flux (W/m^2)
2199 : sst , & ! sea surface temperature (C)
2200 : Tf , & ! freezing temperature (C)
2201 : Tbot , & ! ice bottom surface temperature (deg C)
2202 : Tsnice , & ! snow ice interface temperature (deg C)
2203 : sss , & ! sea surface salinity (ppt)
2204 : meltt , & ! top ice melt (m/step-->cm/day)
2205 : melts , & ! snow melt (m/step-->cm/day)
2206 : meltb , & ! basal ice melt (m/step-->cm/day)
2207 : mlt_onset , & ! day of year that sfc melting begins
2208 : frz_onset ! day of year that freezing begins (congel or frazil)
2209 :
2210 : real (kind=dbl_kind), intent(inout), optional :: &
2211 : fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
2212 : fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2)
2213 : fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2)
2214 : fswthru_idf ! nir dif shortwave penetrating to ocean (W/m^2)
2215 :
2216 : real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
2217 : Qa_iso , & ! isotope specific humidity (kg/kg)
2218 : Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg)
2219 : fiso_atm , & ! isotope deposition rate (kg/m^2 s)
2220 : fiso_ocn , & ! isotope flux to ocean (kg/m^2/s)
2221 : fiso_evap ! isotope evaporation (kg/m^2/s)
2222 :
2223 : real (kind=dbl_kind), optional, intent(in) :: &
2224 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
2225 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
2226 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
2227 :
2228 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
2229 : aicen_init , & ! fractional area of ice
2230 : vicen_init , & ! volume per unit area of ice (m)
2231 : vsnon_init , & ! volume per unit area of snow (m)
2232 : aicen , & ! concentration of ice
2233 : vicen , & ! volume per unit area of ice (m)
2234 : vsnon , & ! volume per unit area of snow (m)
2235 : Tsfc , & ! ice/snow surface temperature, Tsfcn
2236 : alvl , & ! level ice area fraction
2237 : vlvl , & ! level ice volume fraction
2238 : apnd , & ! melt pond area fraction
2239 : hpnd , & ! melt pond depth (m)
2240 : ipnd , & ! melt pond refrozen lid thickness (m)
2241 : iage , & ! volume-weighted ice age
2242 : FY , & ! area-weighted first-year ice area
2243 : fsurfn , & ! net flux to top surface, excluding fcondtop
2244 : fcondtopn , & ! downward cond flux at top surface (W m-2)
2245 : fcondbotn , & ! downward cond flux at bottom surface (W m-2)
2246 : flatn , & ! latent heat flux (W m-2)
2247 : fsensn , & ! sensible heat flux (W m-2)
2248 : fsurfn_f , & ! net flux to top surface, excluding fcondtop
2249 : fcondtopn_f , & ! downward cond flux at top surface (W m-2)
2250 : flatn_f , & ! latent heat flux (W m-2)
2251 : fsensn_f , & ! sensible heat flux (W m-2)
2252 : fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
2253 : fswthrun , & ! SW through ice to ocean (W/m^2)
2254 : fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
2255 : faero_atm , & ! aerosol deposition rate (kg/m^2 s)
2256 : faero_ocn , & ! aerosol flux to ocean (kg/m^2/s)
2257 : dhsn , & ! depth difference for snow on sea ice and pond ice
2258 : ffracn , & ! fraction of fsurfn used to melt ipond
2259 : meltsn , & ! snow melt (m)
2260 : melttn , & ! top ice melt (m)
2261 : meltbn , & ! bottom ice melt (m)
2262 : congeln , & ! congelation ice growth (m)
2263 : snoicen , & ! snow-ice growth (m)
2264 : dsnown ! change in snow thickness (m/step-->cm/day)
2265 :
2266 : real (kind=dbl_kind), optional, dimension(:), intent(inout) :: &
2267 : fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
2268 : fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
2269 : fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
2270 : fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
2271 :
2272 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
2273 : zqsn , & ! snow layer enthalpy (J m-3)
2274 : zqin , & ! ice layer enthalpy (J m-3)
2275 : zSin , & ! internal ice layer salinities
2276 : Sswabsn , & ! SW radiation absorbed in snow layers (W m-2)
2277 : Iswabsn ! SW radiation absorbed in ice layers (W m-2)
2278 :
2279 : real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: &
2280 : aerosno , & ! snow aerosol tracer (kg/m^2)
2281 : aeroice ! ice aerosol tracer (kg/m^2)
2282 :
2283 : real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: &
2284 : isosno , & ! snow isotope tracer (kg/m^2)
2285 : isoice ! ice isotope tracer (kg/m^2)
2286 : !autodocument_end
2287 :
2288 : ! local variables
2289 :
2290 : integer (kind=int_kind) :: &
2291 : n ! category index
2292 :
2293 : real (kind=dbl_kind) :: &
2294 509328 : worka, workb ! temporary variables
2295 :
2296 : ! 2D coupler variables (computed for each category, then aggregated)
2297 : real (kind=dbl_kind) :: &
2298 509328 : fswabsn , & ! shortwave absorbed by ice (W/m^2)
2299 509328 : flwoutn , & ! upward LW at surface (W/m^2)
2300 509328 : evapn , & ! flux of vapor, atmos to ice (kg m-2 s-1)
2301 509328 : evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1)
2302 509328 : evapin , & ! flux of vapor, atmos to ice over ice (kg m-2 s-1)
2303 509328 : freshn , & ! flux of water, ice to ocean (kg/m^2/s)
2304 509328 : fsaltn , & ! flux of salt, ice to ocean (kg/m^2/s)
2305 509328 : fhocnn , & ! fbot corrected for leftover energy (W/m^2)
2306 509328 : strairxn , & ! air/ice zonal stress, (N/m^2)
2307 509328 : strairyn , & ! air/ice meridional stress, (N/m^2)
2308 509328 : Cdn_atm_ratio_n, & ! drag coefficient ratio
2309 509328 : Trefn , & ! air tmp reference level (K)
2310 509328 : Urefn , & ! air speed reference level (m/s)
2311 509328 : Qrefn , & ! air sp hum reference level (kg/kg)
2312 509328 : shcoef , & ! transfer coefficient for sensible heat
2313 509328 : lhcoef , & ! transfer coefficient for latent heat
2314 509328 : rfrac ! water fraction retained for melt ponds
2315 :
2316 : real (kind=dbl_kind), dimension(n_iso) :: &
2317 2260080 : Qrefn_iso , & ! isotope air sp hum reference level (kg/kg)
2318 2260080 : fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s)
2319 3278736 : fiso_evapn ! isotope evaporation (kg/m^2/s)
2320 :
2321 : real (kind=dbl_kind), allocatable, dimension(:,:) :: &
2322 1384704 : l_isosno , & ! local snow isotope tracer (kg/m^2)
2323 1384704 : l_isoice ! local ice isotope tracer (kg/m^2)
2324 :
2325 : real (kind=dbl_kind), allocatable, dimension(:) :: &
2326 1384704 : l_Qa_iso , & ! local isotope specific humidity (kg/kg)
2327 1384704 : l_Qref_iso , & ! local isotope 2m atm reference spec humidity (kg/kg)
2328 1384704 : l_fiso_atm , & ! local isotope deposition rate (kg/m^2 s)
2329 1384704 : l_fiso_ocn , & ! local isotope flux to ocean (kg/m^2/s)
2330 1384704 : l_fiso_evap ! local isotope evaporation (kg/m^2/s)
2331 :
2332 : real (kind=dbl_kind) :: &
2333 509328 : l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg)
2334 509328 : l_H2_16O_ocn, & ! local ocean concentration of H2_16O (kg/kg)
2335 509328 : l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg)
2336 :
2337 : real (kind=dbl_kind) :: &
2338 509328 : l_fswthru_vdr , & ! vis dir SW through ice to ocean (W/m^2)
2339 509328 : l_fswthru_vdf , & ! vis dif SW through ice to ocean (W/m^2)
2340 509328 : l_fswthru_idr , & ! nir dir SW through ice to ocean (W/m^2)
2341 509328 : l_fswthru_idf ! nir dif SW through ice to ocean (W/m^2)
2342 :
2343 : real (kind=dbl_kind), dimension(:), allocatable :: &
2344 1384704 : l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
2345 1384704 : l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
2346 1384704 : l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
2347 1384704 : l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
2348 :
2349 : real (kind=dbl_kind) :: &
2350 509328 : pond ! water retained in ponds (m)
2351 :
2352 : character(len=*),parameter :: subname='(icepack_step_therm1)'
2353 :
2354 : !-----------------------------------------------------------------
2355 : ! allocate local optional arguments
2356 : !-----------------------------------------------------------------
2357 :
2358 1384704 : if (present(isosno) ) then
2359 1384704 : allocate(l_isosno(size(isosno,dim=1),size(isosno,dim=2)))
2360 8844432 : l_isosno = isosno
2361 : else
2362 0 : allocate(l_isosno(1,1))
2363 0 : l_isosno = c0
2364 : endif
2365 :
2366 1384704 : if (present(isoice) ) then
2367 1384704 : allocate(l_isoice(size(isoice,dim=1),size(isoice,dim=2)))
2368 8844432 : l_isoice = isoice
2369 : else
2370 0 : allocate(l_isoice(1,1))
2371 0 : l_isoice = c0
2372 : endif
2373 :
2374 1384704 : if (present(Qa_iso) ) then
2375 1384704 : allocate(l_Qa_iso(size(Qa_iso)))
2376 5538816 : l_Qa_iso = Qa_iso
2377 : else
2378 0 : allocate(l_Qa_iso(1))
2379 0 : l_Qa_iso = c0
2380 : endif
2381 :
2382 1384704 : if (present(Qref_iso) ) then
2383 1384704 : allocate(l_Qref_iso(size(Qref_iso)))
2384 5538816 : l_Qref_iso = Qref_iso
2385 : else
2386 0 : allocate(l_Qref_iso(1))
2387 0 : l_Qref_iso = c0
2388 : endif
2389 :
2390 1384704 : if (present(fiso_atm) ) then
2391 1384704 : allocate(l_fiso_atm(size(fiso_atm)))
2392 5538816 : l_fiso_atm = fiso_atm
2393 : else
2394 0 : allocate(l_fiso_atm(1))
2395 0 : l_fiso_atm = c0
2396 : endif
2397 :
2398 1384704 : if (present(fiso_ocn) ) then
2399 1384704 : allocate(l_fiso_ocn(size(fiso_ocn)))
2400 5538816 : l_fiso_ocn = fiso_ocn
2401 : else
2402 0 : allocate(l_fiso_ocn(1))
2403 0 : l_fiso_ocn = c0
2404 : endif
2405 :
2406 1384704 : if (present(fiso_evap) ) then
2407 1384704 : allocate(l_fiso_evap(size(fiso_evap)))
2408 5538816 : l_fiso_evap = fiso_evap
2409 : else
2410 0 : allocate(l_fiso_evap(1))
2411 0 : l_fiso_evap = c0
2412 : endif
2413 :
2414 1384704 : l_HDO_ocn = c0
2415 1384704 : if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn
2416 :
2417 1384704 : l_H2_16O_ocn = c0
2418 1384704 : if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
2419 :
2420 1384704 : l_H2_18O_ocn = c0
2421 1384704 : if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
2422 :
2423 1384704 : l_fswthru_vdr = c0
2424 1384704 : if (present(fswthru_vdr) ) l_fswthru_vdr = fswthru_vdr
2425 :
2426 1384704 : l_fswthru_vdf = c0
2427 1384704 : if (present(fswthru_vdf) ) l_fswthru_vdf = fswthru_vdf
2428 :
2429 1384704 : l_fswthru_idr = c0
2430 1384704 : if (present(fswthru_idr) ) l_fswthru_idr = fswthru_idr
2431 :
2432 1384704 : l_fswthru_idf = c0
2433 1384704 : if (present(fswthru_idf) ) l_fswthru_idf = fswthru_idf
2434 :
2435 1384704 : allocate(l_fswthrun_vdr(ncat))
2436 1384704 : allocate(l_fswthrun_vdf(ncat))
2437 1384704 : allocate(l_fswthrun_idr(ncat))
2438 1384704 : allocate(l_fswthrun_idf(ncat))
2439 :
2440 7922112 : l_fswthrun_vdr = c0
2441 7922112 : if (present(fswthrun_vdr) ) l_fswthrun_vdr = fswthrun_vdr
2442 :
2443 7922112 : l_fswthrun_vdf = c0
2444 7922112 : if (present(fswthrun_vdf) ) l_fswthrun_vdf = fswthrun_vdf
2445 :
2446 7922112 : l_fswthrun_idr = c0
2447 7922112 : if (present(fswthrun_idr) ) l_fswthrun_idr = fswthrun_idr
2448 :
2449 7922112 : l_fswthrun_idf = c0
2450 7922112 : if (present(fswthrun_idf) ) l_fswthrun_idf = fswthrun_idf
2451 :
2452 : !-----------------------------------------------------------------
2453 : ! Adjust frzmlt to account for ice-ocean heat fluxes since last
2454 : ! call to coupler.
2455 : ! Compute lateral and bottom heat fluxes.
2456 : !-----------------------------------------------------------------
2457 :
2458 : call frzmlt_bottom_lateral (dt, ncat, &
2459 : nilyr, nslyr, &
2460 : aice, frzmlt, &
2461 0 : vicen, vsnon, &
2462 0 : zqin, zqsn, &
2463 : sst, Tf, &
2464 : ustar_min, &
2465 : fbot_xfer_type, &
2466 : strocnxT, strocnyT, &
2467 : Tbot, fbot, &
2468 : rside, Cdn_ocn, &
2469 1384704 : fside)
2470 :
2471 1384704 : if (icepack_warnings_aborted(subname)) return
2472 :
2473 : !-----------------------------------------------------------------
2474 : ! Update the neutral drag coefficients to account for form drag
2475 : ! Oceanic and atmospheric drag coefficients
2476 : !-----------------------------------------------------------------
2477 :
2478 1384704 : if (formdrag) then
2479 0 : call neutral_drag_coeffs (apnd , &
2480 0 : hpnd , ipnd , &
2481 0 : alvl , vlvl , &
2482 : aice , vice, &
2483 0 : vsno , aicen , &
2484 0 : vicen , &
2485 : Cdn_ocn , Cdn_ocn_skin, &
2486 : Cdn_ocn_floe, Cdn_ocn_keel, &
2487 : Cdn_atm , Cdn_atm_skin, &
2488 : Cdn_atm_floe, Cdn_atm_pond, &
2489 : Cdn_atm_rdg , hfreebd , &
2490 : hdraft , hridge , &
2491 : distrdg , hkeel , &
2492 : dkeel , lfloe , &
2493 96528 : dfloe , ncat)
2494 96528 : if (icepack_warnings_aborted(subname)) return
2495 : endif
2496 :
2497 7922112 : do n = 1, ncat
2498 :
2499 6537408 : meltsn (n) = c0
2500 6537408 : melttn (n) = c0
2501 6537408 : meltbn (n) = c0
2502 6537408 : congeln(n) = c0
2503 6537408 : snoicen(n) = c0
2504 6537408 : dsnown (n) = c0
2505 :
2506 6537408 : Trefn = c0
2507 6537408 : Qrefn = c0
2508 7459728 : Qrefn_iso(:) = c0
2509 7459728 : fiso_ocnn(:) = c0
2510 7459728 : fiso_evapn(:) = c0
2511 6537408 : Urefn = c0
2512 6537408 : lhcoef = c0
2513 6537408 : shcoef = c0
2514 6537408 : worka = c0
2515 6537408 : workb = c0
2516 :
2517 6537408 : if (aicen_init(n) > puny) then
2518 :
2519 4420832 : if (calc_Tsfc .or. calc_strair) then
2520 :
2521 : !-----------------------------------------------------------------
2522 : ! Atmosphere boundary layer calculation; compute coefficients
2523 : ! for sensible and latent heat fluxes.
2524 : !
2525 : ! NOTE: The wind stress is computed here for later use if
2526 : ! calc_strair = .true. Otherwise, the wind stress
2527 : ! components are set to the data values.
2528 : !-----------------------------------------------------------------
2529 :
2530 : call icepack_atm_boundary('ice', &
2531 1617332 : Tsfc(n), potT, &
2532 : uatm, vatm, &
2533 : wind, zlvl, &
2534 : Qa, rhoa, &
2535 : strairxn, strairyn, &
2536 : Trefn, Qrefn, &
2537 : worka, workb, &
2538 : lhcoef, shcoef, &
2539 : Cdn_atm, &
2540 : Cdn_atm_ratio_n, &
2541 : Qa_iso=l_Qa_iso, &
2542 0 : Qref_iso=Qrefn_iso, &
2543 : uvel=uvel, vvel=vvel, &
2544 4420832 : Uref=Urefn)
2545 4420832 : if (icepack_warnings_aborted(subname)) return
2546 :
2547 : endif ! calc_Tsfc or calc_strair
2548 :
2549 4420832 : if (.not.(calc_strair)) then
2550 : #ifndef CICE_IN_NEMO
2551 : ! Set to data values (on T points)
2552 0 : strairxn = strax
2553 0 : strairyn = stray
2554 : #else
2555 : ! NEMO wind stress is supplied on u grid, multipied
2556 : ! by ice concentration and set directly in evp, so
2557 : ! strairxT/yT = 0. Zero u-components here for safety.
2558 : strairxn = c0
2559 : strairyn = c0
2560 : #endif
2561 : endif
2562 :
2563 : !-----------------------------------------------------------------
2564 : ! Update ice age
2565 : ! This is further adjusted for freezing in the thermodynamics.
2566 : ! Melting does not alter the ice age.
2567 : !-----------------------------------------------------------------
2568 :
2569 4420832 : if (tr_iage) call increment_age (dt, iage(n))
2570 4420832 : if (icepack_warnings_aborted(subname)) return
2571 4420832 : if (tr_FY) call update_FYarea (dt, &
2572 : lmask_n, lmask_s, &
2573 229760 : yday, FY(n))
2574 4420832 : if (icepack_warnings_aborted(subname)) return
2575 :
2576 : !-----------------------------------------------------------------
2577 : ! Vertical thermodynamics: Heat conduction, growth and melting.
2578 : !-----------------------------------------------------------------
2579 :
2580 4420832 : if (.not.(calc_Tsfc)) then
2581 :
2582 : ! If not calculating surface temperature and fluxes, set
2583 : ! surface fluxes (flatn, fsurfn, and fcondtopn) to be used
2584 : ! in thickness_changes
2585 :
2586 : ! hadgem routine sets fluxes to default values in ice-only mode
2587 0 : call set_sfcflux(aicen (n), &
2588 0 : flatn_f (n), fsensn_f (n), &
2589 0 : fcondtopn_f(n), &
2590 0 : fsurfn_f (n), &
2591 0 : flatn (n), fsensn (n), &
2592 0 : fsurfn (n), &
2593 0 : fcondtopn (n))
2594 0 : if (icepack_warnings_aborted(subname)) return
2595 : endif
2596 :
2597 : call thermo_vertical(nilyr, nslyr, &
2598 1617332 : dt, aicen (n), &
2599 1617332 : vicen (n), vsnon (n), &
2600 1617332 : Tsfc (n), zSin (:,n), &
2601 0 : zqin (:,n), zqsn (:,n), &
2602 1617332 : apnd (n), hpnd (n), &
2603 : tr_pond_topo, &
2604 : flw, potT, &
2605 : Qa, rhoa, &
2606 : fsnow, fpond, &
2607 : fbot, Tbot, &
2608 : Tsnice, sss, &
2609 : lhcoef, shcoef, &
2610 1617332 : fswsfcn (n), fswintn (n), &
2611 0 : Sswabsn(:,n), Iswabsn(:,n), &
2612 1617332 : fsurfn (n), fcondtopn(n), &
2613 1617332 : fcondbotn(n), &
2614 1617332 : fsensn (n), flatn (n), &
2615 : flwoutn, evapn, &
2616 : evapsn, evapin, &
2617 : freshn, fsaltn, &
2618 : fhocnn, &
2619 1617332 : melttn (n), meltsn (n), &
2620 1617332 : meltbn (n), &
2621 1617332 : congeln (n), snoicen (n), &
2622 : mlt_onset, frz_onset, &
2623 1617332 : yday, dsnown (n), &
2624 7655496 : prescribed_ice)
2625 :
2626 4420832 : if (icepack_warnings_aborted(subname)) then
2627 0 : call icepack_warnings_add(subname//' ice: Vertical thermo error: ')
2628 0 : return
2629 : endif
2630 :
2631 : !-----------------------------------------------------------------
2632 : ! Total absorbed shortwave radiation
2633 : !-----------------------------------------------------------------
2634 :
2635 4420832 : fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)
2636 :
2637 : !-----------------------------------------------------------------
2638 : ! Aerosol update
2639 : !-----------------------------------------------------------------
2640 :
2641 4420832 : if (tr_aero) then
2642 : call update_aerosol (dt, &
2643 : nilyr, nslyr, n_aero, &
2644 0 : melttn (n), meltsn (n), &
2645 0 : meltbn (n), congeln (n), &
2646 0 : snoicen (n), fsnow, &
2647 0 : aerosno(:,:,n), aeroice(:,:,n), &
2648 0 : aicen_init (n), vicen_init (n), &
2649 0 : vsnon_init (n), &
2650 0 : vicen (n), vsnon (n), &
2651 0 : aicen (n), &
2652 229765 : faero_atm , faero_ocn)
2653 229765 : if (icepack_warnings_aborted(subname)) return
2654 : endif
2655 :
2656 4420832 : if (tr_iso) then
2657 : call update_isotope (dt = dt, &
2658 : nilyr = nilyr, nslyr = nslyr, &
2659 0 : meltt = melttn(n),melts = meltsn(n), &
2660 0 : meltb = meltbn(n),congel=congeln(n), &
2661 0 : snoice=snoicen(n),evap=evapn, &
2662 0 : fsnow=fsnow, Tsfc=Tsfc(n), &
2663 0 : Qref_iso=Qrefn_iso(:), &
2664 0 : isosno=l_isosno(:,n),isoice=l_isoice(:,n), &
2665 0 : aice_old=aicen_init(n),vice_old=vicen_init(n), &
2666 0 : vsno_old=vsnon_init(n), &
2667 0 : vicen=vicen(n),vsnon=vsnon(n), &
2668 0 : aicen=aicen(n), &
2669 : fiso_atm=l_fiso_atm(:), &
2670 0 : fiso_evapn=fiso_evapn(:), &
2671 0 : fiso_ocnn=fiso_ocnn(:), &
2672 : HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, &
2673 229760 : H2_18O_ocn=l_H2_18O_ocn)
2674 229760 : if (icepack_warnings_aborted(subname)) return
2675 : endif
2676 : endif ! aicen_init
2677 :
2678 : !-----------------------------------------------------------------
2679 : ! Melt ponds
2680 : ! If using tr_pond_cesm, the full calculation is performed here.
2681 : ! If using tr_pond_topo, the rest of the calculation is done after
2682 : ! the surface fluxes are merged, below.
2683 : !-----------------------------------------------------------------
2684 :
2685 : !call ice_timer_start(timer_ponds)
2686 6537408 : if (tr_pond) then
2687 :
2688 5475600 : if (tr_pond_cesm) then
2689 307440 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2690 : call compute_ponds_cesm(dt, hi_min, &
2691 : pndaspect, rfrac, &
2692 0 : melttn(n), meltsn(n), &
2693 : frain, &
2694 0 : aicen (n), vicen (n), &
2695 0 : Tsfc (n), &
2696 307440 : apnd (n), hpnd (n))
2697 307440 : if (icepack_warnings_aborted(subname)) return
2698 :
2699 5168160 : elseif (tr_pond_lvl) then
2700 4378080 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2701 : call compute_ponds_lvl(dt, nilyr, &
2702 : ktherm, &
2703 : hi_min, &
2704 : dpscale, frzpnd, &
2705 : pndaspect, rfrac, &
2706 1845840 : melttn(n), meltsn(n), &
2707 : frain, Tair, &
2708 1845840 : fsurfn(n), &
2709 1845840 : dhsn (n), ffracn(n), &
2710 1845840 : aicen (n), vicen (n), &
2711 1845840 : vsnon (n), &
2712 0 : zqin(:,n), zSin(:,n), &
2713 1845840 : Tsfc (n), alvl (n), &
2714 1845840 : apnd (n), hpnd (n), &
2715 6223920 : ipnd (n))
2716 4378080 : if (icepack_warnings_aborted(subname)) return
2717 :
2718 790080 : elseif (tr_pond_topo) then
2719 790080 : if (aicen_init(n) > puny) then
2720 :
2721 : ! collect liquid water in ponds
2722 : ! assume salt still runs off
2723 457021 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2724 80554 : pond = rfrac/rhofresh * (melttn(n)*rhoi &
2725 80554 : + meltsn(n)*rhos &
2726 457021 : + frain *dt)
2727 :
2728 : ! if pond does not exist, create new pond over full ice area
2729 : ! otherwise increase pond depth without changing pond area
2730 457021 : if (apnd(n) < puny) then
2731 272114 : hpnd(n) = c0
2732 272114 : apnd(n) = c1
2733 : endif
2734 457021 : hpnd(n) = (pond + hpnd(n)*apnd(n)) / apnd(n)
2735 457021 : fpond = fpond + pond * aicen(n) ! m
2736 : endif ! aicen_init
2737 : endif
2738 :
2739 : endif ! tr_pond
2740 : !call ice_timer_stop(timer_ponds)
2741 :
2742 : !-----------------------------------------------------------------
2743 : ! Increment area-weighted fluxes.
2744 : !-----------------------------------------------------------------
2745 :
2746 6537408 : if (aicen_init(n) > puny) &
2747 1617332 : call merge_fluxes (aicen=aicen_init(n), &
2748 : flw=flw, &
2749 : strairxn=strairxn, strairyn=strairyn,&
2750 : Cdn_atm_ratio_n=Cdn_atm_ratio_n, &
2751 1617332 : fsurfn=fsurfn(n), fcondtopn=fcondtopn(n),&
2752 1617332 : fcondbotn=fcondbotn(n), &
2753 1617332 : fsensn=fsensn(n), flatn=flatn(n), &
2754 : fswabsn=fswabsn, flwoutn=flwoutn, &
2755 : evapn=evapn, &
2756 : evapsn=evapsn, evapin=evapin, &
2757 : Trefn=Trefn, Qrefn=Qrefn, &
2758 : freshn=freshn, fsaltn=fsaltn, &
2759 : fhocnn=fhocnn, &
2760 1617332 : fswthrun=fswthrun(n), &
2761 0 : fswthrun_vdr=l_fswthrun_vdr(n), &
2762 0 : fswthrun_vdf=l_fswthrun_vdf(n), &
2763 0 : fswthrun_idr=l_fswthrun_idr(n), &
2764 0 : fswthrun_idf=l_fswthrun_idf(n), &
2765 : strairxT=strairxT, strairyT=strairyT,&
2766 : Cdn_atm_ratio=Cdn_atm_ratio, &
2767 : fsurf=fsurf, fcondtop=fcondtop,&
2768 : fcondbot=fcondbot, &
2769 : fsens=fsens, flat=flat, &
2770 : fswabs=fswabs, flwout=flwout, &
2771 : evap=evap, &
2772 : evaps=evaps, evapi=evapi, &
2773 : Tref=Tref, Qref=Qref, &
2774 : fresh=fresh, fsalt=fsalt, &
2775 : fhocn=fhocn, &
2776 : fswthru=fswthru, &
2777 : fswthru_vdr=l_fswthru_vdr, &
2778 : fswthru_vdf=l_fswthru_vdf, &
2779 : fswthru_idr=l_fswthru_idr, &
2780 : fswthru_idf=l_fswthru_idf, &
2781 1617332 : melttn=melttn (n), meltsn=meltsn(n), &
2782 1617332 : meltbn=meltbn (n), congeln=congeln(n),&
2783 1617332 : snoicen=snoicen(n), &
2784 : meltt=meltt, melts=melts, &
2785 : meltb=meltb, congel=congel, &
2786 : snoice=snoice, &
2787 : Uref=Uref, Urefn=Urefn, &
2788 : Qref_iso=l_Qref_iso, &
2789 0 : Qrefn_iso=Qrefn_iso, &
2790 : fiso_ocn=l_fiso_ocn, &
2791 0 : fiso_ocnn=fiso_ocnn, &
2792 : fiso_evap=l_fiso_evap, &
2793 6038164 : fiso_evapn=fiso_evapn)
2794 :
2795 7922112 : if (icepack_warnings_aborted(subname)) return
2796 :
2797 : enddo ! ncat
2798 :
2799 8844432 : if (present(isosno) ) isosno = l_isosno
2800 8844432 : if (present(isoice) ) isoice = l_isoice
2801 5538816 : if (present(Qa_iso) ) Qa_iso = l_Qa_iso
2802 5538816 : if (present(Qref_iso) ) Qref_iso = l_Qref_iso
2803 5538816 : if (present(fiso_atm) ) fiso_atm = l_fiso_atm
2804 5538816 : if (present(fiso_ocn) ) fiso_ocn = l_fiso_ocn
2805 5538816 : if (present(fiso_evap)) fiso_evap= l_fiso_evap
2806 7922112 : if (present(fswthrun_vdr)) fswthrun_vdr= l_fswthrun_vdr
2807 7922112 : if (present(fswthrun_vdf)) fswthrun_vdf= l_fswthrun_vdf
2808 7922112 : if (present(fswthrun_idr)) fswthrun_idr= l_fswthrun_idr
2809 7922112 : if (present(fswthrun_idf)) fswthrun_idf= l_fswthrun_idf
2810 1384704 : if (present(fswthru_vdr)) fswthru_vdr= l_fswthru_vdr
2811 1384704 : if (present(fswthru_vdf)) fswthru_vdf= l_fswthru_vdf
2812 1384704 : if (present(fswthru_idr)) fswthru_idr= l_fswthru_idr
2813 1384704 : if (present(fswthru_idf)) fswthru_idf= l_fswthru_idf
2814 1384704 : deallocate(l_isosno)
2815 1384704 : deallocate(l_isoice)
2816 1384704 : deallocate(l_Qa_iso)
2817 1384704 : deallocate(l_Qref_iso)
2818 1384704 : deallocate(l_fiso_atm)
2819 1384704 : deallocate(l_fiso_ocn)
2820 1384704 : deallocate(l_fiso_evap)
2821 1384704 : deallocate(l_fswthrun_vdr)
2822 1384704 : deallocate(l_fswthrun_vdf)
2823 1384704 : deallocate(l_fswthrun_idr)
2824 1384704 : deallocate(l_fswthrun_idf)
2825 :
2826 : !-----------------------------------------------------------------
2827 : ! Calculate ponds from the topographic scheme
2828 : !-----------------------------------------------------------------
2829 : !call ice_timer_start(timer_ponds)
2830 1384704 : if (tr_pond_topo) then
2831 : call compute_ponds_topo(dt, ncat, nilyr, &
2832 : ktherm, heat_capacity, &
2833 0 : aice, aicen, &
2834 0 : vice, vicen, &
2835 0 : vsno, vsnon, &
2836 : meltt, &
2837 : fsurf, fpond, &
2838 0 : Tsfc, Tf, &
2839 0 : zqin, zSin, &
2840 158016 : apnd, hpnd, ipnd )
2841 158016 : if (icepack_warnings_aborted(subname)) return
2842 : endif
2843 : !call ice_timer_stop(timer_ponds)
2844 :
2845 1384704 : end subroutine icepack_step_therm1
2846 :
2847 : !=======================================================================
2848 :
2849 : end module icepack_therm_vertical
2850 :
2851 : !=======================================================================
|
