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