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 1599186245 : subroutine thermo_vertical (dt, aicen, &
83 : vicen, vsnon, & ! LCOV_EXCL_LINE
84 0 : Tsf, zSin, & ! LCOV_EXCL_LINE
85 1599186245 : zqin, zqsn, & ! LCOV_EXCL_LINE
86 : apond, hpond, & ! LCOV_EXCL_LINE
87 : flw, potT, & ! LCOV_EXCL_LINE
88 : Qa, rhoa, & ! LCOV_EXCL_LINE
89 : fsnow, fpond, & ! LCOV_EXCL_LINE
90 : fbot, Tbot, & ! LCOV_EXCL_LINE
91 : Tsnice, sss, & ! LCOV_EXCL_LINE
92 0 : sst, rsnw, & ! LCOV_EXCL_LINE
93 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
94 : fswsfc, fswint, & ! LCOV_EXCL_LINE
95 1599186245 : Sswabs, Iswabs, & ! LCOV_EXCL_LINE
96 : fsurfn, fcondtopn, & ! LCOV_EXCL_LINE
97 : fcondbotn, & ! LCOV_EXCL_LINE
98 : fsensn, flatn, & ! LCOV_EXCL_LINE
99 : flwoutn, evapn, & ! LCOV_EXCL_LINE
100 : evapsn, evapin, & ! LCOV_EXCL_LINE
101 : freshn, fsaltn, & ! LCOV_EXCL_LINE
102 : fhocnn, frain, & ! LCOV_EXCL_LINE
103 : meltt, melts, & ! LCOV_EXCL_LINE
104 : meltb, meltsliq, & ! LCOV_EXCL_LINE
105 3198372490 : smice, massice, & ! LCOV_EXCL_LINE
106 3198372490 : smliq, massliq, & ! LCOV_EXCL_LINE
107 : congel, snoice, & ! LCOV_EXCL_LINE
108 : mlt_onset, frz_onset, & ! LCOV_EXCL_LINE
109 : yday, dsnow, & ! LCOV_EXCL_LINE
110 : prescribed_ice)
111 :
112 : real (kind=dbl_kind), intent(in) :: &
113 : dt , & ! time step ! LCOV_EXCL_LINE
114 : frain ! rainfall rate (kg/m2/s)
115 :
116 : ! ice state variables
117 : real (kind=dbl_kind), intent(inout) :: &
118 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
119 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
125 : apond , & ! melt pond area fraction of category ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
134 : zqin , & ! ice layer enthalpy, zqin < 0 (J m-3) ! LCOV_EXCL_LINE
135 : zSin , & ! internal ice layer salinities ! LCOV_EXCL_LINE
136 : rsnw , & ! snow grain radius (10^-6 m) ! LCOV_EXCL_LINE
137 : smice , & ! ice mass tracer in snow (kg/m^3) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
147 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
148 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
149 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
150 : fsnow , & ! snowfall rate (kg m-2 s-1) ! LCOV_EXCL_LINE
151 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
156 : fswint , & ! SW absorbed in ice interior, below surface (W m-2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
166 : Tbot , & ! ice bottom surface temperature (deg C) ! LCOV_EXCL_LINE
167 : sst , & ! sea surface temperature (C) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
173 : evapn , & ! evaporative water flux (kg/m^2/s) ! LCOV_EXCL_LINE
174 : evapsn , & ! evaporative water flux over snow (kg/m^2/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
180 : flatn , & ! latent heat flux (W/m^2) ! LCOV_EXCL_LINE
181 : fsurfn , & ! net flux to top surface, excluding fcondtopn ! LCOV_EXCL_LINE
182 : fcondtopn, & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
188 : fsaltn , & ! salt flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
194 : meltt , & ! top ice melt (m/step-->cm/day) ! LCOV_EXCL_LINE
195 : melts , & ! snow melt (m/step-->cm/day) ! LCOV_EXCL_LINE
196 : meltsliq , & ! snow melt mass (kg/m^2/step-->kg/m^2/day) ! LCOV_EXCL_LINE
197 : meltb , & ! basal ice melt (m/step-->cm/day) ! LCOV_EXCL_LINE
198 : congel , & ! basal ice growth (m/step-->cm/day) ! LCOV_EXCL_LINE
199 : snoice , & ! snow-ice formation (m/step-->cm/day) ! LCOV_EXCL_LINE
200 : dsnow , & ! change in snow thickness (m/step-->cm/day) ! LCOV_EXCL_LINE
201 : mlt_onset, & ! day of year that sfc melting begins ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
214 : dhs ! change in snow thickness
215 :
216 : ! 2D state variables (thickness, temperature)
217 :
218 : real (kind=dbl_kind) :: &
219 : hilyr , & ! ice layer thickness ! LCOV_EXCL_LINE
220 : hslyr , & ! snow layer thickness ! LCOV_EXCL_LINE
221 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
222 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
223 : hsn_new , & ! thickness of new snow (m) ! LCOV_EXCL_LINE
224 : worki , & ! local work array ! LCOV_EXCL_LINE
225 : works ! local work array
226 :
227 : real (kind=dbl_kind), dimension (nilyr) :: &
228 1599186245 : zTin ! internal ice layer temperatures
229 :
230 : real (kind=dbl_kind), dimension (nslyr) :: &
231 1599186245 : 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) ! LCOV_EXCL_LINE
237 : efinal , & ! final energy of melting (J m-2) ! LCOV_EXCL_LINE
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 1599186245 : flwoutn = c0
250 1599186245 : evapn = c0
251 1599186245 : evapsn = c0
252 1599186245 : evapin = c0
253 1599186245 : freshn = c0
254 1599186245 : fsaltn = c0
255 1599186245 : fhocnn = c0
256 1599186245 : fadvocn = c0
257 :
258 1599186245 : meltt = c0
259 1599186245 : meltb = c0
260 1599186245 : melts = c0
261 1599186245 : congel = c0
262 1599186245 : snoice = c0
263 1599186245 : dsnow = c0
264 3381641432 : zTsn(:) = c0
265 11891141500 : zTin(:) = c0
266 1599186245 : meltsliq= c0
267 3381641432 : massice(:) = c0
268 3381641432 : massliq(:) = c0
269 :
270 1599186245 : if (calc_Tsfc) then
271 1575092498 : fsensn = c0
272 1575092498 : flatn = c0
273 1575092498 : fsurfn = c0
274 1575092498 : 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, & ! LCOV_EXCL_LINE
283 : hin, hilyr, & ! LCOV_EXCL_LINE
284 : hsn, hslyr, & ! LCOV_EXCL_LINE
285 : zqin, zTin, & ! LCOV_EXCL_LINE
286 : zqsn, zTsn, & ! LCOV_EXCL_LINE
287 : zSin, & ! LCOV_EXCL_LINE
288 1599186245 : einit )
289 1599186245 : if (icepack_warnings_aborted(subname)) return
290 :
291 : ! Save initial ice and snow thickness (for fresh and fsalt)
292 1599186245 : worki = hin
293 1599186245 : works = hsn
294 :
295 : ! Save initial salt volume for prognostic flux
296 1599186245 : if (saltflux_option == 'prognostic') then
297 0 : saltvol = c0
298 0 : do k=1,nilyr
299 0 : 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 1599186245 : if (ktherm == 2) then
309 :
310 : call temperature_changes_salinity(dt, &
311 : rhoa, flw, & ! LCOV_EXCL_LINE
312 : potT, Qa, & ! LCOV_EXCL_LINE
313 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
314 : fswsfc, fswint, & ! LCOV_EXCL_LINE
315 : Sswabs, Iswabs, & ! LCOV_EXCL_LINE
316 : hilyr, hslyr, & ! LCOV_EXCL_LINE
317 : apond, hpond, & ! LCOV_EXCL_LINE
318 : zqin, zTin, & ! LCOV_EXCL_LINE
319 : zqsn, zTsn, & ! LCOV_EXCL_LINE
320 : zSin, & ! LCOV_EXCL_LINE
321 : Tsf, Tbot, & ! LCOV_EXCL_LINE
322 : sss, & ! LCOV_EXCL_LINE
323 : fsensn, flatn, & ! LCOV_EXCL_LINE
324 : flwoutn, fsurfn, & ! LCOV_EXCL_LINE
325 : fcondtopn, fcondbotn, & ! LCOV_EXCL_LINE
326 : fadvocn, snoice, & ! LCOV_EXCL_LINE
327 1408028216 : smice, smliq)
328 1408028216 : if (icepack_warnings_aborted(subname)) return
329 :
330 : else ! ktherm
331 :
332 : call temperature_changes(dt, &
333 : rhoa, flw, & ! LCOV_EXCL_LINE
334 : potT, Qa, & ! LCOV_EXCL_LINE
335 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
336 : fswsfc, fswint, & ! LCOV_EXCL_LINE
337 : Sswabs, Iswabs, & ! LCOV_EXCL_LINE
338 : hilyr, hslyr, & ! LCOV_EXCL_LINE
339 : zqin, zTin, & ! LCOV_EXCL_LINE
340 : zqsn, zTsn, & ! LCOV_EXCL_LINE
341 : zSin, & ! LCOV_EXCL_LINE
342 : Tsf, Tbot, & ! LCOV_EXCL_LINE
343 : fsensn, flatn, & ! LCOV_EXCL_LINE
344 : flwoutn, fsurfn, & ! LCOV_EXCL_LINE
345 : fcondtopn, fcondbotn, & ! LCOV_EXCL_LINE
346 191158029 : einit )
347 191158029 : if (icepack_warnings_aborted(subname)) return
348 :
349 : endif ! ktherm
350 :
351 : ! mass of ice and liquid water in snow
352 1599186245 : if (snwgrain) then
353 203601768 : massice(:) = smice(:) * hslyr
354 203601768 : massliq(:) = smliq(:) * hslyr
355 : endif
356 :
357 : ! intermediate energy for error check
358 :
359 1599186245 : einter = c0
360 3381641432 : do k = 1, nslyr
361 3381641432 : einter = einter + hslyr * zqsn(k)
362 : enddo ! k
363 11891141500 : do k = 1, nilyr
364 11891141500 : einter = einter + hilyr * zqin(k)
365 : enddo ! k
366 :
367 1599186245 : Tsnice = c0
368 1599186245 : if ((hslyr+hilyr) > puny) then
369 1599186245 : if (hslyr > puny) then
370 1383406330 : Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr)
371 : else
372 215779915 : Tsnice = Tsf
373 : endif
374 : endif
375 :
376 1599186245 : 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, & ! LCOV_EXCL_LINE
386 : hin, hilyr, & ! LCOV_EXCL_LINE
387 : hsn, hslyr, & ! LCOV_EXCL_LINE
388 : zqin, zqsn, & ! LCOV_EXCL_LINE
389 : smice, massice, & ! LCOV_EXCL_LINE
390 : smliq, massliq, & ! LCOV_EXCL_LINE
391 : fbot, Tbot, & ! LCOV_EXCL_LINE
392 : flatn, fsurfn, & ! LCOV_EXCL_LINE
393 : fcondtopn, fcondbotn, & ! LCOV_EXCL_LINE
394 : fsnow, hsn_new, & ! LCOV_EXCL_LINE
395 : fhocnn, evapn, & ! LCOV_EXCL_LINE
396 : evapsn, evapin, & ! LCOV_EXCL_LINE
397 : meltt, melts, & ! LCOV_EXCL_LINE
398 : meltsliq, frain, & ! LCOV_EXCL_LINE
399 : meltb, & ! LCOV_EXCL_LINE
400 : congel, snoice, & ! LCOV_EXCL_LINE
401 : mlt_onset, frz_onset, & ! LCOV_EXCL_LINE
402 : zSin, sss, & ! LCOV_EXCL_LINE
403 : sst, & ! LCOV_EXCL_LINE
404 1599186245 : dsnow, rsnw)
405 1599186245 : 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, & ! LCOV_EXCL_LINE
414 : fhocnn, fswint, & ! LCOV_EXCL_LINE
415 : fsnow, einit, & ! LCOV_EXCL_LINE
416 : einter, efinal, & ! LCOV_EXCL_LINE
417 : fcondtopn, fcondbotn, & ! LCOV_EXCL_LINE
418 1599186245 : fadvocn, fbot )
419 1599186245 : if (icepack_warnings_aborted(subname)) return
420 :
421 : !-----------------------------------------------------------------
422 : ! If prescribed ice, set hi back to old values
423 : !-----------------------------------------------------------------
424 :
425 1599186245 : if (present(prescribed_ice)) then
426 1597131722 : 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 1599186245 : dhi = hin - worki
439 1599186245 : dhs = hsn - works - hsn_new
440 :
441 1599186245 : freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt
442 1599186245 : if (saltflux_option == 'prognostic') then
443 0 : dfsalt = c0
444 0 : do k=1,nilyr
445 0 : dfsalt = dfsalt + rhoi*zSin(k)*hin*p001 / real(nilyr,kind=dbl_kind)
446 : enddo
447 0 : fsaltn = fsaltn - (dfsalt - saltvol) / dt
448 : else
449 1599186245 : fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
450 : endif
451 1599186245 : fhocnn = fhocnn + fadvocn ! for ktherm=2
452 :
453 1599186245 : if (hin == c0) then
454 121640 : 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, & ! LCOV_EXCL_LINE
464 : zqin, zSin, & ! LCOV_EXCL_LINE
465 : zqsn, & ! LCOV_EXCL_LINE
466 : aicen, & ! LCOV_EXCL_LINE
467 1599186245 : vicen, vsnon)
468 1599186245 : 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 1492190040 : subroutine frzmlt_bottom_lateral (dt, &
482 : aice, frzmlt, & ! LCOV_EXCL_LINE
483 0 : vicen, vsnon, & ! LCOV_EXCL_LINE
484 2984380080 : qicen, qsnon, & ! LCOV_EXCL_LINE
485 : sst, Tf, & ! LCOV_EXCL_LINE
486 : ustar_min, & ! LCOV_EXCL_LINE
487 : fbot_xfer_type, & ! LCOV_EXCL_LINE
488 : strocnxT, strocnyT, & ! LCOV_EXCL_LINE
489 : Tbot, fbot, & ! LCOV_EXCL_LINE
490 1492190040 : rsiden, Cdn_ocn, & ! LCOV_EXCL_LINE
491 1492190040 : wlat, aicen, & ! LCOV_EXCL_LINE
492 1492190040 : 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 ! LCOV_EXCL_LINE
499 : frzmlt , & ! freezing/melting potential (W/m^2) ! LCOV_EXCL_LINE
500 : sst , & ! sea surface temperature (C) ! LCOV_EXCL_LINE
501 : Tf , & ! freezing temperature (C) ! LCOV_EXCL_LINE
502 : ustar_min,& ! minimum friction velocity for ice-ocean heat flux ! LCOV_EXCL_LINE
503 : Cdn_ocn , & ! ocean-ice neutral drag coefficient ! LCOV_EXCL_LINE
504 : strocnxT, & ! ice-ocean stress, x-direction ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
512 : vsnon ! snow volume (m)
513 :
514 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
515 : qicen , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
516 : qsnon ! snow layer enthalpy (J m-3)
517 :
518 : real (kind=dbl_kind), intent(out) :: &
519 : Tbot , & ! ice bottom surface temperature (deg C) ! LCOV_EXCL_LINE
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 1492188141 : delta_an , & ! change in the ITD ! LCOV_EXCL_LINE
537 1492190040 : G_radialn ! lateral melt rate in FSD (m/s)
538 :
539 : real (kind=dbl_kind) :: &
540 : rside , & ! ! LCOV_EXCL_LINE
541 : fside ! lateral heat flux (W/m^2)
542 :
543 : integer (kind=int_kind) :: &
544 : n , & ! thickness category index ! LCOV_EXCL_LINE
545 : k ! layer index
546 :
547 : real (kind=dbl_kind) :: &
548 : etot , & ! total energy in column ! LCOV_EXCL_LINE
549 : wlat_loc, & ! lateral melt rate (m/s) ! LCOV_EXCL_LINE
550 : qavg ! average enthalpy in column (approximate)
551 :
552 : real (kind=dbl_kind) :: &
553 : deltaT , & ! SST - Tbot >= 0 ! LCOV_EXCL_LINE
554 : ustar , & ! skin friction velocity for fbot (m/s) ! LCOV_EXCL_LINE
555 : bin1_arealoss, & ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
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 8930659392 : rsiden(:) = c0
578 1492190040 : Tbot = Tf
579 1492190040 : fbot = c0
580 1492190040 : wlat_loc = c0
581 1492190040 : if (present(wlat)) wlat=c0
582 :
583 1492190040 : 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 169279309 : deltaT = max((sst-Tbot),c0)
591 :
592 : ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
593 169279309 : ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow)
594 169279309 : ustar = max (ustar,ustar_min)
595 :
596 169279309 : 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 2309239 : 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 166970070 : cpchr = -cp_ocn*rhow*0.006_dbl_kind
603 : endif
604 :
605 169279309 : fbot = cpchr * deltaT * ustar ! < 0
606 169279309 : 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 169279309 : wlat_loc = m1 * deltaT**m2 ! Maykut & Perovich
618 169279309 : rside = wlat_loc*dt*pi/(floeshape*floediam) ! Steele
619 169279309 : rside = max(c0,min(rside,c1))
620 :
621 169279309 : if (rside == c0) return ! nothing more to do so get out
622 :
623 1016235793 : rsiden(:) = rside
624 :
625 169277410 : if (tr_fsd) then ! alter rsiden now since floes are not of size floediam
626 :
627 11872812 : do n = 1, ncat
628 9894010 : G_radialn(n) = -wlat_loc ! negative
629 :
630 : ! afsdn present check up the calling tree
631 121671560 : 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 9894010 : bin1_arealoss = -afsdn(1,n) / floe_binwidth(1) ! when scaled by *G_radialn(n)*dt*aicen(n)
637 :
638 9894010 : delta_an(n) = c0
639 121671560 : do k = 1, nfsd
640 : ! this is delta_an(n) when scaled by *G_radialn(n)*dt*aicen(n)
641 121671560 : 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 9894010 : delta_an(n) = (delta_an(n) - bin1_arealoss)*G_radialn(n)*dt
646 :
647 9894010 : 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 9894010 : if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n),c1)
654 :
655 11872812 : 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 169277410 : fside = c0
667 1016235793 : do n = 1, ncat
668 :
669 846958383 : etot = c0
670 : ! melting energy/unit area in each column, etot < 0
671 :
672 1754070148 : do k = 1, nslyr
673 1754070148 : etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
674 : enddo
675 :
676 6775667064 : do k = 1, nilyr
677 6775667064 : etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind)
678 : enddo ! nilyr
679 :
680 : ! lateral heat flux, fside < 0
681 1016235793 : 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 169277410 : xtmp = frzmlt/(fbot + fside - puny)
691 169277410 : xtmp = min(xtmp, c1)
692 169277410 : xtmp = max(xtmp, c0)
693 169277410 : fbot = fbot * xtmp
694 :
695 1016235793 : do n = 1, ncat
696 1016235793 : rsiden(n) = rsiden(n) * xtmp ! xtmp is almost always 1 so usually nothing happens here
697 : enddo ! ncat
698 :
699 : endif
700 :
701 1492188141 : 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 1599186245 : subroutine init_vertical_profile(aicen, vicen, &
715 : vsnon, & ! LCOV_EXCL_LINE
716 : hin, hilyr, & ! LCOV_EXCL_LINE
717 : hsn, hslyr, & ! LCOV_EXCL_LINE
718 1599186245 : zqin, zTin, & ! LCOV_EXCL_LINE
719 1599186245 : zqsn, zTsn, & ! LCOV_EXCL_LINE
720 1599186245 : zSin, & ! LCOV_EXCL_LINE
721 : einit )
722 :
723 : real (kind=dbl_kind), intent(in) :: &
724 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
725 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
726 : vsnon ! volume per unit area of snow (m)
727 :
728 : real (kind=dbl_kind), intent(out):: &
729 : hilyr , & ! ice layer thickness ! LCOV_EXCL_LINE
730 : hslyr , & ! snow layer thickness ! LCOV_EXCL_LINE
731 : einit ! initial energy of melting (J m-2)
732 :
733 : real (kind=dbl_kind), intent(out):: &
734 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
735 : hsn ! snow thickness (m)
736 :
737 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
738 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
746 : zTsn ! snow temperature
747 :
748 : ! local variables
749 : real (kind=dbl_kind), dimension(nilyr) :: &
750 1599186245 : Tmlts ! melting temperature
751 :
752 : integer (kind=int_kind) :: &
753 : k ! ice layer index
754 :
755 : real (kind=dbl_kind) :: &
756 : rnslyr, & ! real(nslyr) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
761 : tice_high , & ! flag for zTin > Tmlt ! LCOV_EXCL_LINE
762 : tsno_low , & ! flag for zTsn < Tmin ! LCOV_EXCL_LINE
763 : tice_low ! flag for zTin < Tmin
764 :
765 : character(len=*),parameter :: subname='(init_vertical_profile)'
766 :
767 : !-----------------------------------------------------------------
768 : ! Initialize
769 : !-----------------------------------------------------------------
770 :
771 1599186245 : rnslyr = real(nslyr,kind=dbl_kind)
772 :
773 1599186245 : tsno_high = .false.
774 1599186245 : tice_high = .false.
775 1599186245 : tsno_low = .false.
776 1599186245 : tice_low = .false.
777 :
778 1599186245 : einit = c0
779 :
780 : !-----------------------------------------------------------------
781 : ! Surface temperature, ice and snow thickness
782 : ! Initialize internal energy
783 : !-----------------------------------------------------------------
784 :
785 1599186245 : hin = vicen / aicen
786 1599186245 : hsn = vsnon / aicen
787 1599186245 : hilyr = hin / real(nilyr,kind=dbl_kind)
788 1599186245 : hslyr = hsn / rnslyr
789 :
790 : !-----------------------------------------------------------------
791 : ! Snow enthalpy and maximum allowed snow temperature
792 : !-----------------------------------------------------------------
793 :
794 3381641432 : 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 1782455187 : if (hslyr > hs_min/rnslyr) then
804 : ! zqsn < 0
805 : Tmax = -zqsn(k)*puny*rnslyr / &
806 1455827682 : (rhos*cp_ice*vsnon)
807 : else
808 326627505 : zqsn (k) = -rhos * Lfresh
809 326627505 : Tmax = puny
810 : endif
811 :
812 : !-----------------------------------------------------------------
813 : ! Compute snow temperatures from enthalpies.
814 : ! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0.
815 : !-----------------------------------------------------------------
816 1782455187 : zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice
817 :
818 : !-----------------------------------------------------------------
819 : ! Check for zTsn > Tmax (allowing for roundoff error) and zTsn < Tmin.
820 : !-----------------------------------------------------------------
821 3381641432 : if (zTsn(k) > Tmax) then
822 0 : tsno_high = .true.
823 1782455187 : 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 1599186245 : 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 1599186245 : 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 3381641432 : do k = 1, nslyr
889 :
890 1782455187 : if (zTsn(k) > c0) then ! correct roundoff error
891 15034430 : zTsn(k) = c0
892 15034430 : zqsn(k) = -rhos*Lfresh
893 : endif
894 :
895 : !-----------------------------------------------------------------
896 : ! initial energy per unit area of ice/snow, relative to 0 C
897 : !-----------------------------------------------------------------
898 3381641432 : einit = einit + hslyr*zqsn(k)
899 :
900 : enddo ! nslyr
901 :
902 11891141500 : do k = 1, nilyr
903 :
904 : !---------------------------------------------------------------------
905 : ! Use initial salinity profile for thin ice
906 : !---------------------------------------------------------------------
907 :
908 10291955255 : 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 10291955255 : if (ktherm == 2) then
923 9856197512 : Tmlts(k) = liquidus_temperature_mush(zSin(k))
924 : else
925 435757743 : Tmlts(k) = -zSin(k) * depressT
926 : endif
927 :
928 : !-----------------------------------------------------------------
929 : ! Compute ice temperatures from enthalpies using quadratic formula
930 : !-----------------------------------------------------------------
931 :
932 10291955255 : if (ktherm == 2) then
933 9856197512 : zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k))
934 : else
935 435757743 : zTin(k) = calculate_Tin_from_qin(zqin(k),Tmlts(k))
936 : endif
937 :
938 10291955255 : if (l_brine) then
939 10291955255 : 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 10291955255 : if (zTin(k) > Tmax) then
948 0 : tice_high = .true.
949 10291955255 : 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 10291955255 : 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 10291955255 : 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 10291955255 : if (ktherm /= 2) then
1012 435757743 : if (zTin(k) > c0) then
1013 0 : zTin(k) = c0
1014 0 : zqin(k) = -rhoi*Lfresh
1015 : endif
1016 : endif
1017 :
1018 10291955255 : if (ktherm == 1) then
1019 435757743 : 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 11891141500 : 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 1599186245 : subroutine thickness_changes (dt, yday, &
1044 : efinal, & ! LCOV_EXCL_LINE
1045 : hin, hilyr, & ! LCOV_EXCL_LINE
1046 : hsn, hslyr, & ! LCOV_EXCL_LINE
1047 1599186245 : zqin, zqsn, & ! LCOV_EXCL_LINE
1048 3198372490 : smice, massice, & ! LCOV_EXCL_LINE
1049 3198372490 : smliq, massliq, & ! LCOV_EXCL_LINE
1050 : fbot, Tbot, & ! LCOV_EXCL_LINE
1051 : flatn, fsurfn, & ! LCOV_EXCL_LINE
1052 : fcondtopn, fcondbotn, & ! LCOV_EXCL_LINE
1053 : fsnow, hsn_new, & ! LCOV_EXCL_LINE
1054 : fhocnn, evapn, & ! LCOV_EXCL_LINE
1055 : evapsn, evapin, & ! LCOV_EXCL_LINE
1056 : meltt, melts, & ! LCOV_EXCL_LINE
1057 : meltsliq, frain, & ! LCOV_EXCL_LINE
1058 : meltb, & ! LCOV_EXCL_LINE
1059 : congel, snoice, & ! LCOV_EXCL_LINE
1060 : mlt_onset, frz_onset,& ! LCOV_EXCL_LINE
1061 1599186245 : zSin, sss, & ! LCOV_EXCL_LINE
1062 : sst, & ! LCOV_EXCL_LINE
1063 1599186245 : dsnow, rsnw)
1064 :
1065 : real (kind=dbl_kind), intent(in) :: &
1066 : dt , & ! time step ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1071 : Tbot , & ! ice bottom surface temperature (deg C) ! LCOV_EXCL_LINE
1072 : fsnow , & ! snowfall rate (kg m-2 s-1) ! LCOV_EXCL_LINE
1073 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
1074 : fsurfn , & ! net flux to top surface, excluding fcondtopn ! LCOV_EXCL_LINE
1075 : fcondtopn , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1083 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1084 : rsnw , & ! snow grain radius (10^-6 m) ! LCOV_EXCL_LINE
1085 : smice , & ! ice mass tracer in snow (kg/m^3) ! LCOV_EXCL_LINE
1086 : smliq , & ! liquid water mass tracer in snow (kg/m^3) ! LCOV_EXCL_LINE
1087 : massice , & ! ice mass in snow (kg/m^2) ! LCOV_EXCL_LINE
1088 : massliq ! liquid water mass in snow (kg/m^2)
1089 :
1090 : real (kind=dbl_kind), intent(inout) :: &
1091 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1092 : hslyr ! snow layer thickness (m)
1093 :
1094 : real (kind=dbl_kind), intent(inout) :: &
1095 : meltt , & ! top ice melt (m/step-->cm/day) ! LCOV_EXCL_LINE
1096 : melts , & ! snow melt (m/step-->cm/day) ! LCOV_EXCL_LINE
1097 : meltsliq , & ! snow melt mass (kg/m^2/step-->kg/m^2/day) ! LCOV_EXCL_LINE
1098 : meltb , & ! basal ice melt (m/step-->cm/day) ! LCOV_EXCL_LINE
1099 : congel , & ! basal ice growth (m/step-->cm/day) ! LCOV_EXCL_LINE
1100 : snoice , & ! snow-ice formation (m/step-->cm/day) ! LCOV_EXCL_LINE
1101 : dsnow , & ! snow formation (m/step-->cm/day) ! LCOV_EXCL_LINE
1102 : ! iage , & ! ice age (s) ! LCOV_EXCL_LINE
1103 : mlt_onset , & ! day of year that sfc melting begins ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1115 : evapn , & ! ice/snow mass sublimated/condensed (kg m-2 s-1) ! LCOV_EXCL_LINE
1116 : evapsn , & ! ice/snow mass sublimated/condensed over snow (kg m-2 s-1) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1138 : econ , & ! energy for condensation, < 0 (J m-2) ! LCOV_EXCL_LINE
1139 : etop_mlt , & ! energy for top melting, > 0 (J m-2) ! LCOV_EXCL_LINE
1140 : ebot_mlt , & ! energy for bottom melting, > 0 (J m-2) ! LCOV_EXCL_LINE
1141 : ebot_gro , & ! energy for bottom growth, < 0 (J m-2) ! LCOV_EXCL_LINE
1142 : emlt_atm , & ! total energy of brine, from atmosphere (J m-2) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
1147 : dhi , & ! change in ice thickness ! LCOV_EXCL_LINE
1148 : dhs , & ! change in snow thickness ! LCOV_EXCL_LINE
1149 : Ti , & ! ice temperature ! LCOV_EXCL_LINE
1150 : Ts , & ! snow temperature ! LCOV_EXCL_LINE
1151 : qbot , & ! enthalpy of ice growing at bottom surface (J m-3) ! LCOV_EXCL_LINE
1152 : qsub , & ! energy/unit volume to sublimate ice/snow (J m-3) ! LCOV_EXCL_LINE
1153 : hqtot , & ! sum of h*q for two layers ! LCOV_EXCL_LINE
1154 : wk1 , & ! temporary variable ! LCOV_EXCL_LINE
1155 : zqsnew , & ! enthalpy of new snow (J m-3) ! LCOV_EXCL_LINE
1156 : hstot , & ! snow thickness including new snow (m) ! LCOV_EXCL_LINE
1157 : Tmlts ! melting temperature (deg C)
1158 :
1159 : real (kind=dbl_kind), dimension (nilyr+1) :: &
1160 3198372490 : zi1 , & ! depth of ice layer boundaries (m) ! LCOV_EXCL_LINE
1161 3198372490 : zi2 ! adjusted depths, with equal hilyr (m)
1162 :
1163 : real (kind=dbl_kind), dimension (nslyr+1) :: &
1164 3198372490 : zs1 , & ! depth of snow layer boundaries (m) ! LCOV_EXCL_LINE
1165 3198372490 : zs2 ! adjusted depths, with equal hslyr (m)
1166 :
1167 : real (kind=dbl_kind), dimension (nilyr) :: &
1168 4797558735 : dzi ! ice layer thickness after growth/melting (m)
1169 :
1170 : real (kind=dbl_kind), dimension (nslyr) :: &
1171 3198372490 : dzs ! snow layer thickness after growth/melting (m)
1172 :
1173 : real (kind=dbl_kind), dimension (nilyr) :: &
1174 3198372490 : qm , & ! energy of melting (J m-3) = zqin in BL99 formulation ! LCOV_EXCL_LINE
1175 1599186245 : 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 ! LCOV_EXCL_LINE
1179 : qbotp , & ! enthalpy needed to grow new congelation ice (includes ocean enthalpy) ! LCOV_EXCL_LINE
1180 : qbotw , & ! enthalpy transferred to ocean during congelation freezing ! LCOV_EXCL_LINE
1181 : mass , & ! total mass from snow density tracers (kg/m^2) ! LCOV_EXCL_LINE
1182 : massi , & ! ice mass change factor ! LCOV_EXCL_LINE
1183 : tmp1 ! temporary scalar
1184 :
1185 : character(len=*),parameter :: subname='(thickness_changes)'
1186 :
1187 : !-----------------------------------------------------------------
1188 : ! Initialize
1189 : !-----------------------------------------------------------------
1190 :
1191 1599186245 : dhi = c0
1192 1599186245 : dhs = c0
1193 1599186245 : hsn_new = c0
1194 :
1195 11891141500 : do k = 1, nilyr
1196 11891141500 : dzi(k) = hilyr
1197 : enddo
1198 :
1199 3381641432 : do k = 1, nslyr
1200 3381641432 : dzs(k) = hslyr
1201 : enddo
1202 :
1203 11891141500 : do k = 1, nilyr
1204 10291955255 : if (ktherm == 2) then
1205 9856197512 : qmlt(k) = enthalpy_of_melting(zSin(k))
1206 : else
1207 435757743 : qmlt(k) = c0
1208 : endif
1209 10291955255 : qm(k) = zqin(k) - qmlt(k)
1210 10291955255 : emlt_atm = c0
1211 11891141500 : 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 1599186245 : 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 1599186245 : wk1 = -flatn * dt
1256 1599186245 : esub = max(wk1, c0) ! energy for sublimation, > 0
1257 1599186245 : econ = min(wk1, c0) ! energy for condensation, < 0
1258 :
1259 1599186245 : wk1 = (fsurfn - fcondtopn) * dt
1260 1599186245 : etop_mlt = max(wk1, c0) ! etop_mlt > 0
1261 :
1262 1599186245 : wk1 = (fcondbotn - fbot) * dt
1263 1599186245 : ebot_mlt = max(wk1, c0) ! ebot_mlt > 0
1264 1599186245 : 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 1599186245 : evapn = c0 ! initialize
1274 1599186245 : evapsn = c0 ! initialize
1275 1599186245 : evapin = c0 ! initialize
1276 :
1277 1599186245 : if (hsn > puny) then ! add snow with enthalpy zqsn(1)
1278 1383520165 : dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0
1279 :
1280 : ! assume all condensation becomes ice (no liquid)
1281 1383520165 : massice(1) = massice(1) + dhs*rhos
1282 :
1283 1383520165 : dzs(1) = dzs(1) + dhs
1284 1383520165 : evapn = evapn + dhs*rhos
1285 1383520165 : evapsn = evapsn + dhs*rhos
1286 : else ! add ice with enthalpy zqin(1)
1287 215666080 : dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0
1288 215666080 : dzi(1) = dzi(1) + dhi
1289 215666080 : evapn = evapn + dhi*rhoi
1290 215666080 : evapin = evapin + dhi*rhoi
1291 : ! enthalpy of melt water
1292 215666080 : emlt_atm = emlt_atm - qmlt(1) * dhi
1293 : endif
1294 :
1295 : !--------------------------------------------------------------
1296 : ! Grow ice (bottom)
1297 : !--------------------------------------------------------------
1298 :
1299 1599186245 : if (ktherm == 2) then
1300 :
1301 1408028216 : if (congel_freeze == 'one-step') then
1302 : ! Plante et al., The Cryosphere, 18, 1685-1708, 2024
1303 205470 : qbotm = enthalpy_mush_liquid_fraction(Tbot, phi_i_mushy)
1304 205470 : qbotw = enthalpy_brine(sst)
1305 205470 : qbotp = qbotm - qbotw
1306 205470 : dhi = ebot_gro / qbotp ! dhi > 0
1307 205470 : hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss*phi_i_mushy
1308 205470 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
1309 205470 : emlt_ocn = emlt_ocn - qbotw * dhi
1310 : else ! two-step
1311 1407822746 : qbotm = icepack_enthalpy_mush(Tbot, sss)
1312 1407822746 : qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
1313 1407822746 : qbotw = qbotm - qbotp
1314 1407822746 : dhi = ebot_gro / qbotp ! dhi > 0
1315 1407822746 : hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
1316 1407822746 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
1317 1407822746 : emlt_ocn = emlt_ocn - qbotw * dhi
1318 : endif
1319 :
1320 : else
1321 :
1322 191158029 : Tmlts = -zSin(nilyr) * depressT
1323 :
1324 : ! enthalpy of new ice growing at bottom surface
1325 191158029 : if (l_brine) then
1326 191158029 : qbotmax = -p5*rhoi*Lfresh ! max enthalpy of ice growing at bottom
1327 : qbot = -rhoi * (cp_ice * (Tmlts-Tbot) &
1328 : + Lfresh * (c1-Tmlts/Tbot) & ! LCOV_EXCL_LINE
1329 191158029 : - cp_ocn * Tmlts)
1330 191158029 : qbot = min (qbot, qbotmax) ! in case Tbot is close to Tmlt
1331 : else
1332 0 : qbot = -rhoi * (-cp_ice * Tbot + Lfresh)
1333 : endif
1334 191158029 : dhi = ebot_gro / qbot ! dhi > 0
1335 :
1336 191158029 : hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbot
1337 191158029 : hstot = c0
1338 : endif ! ktherm
1339 :
1340 1599186245 : dzi(nilyr) = dzi(nilyr) + dhi
1341 1599186245 : if (dzi(nilyr) > puny) then
1342 1599186245 : zqin(nilyr) = hqtot / dzi(nilyr)
1343 1599186245 : if (ktherm == 2) then
1344 1408028216 : zSin(nilyr) = hstot / dzi(nilyr)
1345 1408028216 : qmlt(nilyr) = enthalpy_of_melting(zSin(nilyr))
1346 : else
1347 191158029 : qmlt(nilyr) = c0
1348 : endif
1349 1599186245 : 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 1599186245 : congel = congel + dhi
1358 1599186245 : if (dhi > puny .and. frz_onset < puny) &
1359 690824 : frz_onset = yday
1360 :
1361 3381641432 : 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 1782455187 : if (ktherm == 2 .and. zqsn(k) > -rhos * Lfresh) then
1373 : dhs = max(-dzs(k), &
1374 114388255 : -((zqsn(k) + rhos*Lfresh) / (rhos*Lfresh)) * dzs(k))
1375 :
1376 114388255 : mass = massice(k) + massliq(k)
1377 114388255 : massi = c0
1378 114388255 : if (dzs(k) > puny) massi = max(c0, c1 + dhs/dzs(k))
1379 114388255 : massice(k) = massice(k) * massi
1380 114388255 : massliq(k) = mass - massice(k) ! conserve mass
1381 :
1382 114388255 : dzs(k) = dzs(k) + dhs ! dhs < 0
1383 114388255 : zqsn(k) = -rhos * Lfresh
1384 114388255 : melts = melts - dhs
1385 : ! delta E = zqsn(k) + rhos * Lfresh
1386 : endif
1387 :
1388 : !--------------------------------------------------------------
1389 : ! Sublimation of snow (evapn < 0)
1390 : !--------------------------------------------------------------
1391 :
1392 1782455187 : qsub = zqsn(k) - rhos*Lvap ! qsub < 0
1393 1782455187 : dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0
1394 :
1395 1782455187 : mass = massice(k) + massliq(k)
1396 1782455187 : massi = c0
1397 1782455187 : if (dzs(k) > puny) massi = c1 + dhs/dzs(k)
1398 1782455187 : massice(k) = massice(k) * massi
1399 1782455187 : massliq(k) = max(c0, mass + rhos*dhs - massice(k)) ! conserve new total mass
1400 :
1401 1782455187 : dzs(k) = dzs(k) + dhs
1402 1782455187 : esub = esub - dhs*qsub
1403 1782455187 : esub = max(esub, c0) ! in case of roundoff error
1404 1782455187 : evapn = evapn + dhs*rhos
1405 1782455187 : evapsn = evapsn + dhs*rhos
1406 :
1407 : !--------------------------------------------------------------
1408 : ! Melt snow (top)
1409 : !--------------------------------------------------------------
1410 :
1411 1782455187 : dhs = max(-dzs(k), etop_mlt/zqsn(k))
1412 :
1413 1782455187 : mass = massice(k) + massliq(k)
1414 1782455187 : massi = c0
1415 1782455187 : if (dzs(k) > puny) massi = max(c0, c1 + dhs/dzs(k))
1416 1782455187 : massice(k) = massice(k) * massi
1417 1782455187 : massliq(k) = mass - massice(k) ! conserve mass
1418 :
1419 1782455187 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1420 1782455187 : etop_mlt = etop_mlt - dhs*zqsn(k)
1421 1782455187 : etop_mlt = max(etop_mlt, c0) ! in case of roundoff error
1422 :
1423 : ! history diagnostics
1424 1782455187 : if (dhs < -puny .and. mlt_onset < puny) &
1425 499283 : mlt_onset = yday
1426 3381641432 : melts = melts - dhs
1427 :
1428 : enddo ! nslyr
1429 :
1430 11891141500 : do k = 1, nilyr
1431 :
1432 : !--------------------------------------------------------------
1433 : ! Sublimation of ice (evapn < 0)
1434 : !--------------------------------------------------------------
1435 :
1436 10291955255 : qsub = qm(k) - rhoi*Lvap ! qsub < 0
1437 10291955255 : dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0
1438 10291955255 : dzi(k) = dzi(k) + dhi
1439 10291955255 : esub = esub - dhi*qsub
1440 10291955255 : esub = max(esub, c0)
1441 10291955255 : evapn = evapn + dhi*rhoi
1442 10291955255 : evapin = evapin + dhi*rhoi
1443 10291955255 : emlt_ocn = emlt_ocn - qmlt(k) * dhi
1444 :
1445 : !--------------------------------------------------------------
1446 : ! Melt ice (top)
1447 : !--------------------------------------------------------------
1448 :
1449 10291955255 : if (qm(k) < c0) then
1450 10291860748 : dhi = max(-dzi(k), etop_mlt/qm(k))
1451 : else
1452 94507 : qm(k) = c0
1453 94507 : dhi = -dzi(k)
1454 : endif
1455 10291955255 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1456 :
1457 10291955255 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1458 10291955255 : etop_mlt = max(etop_mlt - dhi*qm(k), c0)
1459 :
1460 : ! history diagnostics
1461 10291955255 : if (dhi < -puny .and. mlt_onset < puny) &
1462 107439 : mlt_onset = yday
1463 11891141500 : meltt = meltt - dhi
1464 :
1465 : enddo ! nilyr
1466 :
1467 11891141500 : do k = nilyr, 1, -1
1468 :
1469 : !--------------------------------------------------------------
1470 : ! Melt ice (bottom)
1471 : !--------------------------------------------------------------
1472 :
1473 10291955255 : if (qm(k) < c0) then
1474 10291860748 : dhi = max(-dzi(k), ebot_mlt/qm(k))
1475 : else
1476 94507 : qm(k) = c0
1477 94507 : dhi = -dzi(k)
1478 : endif
1479 10291955255 : emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
1480 :
1481 10291955255 : dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0
1482 10291955255 : ebot_mlt = max(ebot_mlt - dhi*qm(k), c0)
1483 :
1484 : ! history diagnostics
1485 11891141500 : meltb = meltb -dhi
1486 :
1487 : enddo ! nilyr
1488 :
1489 3381641432 : do k = nslyr, 1, -1
1490 :
1491 : !--------------------------------------------------------------
1492 : ! Melt snow (only if all the ice has melted)
1493 : !--------------------------------------------------------------
1494 :
1495 1782455187 : dhs = max(-dzs(k), ebot_mlt/zqsn(k))
1496 :
1497 1782455187 : mass = massice(k) + massliq(k)
1498 1782455187 : massi = c0
1499 1782455187 : if (dzs(k) > puny) massi = max(c0, c1 + dhs/dzs(k))
1500 1782455187 : massice(k) = massice(k) * massi
1501 1782455187 : massliq(k) = mass - massice(k) ! conserve mass
1502 :
1503 1782455187 : dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0
1504 1782455187 : ebot_mlt = ebot_mlt - dhs*zqsn(k)
1505 1782455187 : ebot_mlt = max(ebot_mlt, c0)
1506 3381641432 : 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 1599186245 : + (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 1599186245 : if (fsnow > c0) then
1526 885502418 : hsn_new = fsnow/rhos * dt
1527 885502418 : zqsnew = -rhos*Lfresh
1528 885502418 : hstot = dzs(1) + hsn_new
1529 :
1530 885502418 : if (hstot > c0) then
1531 : zqsn(1) = (dzs(1) * zqsn(1) &
1532 885502418 : + hsn_new * zqsnew) / hstot
1533 885502418 : zqsn(1) = min(zqsn(1), zqsnew) ! avoid roundoff errors
1534 885502418 : dzs(1) = hstot
1535 885502418 : 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 1599186245 : massliq(1) = massliq(1) + frain*dt
1544 :
1545 : !-----------------------------------------------------------------
1546 : ! Find the new ice and snow thicknesses.
1547 : !-----------------------------------------------------------------
1548 :
1549 1599186245 : hin = c0
1550 1599186245 : hsn = c0
1551 :
1552 11891141500 : do k = 1, nilyr
1553 11891141500 : hin = hin + dzi(k)
1554 : enddo ! k
1555 :
1556 3381641432 : do k = 1, nslyr
1557 1782455187 : hsn = hsn + dzs(k)
1558 3381641432 : dsnow = dsnow + dzs(k) - hslyr
1559 : enddo ! k
1560 :
1561 : !-------------------------------------------------------------------
1562 : ! Incorporate new snow for snow grain radius in upper layer
1563 : !-------------------------------------------------------------------
1564 :
1565 1599186245 : if (snwgrain .and. hsn_new > c0) then
1566 24142048 : tmp1 = max(c0, dzs(1) - hsn_new)
1567 : rsnw(1) = (rsnw_fall * hsn_new + rsnw(1) * tmp1) &
1568 24142048 : / max( hsn_new + tmp1, puny)
1569 24142048 : 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 1599186245 : if (ktherm /= 2) &
1577 : call freeboard (snoice, & ! LCOV_EXCL_LINE
1578 : hin, hsn, & ! LCOV_EXCL_LINE
1579 : zqin, zqsn, & ! LCOV_EXCL_LINE
1580 : dzi, dzs, & ! LCOV_EXCL_LINE
1581 : dsnow, & ! LCOV_EXCL_LINE
1582 191158029 : massice, massliq)
1583 1599186245 : if (icepack_warnings_aborted(subname)) return
1584 :
1585 : !-------------------------------------------------------------------
1586 : ! Update snow mass tracers for uneven layers
1587 : !-------------------------------------------------------------------
1588 :
1589 1599186245 : if (snwgrain) then
1590 :
1591 203601768 : do k = 1, nslyr
1592 169668140 : meltsliq = meltsliq + massliq(k) ! used in drain_snow when all snow has melted
1593 203601768 : if (dzs(k) > puny) then
1594 135402094 : smice(k) = massice(k) / dzs(k)
1595 135402094 : smliq(k) = massliq(k) / dzs(k)
1596 : else
1597 34266046 : smice(k) = c0 ! reset to rhos below
1598 34266046 : smliq(k) = c0
1599 34266046 : massice(k) = c0
1600 34266046 : massliq(k) = c0
1601 : endif
1602 : enddo
1603 :
1604 : !-------------------------------------------------------------------
1605 : ! Check for negative snow mass tracers
1606 : !-------------------------------------------------------------------
1607 :
1608 203601768 : do k = 1, nslyr
1609 169668140 : 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 169668140 : 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 203601768 : 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 1599186245 : if (hin > c0) then
1655 1599064605 : hilyr = hin / real(nilyr,kind=dbl_kind)
1656 : else
1657 121640 : hin = c0
1658 121640 : hilyr = c0
1659 : endif
1660 1599186245 : if (hsn > c0) then
1661 1354571852 : hslyr = hsn / real(nslyr,kind=dbl_kind)
1662 : else
1663 244614393 : hsn = c0
1664 244614393 : 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 1599186245 : zi1(1) = c0
1673 1599186245 : zi1(1+nilyr) = hin
1674 :
1675 1599186245 : zi2(1) = c0
1676 1599186245 : zi2(1+nilyr) = hin
1677 :
1678 10291955255 : do k = 1, nilyr-1
1679 8692769010 : zi1(k+1) = zi1(k) + dzi(k)
1680 10141563845 : 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, & ! LCOV_EXCL_LINE
1689 : hilyr, hin, & ! LCOV_EXCL_LINE
1690 1599186245 : zqin)
1691 1599186245 : if (icepack_warnings_aborted(subname)) return
1692 :
1693 1599186245 : if (ktherm == 2) &
1694 : call adjust_enthalpy (nilyr, & ! LCOV_EXCL_LINE
1695 : zi1, zi2, & ! LCOV_EXCL_LINE
1696 : hilyr, hin, & ! LCOV_EXCL_LINE
1697 1408028216 : zSin)
1698 1599186245 : if (icepack_warnings_aborted(subname)) return
1699 :
1700 1599186245 : 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 57502466 : zs1(1) = c0
1708 57502466 : zs1(1+nslyr) = hsn
1709 :
1710 57502466 : zs2(1) = c0
1711 57502466 : zs2(1+nslyr) = hsn
1712 :
1713 240771408 : do k = 1, nslyr-1
1714 183268942 : zs1(k+1) = zs1(k) + dzs(k)
1715 240771408 : 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, & ! LCOV_EXCL_LINE
1725 : hslyr, hsn, & ! LCOV_EXCL_LINE
1726 57502466 : zqsn)
1727 :
1728 57502466 : if (snwgrain) then
1729 : call adjust_enthalpy (nslyr, &
1730 : zs1(:), zs2(:), & ! LCOV_EXCL_LINE
1731 : hslyr, hsn, & ! LCOV_EXCL_LINE
1732 33933628 : rsnw(:))
1733 : call adjust_enthalpy (nslyr, & ! need a routine to adjust
1734 : zs1(:), zs2(:), & ! mass instead of tracer ! LCOV_EXCL_LINE
1735 : hslyr, hsn, & ! LCOV_EXCL_LINE
1736 33933628 : smice(:))
1737 : call adjust_enthalpy (nslyr, &
1738 : zs1(:), zs2(:), & ! LCOV_EXCL_LINE
1739 : hslyr, hsn, & ! LCOV_EXCL_LINE
1740 33933628 : smliq(:))
1741 : ! Update snow mass
1742 203601768 : do k = 1, nslyr
1743 169668140 : massice(k) = smice(k) * hslyr
1744 203601768 : massliq(k) = smliq(k) * hslyr
1745 : enddo
1746 : endif
1747 57502466 : if (icepack_warnings_aborted(subname)) return
1748 :
1749 : endif ! nslyr > 1
1750 :
1751 : !-----------------------------------------------------------------
1752 : ! Remove very thin snow layers (ktherm = 2)
1753 : !-----------------------------------------------------------------
1754 :
1755 1599186245 : if (ktherm == 2) then
1756 1408028216 : if (hsn <= puny .or. hin <= c0) then
1757 349104920 : do k = 1, nslyr
1758 : fhocnn = fhocnn &
1759 185500007 : + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
1760 185500007 : zqsn(k) = -rhos*Lfresh
1761 349104920 : if (snwgrain) then
1762 20934140 : meltsliq = meltsliq + massice(k) ! add to meltponds
1763 20934140 : smice(k) = rhos
1764 20934140 : smliq(k) = c0
1765 20934140 : rsnw(k) = rsnw_fall
1766 : endif
1767 : enddo
1768 163604913 : melts = melts + hsn
1769 163604913 : hsn = c0
1770 163604913 : 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 1599186245 : efinal = -evapn*Lvap
1780 1599186245 : evapn = evapn/dt
1781 1599186245 : evapsn = evapsn/dt
1782 1599186245 : evapin = evapin/dt
1783 :
1784 3381641432 : do k = 1, nslyr
1785 3381641432 : efinal = efinal + hslyr*zqsn(k)
1786 : enddo
1787 :
1788 11891141500 : do k = 1, nilyr
1789 11891141500 : efinal = efinal + hilyr*zqin(k)
1790 : enddo ! k
1791 :
1792 1599186245 : if (ktherm < 2) then
1793 191158029 : emlt_atm = c0
1794 191158029 : emlt_ocn = c0
1795 : endif
1796 :
1797 : ! melt water is no longer zero enthalpy with ktherm=2
1798 1599186245 : fhocnn = fhocnn + emlt_ocn/dt
1799 1599186245 : 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 191158029 : subroutine freeboard (snoice, &
1813 : hin, hsn, & ! LCOV_EXCL_LINE
1814 0 : zqin, zqsn, & ! LCOV_EXCL_LINE
1815 382316058 : dzi, dzs, & ! LCOV_EXCL_LINE
1816 : dsnow, & ! LCOV_EXCL_LINE
1817 191158029 : 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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1836 : dzi , & ! ice layer thicknesses (m) ! LCOV_EXCL_LINE
1837 : dzs , & ! snow layer thicknesses (m) ! LCOV_EXCL_LINE
1838 : massice , & ! total ice mass of snow in each layer (kg/m^2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1848 : dhsn , & ! change in snow thickness (m) ! LCOV_EXCL_LINE
1849 : hqs ! sum of h*q for snow (J m-2)
1850 :
1851 : real (kind=dbl_kind) :: &
1852 : wk1 , & ! temporary variable ! LCOV_EXCL_LINE
1853 : dhs , & ! snow to remove from layer (m) ! LCOV_EXCL_LINE
1854 : mass , & ! total snow mass from tracers (kg/m^2) ! LCOV_EXCL_LINE
1855 : massi ! mass change factor
1856 :
1857 : character(len=*),parameter :: subname='(freeboard)'
1858 :
1859 : !-----------------------------------------------------------------
1860 : ! Determine whether snow lies below freeboard.
1861 : !-----------------------------------------------------------------
1862 :
1863 191158029 : dhin = c0
1864 191158029 : dhsn = c0
1865 191158029 : hqs = c0
1866 :
1867 191158029 : wk1 = hsn - hin*(rhow-rhoi)/rhos ! not yet consistent with smice/smliq
1868 :
1869 191158029 : if (wk1 > puny .and. hsn > puny) then ! snow below freeboard
1870 1346063 : dhsn = min(wk1*rhoi/rhow, hsn) ! snow to remove
1871 1346063 : 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 382316058 : do k = nslyr, 1, -1
1880 382316058 : if (dhin > puny) then
1881 1346063 : dhs = min(dhsn, dzs(k)) ! snow to remove from layer
1882 1346063 : hsn = hsn - dhs
1883 1346063 : dsnow = dsnow - dhs ! new snow
1884 :
1885 : ! remove both ice and liquid from snow to add to ice
1886 1346063 : mass = massice(k) + massliq(k)
1887 1346063 : massi = c0
1888 1346063 : if (dzs(k) > puny) massi = max(c0, c1 - dhs/dzs(k))
1889 1346063 : massice(k) = massice(k) * massi
1890 1346063 : 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 1346063 : dzs(k) = dzs(k) - dhs
1895 1346063 : dhsn = dhsn - dhs
1896 1346063 : dhsn = max(dhsn,c0)
1897 1346063 : 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 191158029 : 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 1346063 : wk1 = dzi(1) + dhin
1911 1346063 : hin = hin + dhin
1912 1346063 : zqin(1) = (dzi(1)*zqin(1) + hqs) / wk1
1913 1346063 : dzi(1) = wk1
1914 :
1915 : ! history diagnostic
1916 1346063 : snoice = snoice + dhin
1917 : endif ! dhin > puny
1918 :
1919 191158029 : 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 1599186245 : subroutine conservation_check_vthermo(dt, &
1931 : fsurfn, flatn, & ! LCOV_EXCL_LINE
1932 : fhocnn, fswint, & ! LCOV_EXCL_LINE
1933 : fsnow, & ! LCOV_EXCL_LINE
1934 : einit, einter, & ! LCOV_EXCL_LINE
1935 : efinal, & ! LCOV_EXCL_LINE
1936 : fcondtopn,fcondbotn, & ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
1944 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
1945 : fhocnn , & ! fbot, corrected for any surplus energy ! LCOV_EXCL_LINE
1946 : fswint , & ! SW absorbed in ice interior, below surface (W m-2) ! LCOV_EXCL_LINE
1947 : fsnow , & ! snowfall rate (kg m-2 s-1) ! LCOV_EXCL_LINE
1948 : fcondtopn , & ! LCOV_EXCL_LINE
1949 : fadvocn , & ! LCOV_EXCL_LINE
1950 : fbot
1951 :
1952 : real (kind=dbl_kind), intent(in) :: &
1953 : einit , & ! initial energy of melting (J m-2) ! LCOV_EXCL_LINE
1954 : einter , & ! intermediate energy of melting (J m-2) ! LCOV_EXCL_LINE
1955 : efinal , & ! final energy of melting (J m-2) ! LCOV_EXCL_LINE
1956 : fcondbotn
1957 :
1958 : ! local variables
1959 :
1960 : real (kind=dbl_kind) :: &
1961 : einp , & ! energy input during timestep (J m-2) ! LCOV_EXCL_LINE
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 1599186245 : - fsnow*Lfresh - fadvocn) * dt
1979 1599186245 : ferr = abs(efinal-einit-einp) / dt
1980 :
1981 1599186245 : 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, & ! LCOV_EXCL_LINE
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 1599186245 : subroutine update_state_vthermo(Tf, Tsf, &
2052 : hin, hsn, & ! LCOV_EXCL_LINE
2053 1599186245 : zqin, zSin, & ! LCOV_EXCL_LINE
2054 1599186245 : zqsn, & ! LCOV_EXCL_LINE
2055 : aicen, vicen, & ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2066 : hsn ! snow thickness (m)
2067 :
2068 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
2069 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
2070 : zSin , & ! ice salinity (ppt) ! LCOV_EXCL_LINE
2071 : zqsn ! snow layer enthalpy (J m-3)
2072 :
2073 : real (kind=dbl_kind), intent(inout) :: &
2074 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
2075 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
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 1599186245 : if (hin <= c0) then
2086 121640 : aicen = c0
2087 121640 : vicen = c0
2088 121640 : vsnon = c0
2089 121640 : Tsf = Tf
2090 973120 : do k = 1, nilyr
2091 973120 : zqin(k) = c0
2092 : enddo
2093 121640 : if (ktherm == 2) then
2094 973120 : do k = 1, nilyr
2095 973120 : zSin(k) = c0
2096 : enddo
2097 : endif
2098 312712 : do k = 1, nslyr
2099 312712 : zqsn(k) = c0
2100 : enddo
2101 : else
2102 : ! aicen is already up to date
2103 1599064605 : vicen = aicen * hin
2104 1599064605 : vsnon = aicen * hsn
2105 : endif
2106 :
2107 1599186245 : 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 1492190040 : subroutine icepack_step_therm1(dt, &
2118 0 : aicen_init , & ! LCOV_EXCL_LINE
2119 1492190040 : vicen_init , vsnon_init , & ! LCOV_EXCL_LINE
2120 0 : aice , aicen , & ! LCOV_EXCL_LINE
2121 0 : vice , vicen , & ! LCOV_EXCL_LINE
2122 1492190040 : vsno , vsnon , & ! LCOV_EXCL_LINE
2123 : uvel , vvel , & ! LCOV_EXCL_LINE
2124 2984380080 : Tsfc , zqsn , & ! LCOV_EXCL_LINE
2125 2984380080 : zqin , zSin , & ! LCOV_EXCL_LINE
2126 2984380080 : alvl , vlvl , & ! LCOV_EXCL_LINE
2127 2984380080 : apnd , hpnd , & ! LCOV_EXCL_LINE
2128 0 : ipnd , & ! LCOV_EXCL_LINE
2129 1492190040 : iage , FY , & ! LCOV_EXCL_LINE
2130 1492190040 : aerosno , aeroice , & ! LCOV_EXCL_LINE
2131 1492190040 : isosno , isoice , & ! LCOV_EXCL_LINE
2132 : uatm , vatm , & ! LCOV_EXCL_LINE
2133 : wind , zlvl , & ! LCOV_EXCL_LINE
2134 : Qa , rhoa , & ! LCOV_EXCL_LINE
2135 1492190040 : Qa_iso , & ! LCOV_EXCL_LINE
2136 : Tair , Tref , & ! LCOV_EXCL_LINE
2137 : Qref , Uref , & ! LCOV_EXCL_LINE
2138 1492190040 : Qref_iso , & ! LCOV_EXCL_LINE
2139 : Cdn_atm_ratio, & ! LCOV_EXCL_LINE
2140 : Cdn_ocn , Cdn_ocn_skin, & ! LCOV_EXCL_LINE
2141 : Cdn_ocn_floe, Cdn_ocn_keel, & ! LCOV_EXCL_LINE
2142 : Cdn_atm , Cdn_atm_skin, & ! LCOV_EXCL_LINE
2143 : Cdn_atm_floe, Cdn_atm_pond, & ! LCOV_EXCL_LINE
2144 : Cdn_atm_rdg , hfreebd , & ! LCOV_EXCL_LINE
2145 : hdraft , hridge , & ! LCOV_EXCL_LINE
2146 : distrdg , hkeel , & ! LCOV_EXCL_LINE
2147 : dkeel , lfloe , & ! LCOV_EXCL_LINE
2148 : dfloe , & ! LCOV_EXCL_LINE
2149 : strax , stray , & ! LCOV_EXCL_LINE
2150 : strairxT , strairyT , & ! LCOV_EXCL_LINE
2151 : potT , sst , & ! LCOV_EXCL_LINE
2152 : sss , Tf , & ! LCOV_EXCL_LINE
2153 : strocnxT , strocnyT , & ! LCOV_EXCL_LINE
2154 : fbot , & ! LCOV_EXCL_LINE
2155 : Tbot , Tsnice , & ! LCOV_EXCL_LINE
2156 1492190040 : frzmlt , rsiden , & ! LCOV_EXCL_LINE
2157 : wlat , & ! LCOV_EXCL_LINE
2158 : fsnow , frain , & ! LCOV_EXCL_LINE
2159 : fpond , fsloss , & ! LCOV_EXCL_LINE
2160 0 : fsurf , fsurfn , & ! LCOV_EXCL_LINE
2161 1492190040 : fcondtop , fcondtopn , & ! LCOV_EXCL_LINE
2162 0 : fcondbot , fcondbotn , & ! LCOV_EXCL_LINE
2163 1492190040 : fswsfcn , fswintn , & ! LCOV_EXCL_LINE
2164 1492190040 : fswthrun , & ! LCOV_EXCL_LINE
2165 1492190040 : fswthrun_vdr, & ! LCOV_EXCL_LINE
2166 1492190040 : fswthrun_vdf, & ! LCOV_EXCL_LINE
2167 1492190040 : fswthrun_idr, & ! LCOV_EXCL_LINE
2168 1492190040 : fswthrun_idf, & ! LCOV_EXCL_LINE
2169 : fswabs , & ! LCOV_EXCL_LINE
2170 : flwout , & ! LCOV_EXCL_LINE
2171 1492190040 : Sswabsn , Iswabsn , & ! LCOV_EXCL_LINE
2172 : flw , & ! LCOV_EXCL_LINE
2173 1492190040 : fsens , fsensn , & ! LCOV_EXCL_LINE
2174 1492190040 : flat , flatn , & ! LCOV_EXCL_LINE
2175 : evap , & ! LCOV_EXCL_LINE
2176 : evaps , evapi , & ! LCOV_EXCL_LINE
2177 : fresh , fsalt , & ! LCOV_EXCL_LINE
2178 : fhocn , & ! LCOV_EXCL_LINE
2179 : fswthru , & ! LCOV_EXCL_LINE
2180 : fswthru_vdr , & ! LCOV_EXCL_LINE
2181 : fswthru_vdf , & ! LCOV_EXCL_LINE
2182 : fswthru_idr , & ! LCOV_EXCL_LINE
2183 : fswthru_idf , & ! LCOV_EXCL_LINE
2184 1492190040 : flatn_f , fsensn_f , & ! LCOV_EXCL_LINE
2185 1492190040 : fsurfn_f , fcondtopn_f , & ! LCOV_EXCL_LINE
2186 1492190040 : faero_atm , faero_ocn , & ! LCOV_EXCL_LINE
2187 1492190040 : fiso_atm , fiso_ocn , & ! LCOV_EXCL_LINE
2188 1492190040 : fiso_evap , & ! LCOV_EXCL_LINE
2189 : HDO_ocn , H2_16O_ocn , & ! LCOV_EXCL_LINE
2190 : H2_18O_ocn , & ! LCOV_EXCL_LINE
2191 1492190040 : dhsn , ffracn , & ! LCOV_EXCL_LINE
2192 0 : meltt , melttn , & ! LCOV_EXCL_LINE
2193 1492190040 : meltb , meltbn , & ! LCOV_EXCL_LINE
2194 1492190040 : melts , meltsn , & ! LCOV_EXCL_LINE
2195 1492190040 : congel , congeln , & ! LCOV_EXCL_LINE
2196 1492190040 : snoice , snoicen , & ! LCOV_EXCL_LINE
2197 1492190040 : dsnow , dsnown , & ! LCOV_EXCL_LINE
2198 1492190040 : meltsliq , meltsliqn , & ! LCOV_EXCL_LINE
2199 1492190040 : rsnwn , & ! LCOV_EXCL_LINE
2200 1492190040 : smicen , smliqn , & ! LCOV_EXCL_LINE
2201 : lmask_n , lmask_s , & ! LCOV_EXCL_LINE
2202 : mlt_onset , frz_onset , & ! LCOV_EXCL_LINE
2203 : yday , prescribed_ice, & ! LCOV_EXCL_LINE
2204 1492190040 : zlvs , afsdn)
2205 :
2206 : real (kind=dbl_kind), intent(in) :: &
2207 : dt , & ! time step ! LCOV_EXCL_LINE
2208 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
2209 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
2210 : strax , & ! wind stress components (N/m^2) ! LCOV_EXCL_LINE
2211 : stray , & ! ! LCOV_EXCL_LINE
2212 : yday ! day of year
2213 :
2214 : logical (kind=log_kind), intent(in) :: &
2215 : lmask_n , & ! northern hemisphere mask ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
2223 : vice , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
2224 : vsno , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
2225 : zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) ! LCOV_EXCL_LINE
2226 : uatm , & ! wind velocity components (m/s) ! LCOV_EXCL_LINE
2227 : vatm , & ! (m/s) ! LCOV_EXCL_LINE
2228 : wind , & ! wind speed (m/s) ! LCOV_EXCL_LINE
2229 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
2230 : Tair , & ! air temperature (K) ! LCOV_EXCL_LINE
2231 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
2232 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
2233 : frain , & ! rainfall rate (kg/m^2 s) ! LCOV_EXCL_LINE
2234 : fsnow , & ! snowfall rate (kg/m^2 s) ! LCOV_EXCL_LINE
2235 : fpond , & ! fresh water flux to ponds (kg/m^2/s) ! LCOV_EXCL_LINE
2236 : fresh , & ! fresh water flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
2237 : fsalt , & ! salt flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
2238 : fhocn , & ! net heat flux to ocean (W/m^2) ! LCOV_EXCL_LINE
2239 : fswthru , & ! shortwave penetrating to ocean (W/m^2) ! LCOV_EXCL_LINE
2240 : fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2) ! LCOV_EXCL_LINE
2241 : fcondtop , & ! top surface conductive flux (W/m^2) ! LCOV_EXCL_LINE
2242 : fcondbot , & ! bottom surface conductive flux (W/m^2) ! LCOV_EXCL_LINE
2243 : fsens , & ! sensible heat flux (W/m^2) ! LCOV_EXCL_LINE
2244 : flat , & ! latent heat flux (W/m^2) ! LCOV_EXCL_LINE
2245 : fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2) ! LCOV_EXCL_LINE
2246 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
2247 : flwout , & ! outgoing longwave radiation (W/m^2) ! LCOV_EXCL_LINE
2248 : evap , & ! evaporative water flux (kg/m^2/s) ! LCOV_EXCL_LINE
2249 : evaps , & ! evaporative water flux over snow(kg/m^2/s) ! LCOV_EXCL_LINE
2250 : evapi , & ! evaporative water flux over ice (kg/m^2/s) ! LCOV_EXCL_LINE
2251 : congel , & ! basal ice growth (m/step-->cm/day) ! LCOV_EXCL_LINE
2252 : snoice , & ! snow-ice formation (m/step-->cm/day) ! LCOV_EXCL_LINE
2253 : Tref , & ! 2m atm reference temperature (K) ! LCOV_EXCL_LINE
2254 : Qref , & ! 2m atm reference spec humidity (kg/kg) ! LCOV_EXCL_LINE
2255 : Uref , & ! 10m atm reference wind speed (m/s) ! LCOV_EXCL_LINE
2256 : Cdn_atm , & ! atm drag coefficient ! LCOV_EXCL_LINE
2257 : Cdn_ocn , & ! ocn drag coefficient ! LCOV_EXCL_LINE
2258 : hfreebd , & ! freeboard (m) ! LCOV_EXCL_LINE
2259 : hdraft , & ! draft of ice + snow column (Stoessel1993) ! LCOV_EXCL_LINE
2260 : hridge , & ! ridge height ! LCOV_EXCL_LINE
2261 : distrdg , & ! distance between ridges ! LCOV_EXCL_LINE
2262 : hkeel , & ! keel depth ! LCOV_EXCL_LINE
2263 : dkeel , & ! distance between keels ! LCOV_EXCL_LINE
2264 : lfloe , & ! floe length ! LCOV_EXCL_LINE
2265 : dfloe , & ! distance between floes ! LCOV_EXCL_LINE
2266 : Cdn_atm_skin, & ! neutral skin drag coefficient ! LCOV_EXCL_LINE
2267 : Cdn_atm_floe, & ! neutral floe edge drag coefficient ! LCOV_EXCL_LINE
2268 : Cdn_atm_pond, & ! neutral pond edge drag coefficient ! LCOV_EXCL_LINE
2269 : Cdn_atm_rdg , & ! neutral ridge drag coefficient ! LCOV_EXCL_LINE
2270 : Cdn_ocn_skin, & ! skin drag coefficient ! LCOV_EXCL_LINE
2271 : Cdn_ocn_floe, & ! floe edge drag coefficient ! LCOV_EXCL_LINE
2272 : Cdn_ocn_keel, & ! keel drag coefficient ! LCOV_EXCL_LINE
2273 : Cdn_atm_ratio,& ! ratio drag atm / neutral drag atm ! LCOV_EXCL_LINE
2274 : strairxT , & ! stress on ice by air, x-direction ! LCOV_EXCL_LINE
2275 : strairyT , & ! stress on ice by air, y-direction ! LCOV_EXCL_LINE
2276 : strocnxT , & ! ice-ocean stress, x-direction ! LCOV_EXCL_LINE
2277 : strocnyT , & ! ice-ocean stress, y-direction ! LCOV_EXCL_LINE
2278 : fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) ! LCOV_EXCL_LINE
2279 : frzmlt , & ! freezing/melting potential (W/m^2) ! LCOV_EXCL_LINE
2280 : sst , & ! sea surface temperature (C) ! LCOV_EXCL_LINE
2281 : Tf , & ! freezing temperature (C) ! LCOV_EXCL_LINE
2282 : Tbot , & ! ice bottom surface temperature (deg C) ! LCOV_EXCL_LINE
2283 : Tsnice , & ! snow ice interface temperature (deg C) ! LCOV_EXCL_LINE
2284 : sss , & ! sea surface salinity (ppt) ! LCOV_EXCL_LINE
2285 : meltt , & ! top ice melt (m/step-->cm/day) ! LCOV_EXCL_LINE
2286 : melts , & ! snow melt (m/step-->cm/day) ! LCOV_EXCL_LINE
2287 : meltb , & ! basal ice melt (m/step-->cm/day) ! LCOV_EXCL_LINE
2288 : mlt_onset , & ! day of year that sfc melting begins ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2296 : fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2) ! LCOV_EXCL_LINE
2297 : fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2) ! LCOV_EXCL_LINE
2298 : fswthru_idf , & ! nir dif shortwave penetrating to ocean (W/m^2) ! LCOV_EXCL_LINE
2299 : dsnow , & ! change in snow depth (m/step-->cm/day) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2307 : Qref_iso , & ! isotope 2m atm ref spec humidity (kg/kg) ! LCOV_EXCL_LINE
2308 : fiso_atm , & ! isotope deposition rate (kg/m^2 s) ! LCOV_EXCL_LINE
2309 : fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2317 : smicen , & ! tracer for mass of ice in snow (kg/m^3) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2322 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) ! LCOV_EXCL_LINE
2323 : H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
2331 : vicen_init , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
2332 : vsnon_init , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
2333 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
2334 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
2335 : vsnon , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
2336 : Tsfc , & ! ice/snow surface temperature, Tsfcn ! LCOV_EXCL_LINE
2337 : alvl , & ! level ice area fraction ! LCOV_EXCL_LINE
2338 : vlvl , & ! level ice volume fraction ! LCOV_EXCL_LINE
2339 : apnd , & ! melt pond area fraction tracer ! LCOV_EXCL_LINE
2340 : hpnd , & ! melt pond depth (m) ! LCOV_EXCL_LINE
2341 : ipnd , & ! melt pond refrozen lid thickness (m) ! LCOV_EXCL_LINE
2342 : iage , & ! volume-weighted ice age ! LCOV_EXCL_LINE
2343 : FY , & ! area-weighted first-year ice area ! LCOV_EXCL_LINE
2344 : rsiden , & ! fraction of ice that melts laterally ! LCOV_EXCL_LINE
2345 : fsurfn , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
2346 : fcondtopn , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
2347 : fcondbotn , & ! downward cond flux at bottom surface (W m-2) ! LCOV_EXCL_LINE
2348 : flatn , & ! latent heat flux (W m-2) ! LCOV_EXCL_LINE
2349 : fsensn , & ! sensible heat flux (W m-2) ! LCOV_EXCL_LINE
2350 : fsurfn_f , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
2351 : fcondtopn_f , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
2352 : flatn_f , & ! latent heat flux (W m-2) ! LCOV_EXCL_LINE
2353 : fsensn_f , & ! sensible heat flux (W m-2) ! LCOV_EXCL_LINE
2354 : fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
2355 : fswintn , & ! SW absorbed in ice interior, below surface (W m-2) ! LCOV_EXCL_LINE
2356 : faero_atm , & ! aerosol deposition rate (kg/m^2 s) ! LCOV_EXCL_LINE
2357 : faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
2358 : dhsn , & ! depth difference for snow on sea ice and pond ice ! LCOV_EXCL_LINE
2359 : ffracn , & ! fraction of fsurfn used to melt ipond ! LCOV_EXCL_LINE
2360 : meltsn , & ! snow melt (m) ! LCOV_EXCL_LINE
2361 : melttn , & ! top ice melt (m) ! LCOV_EXCL_LINE
2362 : meltbn , & ! bottom ice melt (m) ! LCOV_EXCL_LINE
2363 : congeln , & ! congelation ice growth (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2374 : fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) ! LCOV_EXCL_LINE
2375 : fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2380 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
2381 : zSin , & ! internal ice layer salinities ! LCOV_EXCL_LINE
2382 : Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
2398 : n ! category index
2399 :
2400 : real (kind=dbl_kind) :: &
2401 : worka , & ! temporary variables ! LCOV_EXCL_LINE
2402 : workb , & ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2408 : flwoutn , & ! upward LW at surface (W/m^2) ! LCOV_EXCL_LINE
2409 : evapn , & ! flux of vapor, atmos to ice (kg m-2 s-1) ! LCOV_EXCL_LINE
2410 : evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1) ! LCOV_EXCL_LINE
2411 : evapin , & ! flux of vapor, atmos to ice over ice (kg m-2 s-1) ! LCOV_EXCL_LINE
2412 : freshn , & ! flux of water, ice to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
2413 : fsaltn , & ! flux of salt, ice to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
2414 : fhocnn , & ! fbot corrected for leftover energy (W/m^2) ! LCOV_EXCL_LINE
2415 : strairxn , & ! air/ice zonal stress, (N/m^2) ! LCOV_EXCL_LINE
2416 : strairyn , & ! air/ice meridional stress, (N/m^2) ! LCOV_EXCL_LINE
2417 : Cdn_atm_ratio_n, & ! drag coefficient ratio ! LCOV_EXCL_LINE
2418 : Trefn , & ! air tmp reference level (K) ! LCOV_EXCL_LINE
2419 : Urefn , & ! air speed reference level (m/s) ! LCOV_EXCL_LINE
2420 : Qrefn , & ! air sp hum reference level (kg/kg) ! LCOV_EXCL_LINE
2421 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
2422 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
2423 : rfrac ! water fraction retained for melt ponds
2424 :
2425 : real (kind=dbl_kind), dimension(nslyr,ncat) :: &
2426 2984380080 : massicen , & ! mass of ice in snow (kg/m^2) ! LCOV_EXCL_LINE
2427 2984380080 : massliqn ! mass of liquid in snow (kg/m^2)
2428 :
2429 : real (kind=dbl_kind), dimension(n_iso) :: &
2430 2984380080 : Qrefn_iso , & ! isotope air sp hum reference level (kg/kg) ! LCOV_EXCL_LINE
2431 2984380080 : fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
2432 2984380080 : fiso_evapn ! isotope evaporation (kg/m^2/s)
2433 :
2434 : real (kind=dbl_kind), dimension(nslyr) :: &
2435 2984380080 : rsnw , & ! snow grain radius (10^-6 m) ! LCOV_EXCL_LINE
2436 2984380080 : smice , & ! tracer for mass of ice in snow (kg/m^3) ! LCOV_EXCL_LINE
2437 2984380080 : smliq ! tracer for mass of liquid in snow (kg/m^3)
2438 :
2439 : real (kind=dbl_kind), dimension(ncat) :: &
2440 4476570120 : apond , & ! melt pond area fraction of category ! LCOV_EXCL_LINE
2441 1492190040 : 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) ! LCOV_EXCL_LINE
2445 : l_fswthrun_vdf, & ! vis dif SW local n ice to ocean (W/m^2) ! LCOV_EXCL_LINE
2446 : l_fswthrun_idr, & ! nir dir SW local n ice to ocean (W/m^2) ! LCOV_EXCL_LINE
2447 : l_fswthrun_idf, & ! nir dif SW local n ice to ocean (W/m^2) ! LCOV_EXCL_LINE
2448 : l_dsnow, & ! local snow change ! LCOV_EXCL_LINE
2449 : l_dsnown, & ! local snow change category ! LCOV_EXCL_LINE
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 1492190040 : if (icepack_chkoptargflag(first_call)) then
2465 6403 : if (tr_iso) then
2466 207 : if (.not.(present(isosno) .and. present(isoice) .and. &
2467 : present(fiso_atm) .and. present(fiso_ocn) .and. & ! LCOV_EXCL_LINE
2468 : present(fiso_evap).and. & ! LCOV_EXCL_LINE
2469 : present(Qa_iso) .and. present(Qref_iso) .and. & ! LCOV_EXCL_LINE
2470 : present(HDO_ocn) .and. present(H2_16O_ocn) .and. & ! LCOV_EXCL_LINE
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 6403 : if (snwgrain) then
2478 429 : 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. & ! LCOV_EXCL_LINE
2487 6403 : (present(fswthru_idr) .and. .not.present(fswthrun_idr)) .or. & ! LCOV_EXCL_LINE
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 6403 : 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 6403 : if (tr_fsd) then
2500 295 : 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 295 : 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 3225553392 : rsnw (:) = c0
2518 3225553392 : smice(:) = c0
2519 3225553392 : smliq(:) = c0
2520 :
2521 1492190040 : l_meltsliq = c0
2522 8930659392 : l_meltsliqn = c0
2523 1492190040 : l_dsnow = c0
2524 1492190040 : if (present(dsnow)) l_dsnow = dsnow
2525 :
2526 : ! solid and liquid components of snow mass
2527 17633407080 : massicen(:,:) = c0
2528 17633407080 : massliqn(:,:) = c0
2529 :
2530 : !-----------------------------------------------------------------
2531 : ! Initialize pond area fractions
2532 : !-----------------------------------------------------------------
2533 8930659392 : do n= 1, ncat
2534 8930659392 : if (tr_pond_lvl) then
2535 7046304408 : apond(n) = apnd(n) * alvl(n)
2536 : else
2537 392164944 : apond(n) = apnd(n)
2538 : endif
2539 : enddo
2540 :
2541 : !-----------------------------------------------------------------
2542 : ! Initialize rate of snow loss to leads
2543 : !-----------------------------------------------------------------
2544 :
2545 1492190040 : if (present(fsloss)) &
2546 1490268600 : fsloss = fsnow * (c1 - aice)
2547 :
2548 : !-----------------------------------------------------------------
2549 : ! snow redistribution using snwlvlfac: precip factor
2550 : !-----------------------------------------------------------------
2551 :
2552 1492190040 : if (trim(snwredist) == 'bulk') then
2553 192144 : worka = c0
2554 192144 : if (aice > puny) then
2555 261546 : do n = 1, ncat
2556 261546 : worka = worka + alvl(n)*aicen(n)
2557 : enddo
2558 43591 : worka = worka * (snwlvlfac/(c1+snwlvlfac)) / aice
2559 : endif
2560 192144 : if (present(fsloss)) &
2561 192144 : fsloss = fsloss + fsnow* worka
2562 192144 : 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 1492190040 : if (formdrag) then
2571 : call neutral_drag_coeffs (apond , &
2572 : alvl , vlvl , & ! LCOV_EXCL_LINE
2573 : aice , vice, & ! LCOV_EXCL_LINE
2574 : vsno , aicen , & ! LCOV_EXCL_LINE
2575 : vicen , & ! LCOV_EXCL_LINE
2576 : Cdn_ocn , Cdn_ocn_skin, & ! LCOV_EXCL_LINE
2577 : Cdn_ocn_floe, Cdn_ocn_keel, & ! LCOV_EXCL_LINE
2578 : Cdn_atm , Cdn_atm_skin, & ! LCOV_EXCL_LINE
2579 : Cdn_atm_floe, Cdn_atm_pond, & ! LCOV_EXCL_LINE
2580 : Cdn_atm_rdg , hfreebd , & ! LCOV_EXCL_LINE
2581 : hdraft , hridge , & ! LCOV_EXCL_LINE
2582 : distrdg , hkeel , & ! LCOV_EXCL_LINE
2583 : dkeel , lfloe , & ! LCOV_EXCL_LINE
2584 14602944 : dfloe)
2585 14602944 : 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, & ! LCOV_EXCL_LINE
2596 : vicen, vsnon, & ! LCOV_EXCL_LINE
2597 : zqin, zqsn, & ! LCOV_EXCL_LINE
2598 : sst, Tf, & ! LCOV_EXCL_LINE
2599 : ustar_min, & ! LCOV_EXCL_LINE
2600 : fbot_xfer_type, & ! LCOV_EXCL_LINE
2601 : strocnxT, strocnyT, & ! LCOV_EXCL_LINE
2602 : Tbot, fbot, & ! LCOV_EXCL_LINE
2603 : rsiden, Cdn_ocn, & ! LCOV_EXCL_LINE
2604 : wlat, aicen, & ! LCOV_EXCL_LINE
2605 1492190040 : afsdn)
2606 :
2607 1492190040 : if (icepack_warnings_aborted(subname)) return
2608 :
2609 8930659392 : do n = 1, ncat
2610 7438469352 : meltsn (n) = c0
2611 7438469352 : melttn (n) = c0
2612 7438469352 : meltbn (n) = c0
2613 7438469352 : congeln(n) = c0
2614 7438469352 : snoicen(n) = c0
2615 7438469352 : l_dsnown = c0
2616 :
2617 7438469352 : Trefn = c0
2618 7438469352 : Qrefn = c0
2619 7657513512 : Qrefn_iso(:) = c0
2620 7657513512 : fiso_ocnn(:) = c0
2621 7657513512 : fiso_evapn(:) = c0
2622 7438469352 : Urefn = c0
2623 7438469352 : lhcoef = c0
2624 7438469352 : shcoef = c0
2625 7438469352 : worka = c0
2626 7438469352 : workb = c0
2627 :
2628 7438469352 : if (aicen_init(n) > puny) then
2629 :
2630 1599186245 : 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, & ! LCOV_EXCL_LINE
2643 : uatm, vatm, & ! LCOV_EXCL_LINE
2644 : wind, zlvl, & ! LCOV_EXCL_LINE
2645 : Qa, rhoa, & ! LCOV_EXCL_LINE
2646 : strairxn, strairyn, & ! LCOV_EXCL_LINE
2647 : Trefn, Qrefn, & ! LCOV_EXCL_LINE
2648 : worka, workb, & ! LCOV_EXCL_LINE
2649 : lhcoef, shcoef, & ! LCOV_EXCL_LINE
2650 : Cdn_atm, & ! LCOV_EXCL_LINE
2651 : Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
2652 : Qa_iso=Qa_iso, & ! LCOV_EXCL_LINE
2653 : Qref_iso=Qrefn_iso, & ! LCOV_EXCL_LINE
2654 : uvel=uvel, vvel=vvel, & ! LCOV_EXCL_LINE
2655 1599186245 : Uref=Urefn, zlvs=zlvs)
2656 1599186245 : if (icepack_warnings_aborted(subname)) return
2657 :
2658 : endif ! calc_Tsfc or calc_strair
2659 :
2660 1599186245 : 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 1599186245 : if (tr_iage) call increment_age (dt, iage(n))
2681 1599186245 : if (icepack_warnings_aborted(subname)) return
2682 1599186245 : if (tr_FY) call update_FYarea (dt, &
2683 : lmask_n, lmask_s, & ! LCOV_EXCL_LINE
2684 1434232947 : yday, FY(n))
2685 1599186245 : if (icepack_warnings_aborted(subname)) return
2686 :
2687 : !-----------------------------------------------------------------
2688 : ! Vertical thermodynamics: Heat conduction, growth and melting.
2689 : !-----------------------------------------------------------------
2690 :
2691 1599186245 : 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), & ! LCOV_EXCL_LINE
2700 : fcondtopn_f(n), & ! LCOV_EXCL_LINE
2701 : fsurfn_f (n), & ! LCOV_EXCL_LINE
2702 : flatn (n), fsensn (n), & ! LCOV_EXCL_LINE
2703 : fsurfn (n), & ! LCOV_EXCL_LINE
2704 24093747 : fcondtopn (n))
2705 24093747 : if (icepack_warnings_aborted(subname)) return
2706 : endif
2707 :
2708 1599186245 : if (snwgrain) then
2709 203601768 : rsnw (:) = rsnwn (:,n)
2710 203601768 : smice(:) = smicen(:,n)
2711 203601768 : smliq(:) = smliqn(:,n)
2712 : endif
2713 :
2714 : call thermo_vertical(dt=dt, aicen=aicen (n), &
2715 : vicen=vicen (n), vsnon=vsnon (n), & ! LCOV_EXCL_LINE
2716 : Tsf=Tsfc (n), zSin=zSin (:,n), & ! LCOV_EXCL_LINE
2717 : zqin=zqin (:,n), zqsn=zqsn (:,n), & ! LCOV_EXCL_LINE
2718 : apond=apond (n), hpond=hpnd (n), & ! LCOV_EXCL_LINE
2719 : flw=flw, potT=potT, & ! LCOV_EXCL_LINE
2720 : Qa=Qa, rhoa=rhoa, & ! LCOV_EXCL_LINE
2721 : fsnow=fsnow, fpond=fpond, & ! LCOV_EXCL_LINE
2722 : fbot=fbot, Tbot=Tbot, & ! LCOV_EXCL_LINE
2723 : Tsnice=Tsnice, sss=sss, & ! LCOV_EXCL_LINE
2724 : sst=sst, rsnw=rsnw, & ! LCOV_EXCL_LINE
2725 : lhcoef=lhcoef, shcoef=shcoef, & ! LCOV_EXCL_LINE
2726 : fswsfc=fswsfcn (n), fswint=fswintn (n), & ! LCOV_EXCL_LINE
2727 : Sswabs=Sswabsn(:,n), Iswabs=Iswabsn (:,n), & ! LCOV_EXCL_LINE
2728 : fsurfn=fsurfn (n), fcondtopn=fcondtopn (n), & ! LCOV_EXCL_LINE
2729 : fcondbotn=fcondbotn (n), & ! LCOV_EXCL_LINE
2730 : fsensn=fsensn (n), flatn=flatn (n), & ! LCOV_EXCL_LINE
2731 : flwoutn=flwoutn, evapn=evapn, & ! LCOV_EXCL_LINE
2732 : evapsn=evapsn, evapin=evapin, & ! LCOV_EXCL_LINE
2733 : freshn=freshn, fsaltn=fsaltn, & ! LCOV_EXCL_LINE
2734 : fhocnn=fhocnn, frain=frain, & ! LCOV_EXCL_LINE
2735 : meltt=melttn (n), melts=meltsn (n), & ! LCOV_EXCL_LINE
2736 : meltb=meltbn (n), meltsliq=l_meltsliqn(n), & ! LCOV_EXCL_LINE
2737 : smice=smice, massice=massicen (:,n), & ! LCOV_EXCL_LINE
2738 : smliq=smliq, massliq=massliqn (:,n), & ! LCOV_EXCL_LINE
2739 : congel=congeln (n), snoice=snoicen (n), & ! LCOV_EXCL_LINE
2740 : mlt_onset=mlt_onset, frz_onset=frz_onset , & ! LCOV_EXCL_LINE
2741 : yday=yday, dsnow=l_dsnown , & ! LCOV_EXCL_LINE
2742 1599186245 : prescribed_ice=prescribed_ice)
2743 :
2744 1599186245 : 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 1599186245 : if (tr_pond_lvl) then
2752 1404007548 : if (alvl(n) > puny) then
2753 1347129798 : apnd(n) = max(apond(n) / alvl(n), c1)
2754 : else
2755 56877750 : apnd(n) = c0
2756 : endif
2757 : else
2758 195178697 : apnd(n) = apond(n)
2759 : endif
2760 :
2761 1599186245 : if (snwgrain) then
2762 203601768 : rsnwn (:,n) = rsnw (:)
2763 203601768 : smicen(:,n) = smice(:)
2764 203601768 : smliqn(:,n) = smliq(:)
2765 : endif
2766 :
2767 : !-----------------------------------------------------------------
2768 : ! Total absorbed shortwave radiation
2769 : !-----------------------------------------------------------------
2770 :
2771 1599186245 : fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)
2772 :
2773 : !-----------------------------------------------------------------
2774 : ! Aerosol update
2775 : !-----------------------------------------------------------------
2776 :
2777 1599186245 : if (tr_aero) then
2778 : call update_aerosol (dt, &
2779 : melttn (n), meltsn (n), & ! LCOV_EXCL_LINE
2780 : meltbn (n), congeln (n), & ! LCOV_EXCL_LINE
2781 : snoicen (n), fsnow, & ! LCOV_EXCL_LINE
2782 : aerosno(:,:,n), aeroice(:,:,n), & ! LCOV_EXCL_LINE
2783 : aicen_init (n), vicen_init (n), & ! LCOV_EXCL_LINE
2784 : vsnon_init (n), & ! LCOV_EXCL_LINE
2785 : vicen (n), vsnon (n), & ! LCOV_EXCL_LINE
2786 : aicen (n), & ! LCOV_EXCL_LINE
2787 43852465 : faero_atm , faero_ocn)
2788 43852465 : if (icepack_warnings_aborted(subname)) return
2789 : endif
2790 :
2791 1599186245 : if (tr_iso) then
2792 : call update_isotope (dt = dt, &
2793 : meltt = melttn(n),melts = meltsn(n), & ! LCOV_EXCL_LINE
2794 : meltb = meltbn(n),congel=congeln(n), & ! LCOV_EXCL_LINE
2795 : snoice=snoicen(n),evap=evapn, & ! LCOV_EXCL_LINE
2796 : fsnow=fsnow, Tsfc=Tsfc(n), & ! LCOV_EXCL_LINE
2797 : Qref_iso=Qrefn_iso(:), & ! LCOV_EXCL_LINE
2798 : isosno=isosno(:,n),isoice=isoice(:,n), & ! LCOV_EXCL_LINE
2799 : aice_old=aicen_init(n), & ! LCOV_EXCL_LINE
2800 : vice_old=vicen_init(n), & ! LCOV_EXCL_LINE
2801 : vsno_old=vsnon_init(n), & ! LCOV_EXCL_LINE
2802 : vicen=vicen(n), & ! LCOV_EXCL_LINE
2803 : vsnon=vsnon(n), & ! LCOV_EXCL_LINE
2804 : aicen=aicen(n), & ! LCOV_EXCL_LINE
2805 : fiso_atm=fiso_atm(:), & ! LCOV_EXCL_LINE
2806 : fiso_evapn=fiso_evapn(:), & ! LCOV_EXCL_LINE
2807 : fiso_ocnn=fiso_ocnn(:), & ! LCOV_EXCL_LINE
2808 : HDO_ocn=HDO_ocn, & ! LCOV_EXCL_LINE
2809 : H2_16O_ocn=H2_16O_ocn, & ! LCOV_EXCL_LINE
2810 15774367 : H2_18O_ocn=H2_18O_ocn)
2811 15774367 : if (icepack_warnings_aborted(subname)) return
2812 : endif
2813 :
2814 : endif ! aicen_init
2815 :
2816 7438469352 : if (snwgrain) then
2817 : call drain_snow (vsnon = vsnon(n), &
2818 : aicen = aicen(n), & ! LCOV_EXCL_LINE
2819 : massice = massicen(:,n), & ! LCOV_EXCL_LINE
2820 : massliq = massliqn(:,n), & ! LCOV_EXCL_LINE
2821 263998560 : meltsliq = l_meltsliqn(n))
2822 263998560 : 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 7438469352 : if (tr_pond) then
2833 :
2834 7151215032 : if (tr_pond_lvl) then
2835 7046304408 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2836 : call compute_ponds_lvl (dt=dt, &
2837 : rfrac=rfrac, & ! LCOV_EXCL_LINE
2838 : meltt=melttn (n), & ! LCOV_EXCL_LINE
2839 : melts=meltsn (n), & ! LCOV_EXCL_LINE
2840 : frain=frain, & ! LCOV_EXCL_LINE
2841 : Tair=Tair, & ! LCOV_EXCL_LINE
2842 : fsurfn=fsurfn(n), & ! LCOV_EXCL_LINE
2843 : dhs=dhsn (n), & ! LCOV_EXCL_LINE
2844 : ffrac=ffracn (n), & ! LCOV_EXCL_LINE
2845 : aicen=aicen (n), & ! LCOV_EXCL_LINE
2846 : vicen=vicen (n), & ! LCOV_EXCL_LINE
2847 : vsnon=vsnon (n), & ! LCOV_EXCL_LINE
2848 : qicen=zqin (:,n), & ! LCOV_EXCL_LINE
2849 : sicen=zSin (:,n), & ! LCOV_EXCL_LINE
2850 : Tsfcn=Tsfc (n), & ! LCOV_EXCL_LINE
2851 : alvl=alvl (n), & ! LCOV_EXCL_LINE
2852 : apnd=apnd (n), & ! LCOV_EXCL_LINE
2853 : hpnd=hpnd (n), & ! LCOV_EXCL_LINE
2854 : ipnd=ipnd (n), & ! LCOV_EXCL_LINE
2855 7046304408 : meltsliqn=l_meltsliqn(n))
2856 7046304408 : if (icepack_warnings_aborted(subname)) return
2857 :
2858 104910624 : elseif (tr_pond_topo) then
2859 104910624 : if (aicen_init(n) > puny) then
2860 :
2861 : ! collect liquid water in ponds
2862 : ! assume salt still runs off
2863 24093747 : rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
2864 24093747 : 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 & ! LCOV_EXCL_LINE
2870 24093747 : + 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 24093747 : if (apnd(n) < puny) then
2876 20393084 : hpnd(n) = c0
2877 20393084 : apnd(n) = c1
2878 : endif
2879 24093747 : hpnd(n) = (pond + hpnd(n)*apnd(n)) / apnd(n)
2880 24093747 : 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 7438469352 : if (aicen_init(n) > puny) then
2892 :
2893 1599186245 : l_fswthrun_vdr = c0
2894 1599186245 : l_fswthrun_vdf = c0
2895 1599186245 : l_fswthrun_idr = c0
2896 1599186245 : l_fswthrun_idf = c0
2897 1599186245 : if (present(fswthrun_vdr)) l_fswthrun_vdr = fswthrun_vdr(n)
2898 1599186245 : if (present(fswthrun_vdf)) l_fswthrun_vdf = fswthrun_vdf(n)
2899 1599186245 : if (present(fswthrun_idr)) l_fswthrun_idr = fswthrun_idr(n)
2900 1599186245 : if (present(fswthrun_idf)) l_fswthrun_idf = fswthrun_idf(n)
2901 :
2902 : call merge_fluxes (aicen=aicen_init(n), &
2903 : flw=flw, & ! LCOV_EXCL_LINE
2904 : strairxn=strairxn, strairyn=strairyn,& ! LCOV_EXCL_LINE
2905 : Cdn_atm_ratio_n=Cdn_atm_ratio_n, & ! LCOV_EXCL_LINE
2906 : fsurfn=fsurfn(n), fcondtopn=fcondtopn(n),& ! LCOV_EXCL_LINE
2907 : fcondbotn=fcondbotn(n), & ! LCOV_EXCL_LINE
2908 : fsensn=fsensn(n), flatn=flatn(n), & ! LCOV_EXCL_LINE
2909 : fswabsn=fswabsn, flwoutn=flwoutn, & ! LCOV_EXCL_LINE
2910 : evapn=evapn, & ! LCOV_EXCL_LINE
2911 : evapsn=evapsn, evapin=evapin, & ! LCOV_EXCL_LINE
2912 : Trefn=Trefn, Qrefn=Qrefn, & ! LCOV_EXCL_LINE
2913 : freshn=freshn, fsaltn=fsaltn, & ! LCOV_EXCL_LINE
2914 : fhocnn=fhocnn, & ! LCOV_EXCL_LINE
2915 : fswthrun=fswthrun(n), & ! LCOV_EXCL_LINE
2916 : fswthrun_vdr=l_fswthrun_vdr, & ! LCOV_EXCL_LINE
2917 : fswthrun_vdf=l_fswthrun_vdf, & ! LCOV_EXCL_LINE
2918 : fswthrun_idr=l_fswthrun_idr, & ! LCOV_EXCL_LINE
2919 : fswthrun_idf=l_fswthrun_idf, & ! LCOV_EXCL_LINE
2920 : strairxT=strairxT, strairyT=strairyT,& ! LCOV_EXCL_LINE
2921 : Cdn_atm_ratio=Cdn_atm_ratio, & ! LCOV_EXCL_LINE
2922 : fsurf=fsurf, fcondtop=fcondtop,& ! LCOV_EXCL_LINE
2923 : fcondbot=fcondbot, & ! LCOV_EXCL_LINE
2924 : fsens=fsens, flat=flat, & ! LCOV_EXCL_LINE
2925 : fswabs=fswabs, flwout=flwout, & ! LCOV_EXCL_LINE
2926 : evap=evap, & ! LCOV_EXCL_LINE
2927 : evaps=evaps, evapi=evapi, & ! LCOV_EXCL_LINE
2928 : Tref=Tref, Qref=Qref, & ! LCOV_EXCL_LINE
2929 : fresh=fresh, fsalt=fsalt, & ! LCOV_EXCL_LINE
2930 : fhocn=fhocn, & ! LCOV_EXCL_LINE
2931 : fswthru=fswthru, & ! LCOV_EXCL_LINE
2932 : fswthru_vdr=fswthru_vdr, & ! LCOV_EXCL_LINE
2933 : fswthru_vdf=fswthru_vdf, & ! LCOV_EXCL_LINE
2934 : fswthru_idr=fswthru_idr, & ! LCOV_EXCL_LINE
2935 : fswthru_idf=fswthru_idf, & ! LCOV_EXCL_LINE
2936 : melttn=melttn (n), meltsn=meltsn(n), & ! LCOV_EXCL_LINE
2937 : meltbn=meltbn (n), congeln=congeln(n),& ! LCOV_EXCL_LINE
2938 : meltt=meltt, melts=melts, & ! LCOV_EXCL_LINE
2939 : meltb=meltb, snoicen=snoicen(n),& ! LCOV_EXCL_LINE
2940 : dsnow=l_dsnow, dsnown=l_dsnown, & ! LCOV_EXCL_LINE
2941 : congel=congel, snoice=snoice, & ! LCOV_EXCL_LINE
2942 : meltsliq=l_meltsliq, & ! LCOV_EXCL_LINE
2943 : meltsliqn=l_meltsliqn(n), & ! LCOV_EXCL_LINE
2944 : Uref=Uref, Urefn=Urefn, & ! LCOV_EXCL_LINE
2945 : Qref_iso=Qref_iso, & ! LCOV_EXCL_LINE
2946 : Qrefn_iso=Qrefn_iso, & ! LCOV_EXCL_LINE
2947 : fiso_ocn=fiso_ocn, & ! LCOV_EXCL_LINE
2948 : fiso_ocnn=fiso_ocnn, & ! LCOV_EXCL_LINE
2949 : fiso_evap=fiso_evap, & ! LCOV_EXCL_LINE
2950 1599186245 : fiso_evapn=fiso_evapn)
2951 :
2952 1599186245 : if (icepack_warnings_aborted(subname)) return
2953 :
2954 : endif
2955 :
2956 8930659392 : if (present(dsnown )) dsnown(n) = l_dsnown
2957 :
2958 : enddo ! ncat
2959 :
2960 : !-----------------------------------------------------------------
2961 : ! reload snow mass tracers
2962 : !-----------------------------------------------------------------
2963 :
2964 1492190040 : if (snwgrain) then
2965 316798272 : do n = 1, ncat
2966 316798272 : if (vsnon(n) > puny) then
2967 28904011 : workc = real(nslyr, dbl_kind) * aicen(n) / vsnon(n)
2968 173424066 : do k = 1, nslyr
2969 144520055 : smicen(k,n) = massicen(k,n) * workc
2970 173424066 : smliqn(k,n) = massliqn(k,n) * workc
2971 : enddo
2972 : else ! reset to default values
2973 1410567294 : do k = 1, nslyr
2974 1175472745 : smicen(k,n) = rhos
2975 1410567294 : smliqn(k,n) = c0
2976 : enddo
2977 : endif
2978 : enddo
2979 : endif
2980 :
2981 8921052192 : if (present(meltsliqn )) meltsliqn = l_meltsliqn
2982 1492190040 : if (present(meltsliq )) meltsliq = l_meltsliq
2983 1492190040 : if (present(dsnow )) dsnow = l_dsnow
2984 :
2985 : !-----------------------------------------------------------------
2986 : ! Calculate ponds from the topographic scheme
2987 : !-----------------------------------------------------------------
2988 : !call ice_timer_start(timer_ponds)
2989 1492190040 : if (tr_pond_topo) then
2990 : call compute_ponds_topo(dt, &
2991 : ktherm, & ! LCOV_EXCL_LINE
2992 : aice, aicen, & ! LCOV_EXCL_LINE
2993 : vice, vicen, & ! LCOV_EXCL_LINE
2994 : vsno, vsnon, & ! LCOV_EXCL_LINE
2995 : meltt, & ! LCOV_EXCL_LINE
2996 : fsurf, fpond, & ! LCOV_EXCL_LINE
2997 : Tsfc, Tf, & ! LCOV_EXCL_LINE
2998 : zqin, zSin, & ! LCOV_EXCL_LINE
2999 17485104 : apnd, hpnd, ipnd )
3000 17485104 : if (icepack_warnings_aborted(subname)) return
3001 : endif
3002 : !call ice_timer_stop(timer_ponds)
3003 :
3004 1492190040 : first_call = .false.
3005 :
3006 : end subroutine icepack_step_therm1
3007 :
3008 : !=======================================================================
3009 :
3010 : end module icepack_therm_vertical
3011 :
3012 : !=======================================================================
|
