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