Line data Source code
1 : !=======================================================================
2 :
3 : module icepack_therm_mushy
4 :
5 : use icepack_kinds
6 : use icepack_parameters, only: c0, c1, c2, c8, c10
7 : use icepack_parameters, only: p01, p05, p1, p2, p5, pi, bignum, puny
8 : use icepack_parameters, only: viscosity_dyn, rhow, rhoi, rhos, cp_ocn, cp_ice, Lfresh, gravit, rhofresh
9 : use icepack_parameters, only: hs_min, snwgrain
10 : use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode, tscale_pnd_drain
11 : use icepack_parameters, only: aspect_rapid_mode, dSdt_slow_mode, phi_c_slow_mode
12 : use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp
13 : use icepack_tracers, only: nilyr, nslyr, tr_pond
14 : use icepack_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, icepack_enthalpy_snow
15 : use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction
16 : use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction
17 : use icepack_mushy_physics, only: temperature_snow, temperature_mush_liquid_fraction
18 : use icepack_mushy_physics, only: liquidus_brine_salinity_mush, liquidus_temperature_mush
19 : use icepack_mushy_physics, only: conductivity_mush_array, conductivity_snow_array
20 : use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf
21 : use icepack_therm_shared, only: ferrmax
22 : use icepack_warnings, only: warnstr, icepack_warnings_add
23 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
24 :
25 : implicit none
26 :
27 : private
28 : public :: &
29 : temperature_changes_salinity, & ! LCOV_EXCL_LINE
30 : permeability
31 :
32 : real(kind=dbl_kind), parameter :: &
33 : dTemp_errmax = 5.0e-4_dbl_kind ! max allowed change in temperature
34 : ! between iterations
35 :
36 : !=======================================================================
37 :
38 : contains
39 :
40 : !=======================================================================
41 :
42 1408028216 : subroutine temperature_changes_salinity(dt, &
43 : rhoa, flw, & ! LCOV_EXCL_LINE
44 : potT, Qa, & ! LCOV_EXCL_LINE
45 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
46 : fswsfc, fswint, & ! LCOV_EXCL_LINE
47 1408028216 : Sswabs, Iswabs, & ! LCOV_EXCL_LINE
48 : hilyr, hslyr, & ! LCOV_EXCL_LINE
49 : apond, hpond, & ! LCOV_EXCL_LINE
50 2816056432 : zqin, zTin, & ! LCOV_EXCL_LINE
51 2816056432 : zqsn, zTsn, & ! LCOV_EXCL_LINE
52 1408028216 : zSin, & ! LCOV_EXCL_LINE
53 : Tsf, Tbot, & ! LCOV_EXCL_LINE
54 : sss, & ! LCOV_EXCL_LINE
55 : fsensn, flatn, & ! LCOV_EXCL_LINE
56 : flwoutn, fsurfn, & ! LCOV_EXCL_LINE
57 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
58 : fadvheat, snoice, & ! LCOV_EXCL_LINE
59 1408028216 : smice, smliq)
60 :
61 : ! solve the enthalpy and bulk salinity of the ice for a single column
62 :
63 : real (kind=dbl_kind), intent(in) :: &
64 : dt ! time step (s)
65 :
66 : real (kind=dbl_kind), intent(in) :: &
67 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
68 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
69 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
70 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
71 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
72 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
73 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
74 : sss ! sea surface salinity (PSU)
75 :
76 : real (kind=dbl_kind), intent(inout) :: &
77 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
78 : fswint ! SW absorbed in ice interior below surface (W m-2)
79 :
80 : real (kind=dbl_kind), intent(inout) :: &
81 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
82 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
83 : apond , & ! melt pond area fraction of category ! LCOV_EXCL_LINE
84 : hpond ! melt pond depth (m)
85 :
86 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
87 : Sswabs , & ! SW radiation absorbed in snow layers (W m-2) ! LCOV_EXCL_LINE
88 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
89 : smice , & ! ice mass tracer in snow (kg/m^3) ! LCOV_EXCL_LINE
90 : smliq ! liquid water mass tracer in snow (kg/m^3)
91 :
92 : real (kind=dbl_kind), intent(inout):: &
93 : fsurfn , & ! net flux to top surface, excluding fcondtopn ! LCOV_EXCL_LINE
94 : fcondtop , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
95 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
96 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
97 : flwoutn ! upward LW at surface (W m-2)
98 :
99 : real (kind=dbl_kind), intent(out):: &
100 : fcondbot , & ! downward cond flux at bottom surface (W m-2) ! LCOV_EXCL_LINE
101 : fadvheat , & ! flow of heat to ocean due to advection (W m-2) ! LCOV_EXCL_LINE
102 : snoice ! snow ice formation
103 :
104 : real (kind=dbl_kind), intent(inout):: &
105 : Tsf ! ice/snow surface temperature (C)
106 :
107 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
108 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
109 : zTin , & ! internal ice layer temperatures ! LCOV_EXCL_LINE
110 : zSin , & ! internal ice layer salinities ! LCOV_EXCL_LINE
111 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
112 : zTsn ! internal snow layer temperatures
113 :
114 : ! local variables
115 : real(kind=dbl_kind), dimension(1:nilyr) :: &
116 2816056432 : zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
117 2816056432 : zSin0 , & ! internal ice layer salinities (ppt) at start of timestep ! LCOV_EXCL_LINE
118 2816056432 : phi , & ! liquid fraction ! LCOV_EXCL_LINE
119 2816056432 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
120 1408028216 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
121 :
122 : real(kind=dbl_kind), dimension(1:nilyr+1) :: &
123 2816056432 : Sbr , & ! brine salinity (ppt) ! LCOV_EXCL_LINE
124 2816056432 : qbr ! brine enthalpy (J m-3)
125 :
126 : real(kind=dbl_kind), dimension(0:nilyr) :: &
127 2816056432 : q ! upward interface vertical Darcy flow (m s-1)
128 :
129 : real(kind=dbl_kind), dimension(1:nslyr) :: &
130 2816056432 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
131 1408028216 : ks ! snow conductivity (W m-1 K-1)
132 :
133 : real(kind=dbl_kind) :: &
134 : Tsf0 , & ! ice/snow surface temperature (C) at start of timestep ! LCOV_EXCL_LINE
135 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
136 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
137 : hslyr_min , & ! minimum snow layer thickness (m) ! LCOV_EXCL_LINE
138 : w , & ! vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
139 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
140 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
141 : Spond ! melt pond salinity (ppt)
142 :
143 : real(kind=dbl_kind) :: Tmlt, Iswabs_tmp, dt_rhoi_hlyr, ci, dswabs
144 :
145 : integer(kind=int_kind) :: &
146 : k ! ice/snow layer index
147 :
148 : logical(kind=log_kind) :: &
149 : lsnow ! snow presence: T: has snow, F: no snow
150 :
151 : character(len=*),parameter :: subname='(temperature_changes_salinity)'
152 :
153 1408028216 : fadvheat = c0
154 1408028216 : snoice = c0
155 :
156 1408028216 : Tsf0 = Tsf
157 2999325374 : zqsn0 = zqsn
158 11264225728 : zqin0 = zqin
159 11264225728 : zSin0 = zSin
160 :
161 1408028216 : Spond = c0
162 1408028216 : qpond = enthalpy_brine(c0)
163 :
164 1408028216 : hslyr_min = hs_min / real(nslyr, dbl_kind)
165 :
166 1408028216 : lsnow = (hslyr > hslyr_min)
167 :
168 1408028216 : hin = hilyr * real(nilyr,dbl_kind)
169 :
170 1408028216 : qocn = enthalpy_brine(Tbot)
171 :
172 1408028216 : if (lsnow) then
173 1211582191 : hsn = hslyr * real(nslyr,dbl_kind)
174 : else
175 196446025 : hsn = c0
176 : endif
177 :
178 11264225728 : do k = 1, nilyr
179 11264225728 : phi(k) = icepack_mushy_liquid_fraction(icepack_mushy_temperature_mush(zqin(k),zSin(k)),zSin(k))
180 : enddo ! k
181 :
182 : ! calculate vertical bulk darcy flow
183 : call flushing_velocity(zTin, phi, &
184 : hin, hsn, & ! LCOV_EXCL_LINE
185 : hilyr, & ! LCOV_EXCL_LINE
186 : hpond, apond, & ! LCOV_EXCL_LINE
187 1408028216 : dt, w)
188 1408028216 : if (icepack_warnings_aborted(subname)) return
189 :
190 : ! calculate quantities related to drainage
191 : call explicit_flow_velocities(zSin, &
192 : zTin, Tsf, & ! LCOV_EXCL_LINE
193 : Tbot, q, & ! LCOV_EXCL_LINE
194 : dSdt, Sbr, & ! LCOV_EXCL_LINE
195 : qbr, dt, & ! LCOV_EXCL_LINE
196 : sss, qocn, & ! LCOV_EXCL_LINE
197 1408028216 : hilyr, hin)
198 1408028216 : if (icepack_warnings_aborted(subname)) return
199 :
200 : ! calculate the conductivities
201 1408028216 : call conductivity_mush_array(zqin0, zSin0, km)
202 1408028216 : if (icepack_warnings_aborted(subname)) return
203 :
204 : !-----------------------------------------------------------------
205 : ! Check for excessive absorbed solar radiation that may result in
206 : ! temperature overshoots. Convergence is particularly difficult
207 : ! if the starting temperature is already very close to the melting
208 : ! temperature and extra energy is added. In that case, or if the
209 : ! amount of energy absorbed is greater than the amount needed to
210 : ! melt through a given fraction of a layer, we put the extra
211 : ! energy into the surface.
212 : ! NOTE: This option is not available if the atmosphere model
213 : ! has already computed fsurf. (Unless we adjust fsurf here)
214 : !-----------------------------------------------------------------
215 : !mclaren: Should there be an if calc_Tsfc statement here then??
216 :
217 1408028216 : dswabs = c0
218 1408028216 : if (sw_redist) then
219 0 : dt_rhoi_hlyr = dt / (rhoi*hilyr)
220 0 : do k = 1, nilyr
221 0 : Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc
222 0 : Tmlt = liquidus_temperature_mush(zSin(k))
223 :
224 0 : if (zTin(k) <= Tmlt - sw_dtemp) then
225 0 : ci = cp_ice - Lfresh * Tmlt / (zTin(k)**2)
226 : Iswabs_tmp = min(Iswabs(k), &
227 0 : sw_frac*(Tmlt-zTin(k))*ci/dt_rhoi_hlyr)
228 : endif
229 0 : if (Iswabs_tmp < puny) Iswabs_tmp = c0
230 0 : dswabs = dswabs + min(Iswabs(k) - Iswabs_tmp, fswint)
231 0 : Iswabs(k) = Iswabs_tmp
232 : enddo
233 : endif
234 1408028216 : if (.not. lsnow) then ! hs <= hs_min
235 430390052 : dswabs = dswabs + sum(Sswabs(:))
236 : endif
237 1408028216 : fswsfc = fswsfc + dswabs
238 1408028216 : fswint = fswint - dswabs
239 :
240 1408028216 : if (lsnow) then
241 : ! case with snow
242 :
243 : ! calculate the snow conductivities
244 1211582191 : call conductivity_snow_array(ks)
245 1211582191 : if (icepack_warnings_aborted(subname)) return
246 :
247 : ! run the two stage solver
248 : call two_stage_solver_snow(Tsf, Tsf0, &
249 : zqsn, zqsn0, & ! LCOV_EXCL_LINE
250 : zqin, zqin0, & ! LCOV_EXCL_LINE
251 : zSin, zSin0, & ! LCOV_EXCL_LINE
252 : zTsn, & ! LCOV_EXCL_LINE
253 : zTin, & ! LCOV_EXCL_LINE
254 : phi, Tbot, & ! LCOV_EXCL_LINE
255 : km, ks, & ! LCOV_EXCL_LINE
256 : q, dSdt, & ! LCOV_EXCL_LINE
257 : w, dt, & ! LCOV_EXCL_LINE
258 : fswint, fswsfc, & ! LCOV_EXCL_LINE
259 : rhoa, flw, & ! LCOV_EXCL_LINE
260 : potT, Qa, & ! LCOV_EXCL_LINE
261 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
262 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
263 : qpond, qocn, & ! LCOV_EXCL_LINE
264 : Spond, sss, & ! LCOV_EXCL_LINE
265 : hilyr, hslyr, & ! LCOV_EXCL_LINE
266 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
267 : fadvheat, & ! LCOV_EXCL_LINE
268 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
269 1211582191 : flatn, fsurfn )
270 :
271 1211582191 : if (icepack_warnings_aborted(subname)) then
272 0 : write(warnstr,*) subname, "temperature_changes_salinity: Picard solver non-convergence (snow)"
273 0 : call icepack_warnings_add(warnstr)
274 0 : return
275 : endif
276 :
277 : ! given the updated enthalpy and bulk salinity calculate other quantities
278 2568935322 : do k = 1, nslyr
279 2568935322 : zTsn(k) = temperature_snow(zqsn(k))
280 : enddo ! k
281 :
282 9692657528 : do k = 1, nilyr
283 8481075337 : zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
284 8481075337 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
285 9692657528 : qbr(k) = enthalpy_brine(zTin(k))
286 : enddo ! k
287 :
288 : else
289 : ! case without snow
290 :
291 : ! run the two stage solver
292 : call two_stage_solver_nosnow(Tsf, Tsf0, &
293 : zqsn, & ! LCOV_EXCL_LINE
294 : zqin, zqin0, & ! LCOV_EXCL_LINE
295 : zSin, zSin0, & ! LCOV_EXCL_LINE
296 : zTsn, & ! LCOV_EXCL_LINE
297 : zTin, & ! LCOV_EXCL_LINE
298 : phi, Tbot, & ! LCOV_EXCL_LINE
299 : km, ks, & ! LCOV_EXCL_LINE
300 : q, dSdt, & ! LCOV_EXCL_LINE
301 : w, dt, & ! LCOV_EXCL_LINE
302 : fswint, fswsfc, & ! LCOV_EXCL_LINE
303 : rhoa, flw, & ! LCOV_EXCL_LINE
304 : potT, Qa, & ! LCOV_EXCL_LINE
305 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
306 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
307 : qpond, qocn, & ! LCOV_EXCL_LINE
308 : Spond, sss, & ! LCOV_EXCL_LINE
309 : hilyr, hslyr, & ! LCOV_EXCL_LINE
310 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
311 : fadvheat, & ! LCOV_EXCL_LINE
312 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
313 196446025 : flatn, fsurfn )
314 :
315 196446025 : if (icepack_warnings_aborted(subname)) then
316 0 : write(warnstr,*) subname, "temperature_changes_salinity: Picard solver non-convergence (no snow)"
317 0 : call icepack_warnings_add(warnstr)
318 0 : return
319 : endif
320 :
321 : ! given the updated enthalpy and bulk salinity calculate other quantities
322 1571568200 : do k = 1, nilyr
323 1375122175 : zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
324 1375122175 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
325 1571568200 : qbr(k) = enthalpy_brine(zTin(k))
326 : enddo ! k
327 :
328 : endif
329 :
330 : ! drain ponds from flushing
331 1408028216 : call flush_pond(w, hpond, apond, dt)
332 1408028216 : if (icepack_warnings_aborted(subname)) return
333 :
334 : ! flood snow ice
335 : call flood_ice(hsn, hin, &
336 : hslyr, hilyr, & ! LCOV_EXCL_LINE
337 : zqsn, zqin, & ! LCOV_EXCL_LINE
338 : phi, dt, & ! LCOV_EXCL_LINE
339 : zSin, Sbr, & ! LCOV_EXCL_LINE
340 : sss, qocn, & ! LCOV_EXCL_LINE
341 : smice, smliq, & ! LCOV_EXCL_LINE
342 1408028216 : snoice, fadvheat)
343 1408028216 : if (icepack_warnings_aborted(subname)) return
344 :
345 : end subroutine temperature_changes_salinity
346 :
347 : !=======================================================================
348 :
349 1211582191 : subroutine two_stage_solver_snow(Tsf, Tsf0, &
350 2423164382 : zqsn, zqsn0, & ! LCOV_EXCL_LINE
351 2423164382 : zqin, zqin0, & ! LCOV_EXCL_LINE
352 2423164382 : zSin, zSin0, & ! LCOV_EXCL_LINE
353 0 : zTsn, & ! LCOV_EXCL_LINE
354 1211582191 : zTin, & ! LCOV_EXCL_LINE
355 1211582191 : phi, Tbot, & ! LCOV_EXCL_LINE
356 1211582191 : km, ks, & ! LCOV_EXCL_LINE
357 1211582191 : q, dSdt, & ! LCOV_EXCL_LINE
358 : w, dt, & ! LCOV_EXCL_LINE
359 : fswint, fswsfc, & ! LCOV_EXCL_LINE
360 : rhoa, flw, & ! LCOV_EXCL_LINE
361 : potT, Qa, & ! LCOV_EXCL_LINE
362 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
363 1211582191 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
364 : qpond, qocn, & ! LCOV_EXCL_LINE
365 : Spond, sss, & ! LCOV_EXCL_LINE
366 : hilyr, hslyr, & ! LCOV_EXCL_LINE
367 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
368 : fadvheat, & ! LCOV_EXCL_LINE
369 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
370 : flatn, fsurfn )
371 :
372 : ! solve the vertical temperature and salt change for case with snow
373 : ! 1) determine what type of surface condition existed previously - cold or melting
374 : ! 2) solve the system assuming this condition persists
375 : ! 3) check the consistency of the surface condition of the solution
376 : ! 4) If the surface condition is inconsistent resolve for the other surface condition
377 : ! 5) If neither solution is consistent the resolve the inconsistency
378 :
379 : real(kind=dbl_kind), intent(inout) :: &
380 : Tsf ! snow surface temperature (C)
381 :
382 : real(kind=dbl_kind), intent(out) :: &
383 : fcondtop , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
384 : fcondbot , & ! downward cond flux at bottom surface (W m-2) ! LCOV_EXCL_LINE
385 : flwoutn , & ! upward LW at surface (W m-2) ! LCOV_EXCL_LINE
386 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
387 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
388 : fsurfn , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
389 : fadvheat ! flow of heat to ocean due to advection (W m-2)
390 :
391 : real(kind=dbl_kind), intent(in) :: &
392 : Tsf0 ! snow surface temperature (C) at beginning of timestep
393 :
394 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
395 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
396 : zTsn ! snow layer temperature (C)
397 :
398 : real(kind=dbl_kind), dimension(:), intent(in) :: &
399 : zqsn0 , & ! snow layer enthalpy (J m-3) at beginning of timestep ! LCOV_EXCL_LINE
400 : ks , & ! snow conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
401 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
402 :
403 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
404 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
405 : zSin , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
406 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
407 : phi ! ice layer liquid fraction
408 :
409 : real(kind=dbl_kind), dimension(:), intent(in) :: &
410 : zqin0 , & ! ice layer enthalpy (J m-3) at beginning of timestep ! LCOV_EXCL_LINE
411 : zSin0 , & ! ice layer bulk salinity (ppt) at beginning of timestep ! LCOV_EXCL_LINE
412 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
413 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
414 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
415 :
416 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
417 : q ! upward interface vertical Darcy flow (m s-1)
418 :
419 : real(kind=dbl_kind), intent(in) :: &
420 : dt , & ! time step (s) ! LCOV_EXCL_LINE
421 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
422 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
423 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
424 : fswint , & ! SW absorbed in ice interior below surface (W m-2) ! LCOV_EXCL_LINE
425 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
426 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
427 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
428 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
429 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
430 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
431 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
432 : w , & ! vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
433 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
434 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
435 : Spond , & ! melt pond salinity (ppt) ! LCOV_EXCL_LINE
436 : sss ! sea surface salinity (PSU)
437 :
438 : real(kind=dbl_kind) :: &
439 : fcondtop1 , & ! first stage downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
440 : fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
441 : Tsf1 ! first stage ice surface temperature (C)
442 :
443 : character(len=*),parameter :: subname='(two_stage_solver_snow)'
444 :
445 :
446 : ! determine if surface is initially cold or melting
447 1211582191 : if (Tsf < c0) then
448 :
449 : ! initially cold
450 :
451 : ! solve the system for cold and snow
452 : call picard_solver(.true., .true., &
453 : Tsf, zqsn, & ! LCOV_EXCL_LINE
454 : zqin, zSin, & ! LCOV_EXCL_LINE
455 : zTin, zTsn, & ! LCOV_EXCL_LINE
456 : phi, dt, & ! LCOV_EXCL_LINE
457 : hilyr, hslyr, & ! LCOV_EXCL_LINE
458 : km, ks, & ! LCOV_EXCL_LINE
459 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
460 : Tbot, & ! LCOV_EXCL_LINE
461 : fswint, fswsfc, & ! LCOV_EXCL_LINE
462 : rhoa, flw, & ! LCOV_EXCL_LINE
463 : potT, Qa, & ! LCOV_EXCL_LINE
464 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
465 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
466 : fadvheat, & ! LCOV_EXCL_LINE
467 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
468 : flatn, fsurfn, & ! LCOV_EXCL_LINE
469 : qpond, qocn, & ! LCOV_EXCL_LINE
470 : Spond, sss, & ! LCOV_EXCL_LINE
471 : q, dSdt, & ! LCOV_EXCL_LINE
472 1157663616 : w )
473 :
474 : ! halt if solver failed
475 2369245807 : if (icepack_warnings_aborted(subname)) return
476 :
477 : ! check if solution is consistent - surface should still be cold
478 1157663616 : if (Tsf < dTemp_errmax) then
479 :
480 : ! solution is consistent - have solution so finish
481 1047438353 : return
482 :
483 : else
484 :
485 : ! solution is inconsistent - surface is warmer than melting
486 : ! resolve assuming surface is melting
487 110225263 : Tsf1 = Tsf
488 :
489 : ! reset the solution to initial values
490 110225263 : Tsf = c0
491 231357446 : zqsn = zqsn0
492 881802104 : zqin = zqin0
493 881802104 : zSin = zSin0
494 :
495 : ! solve the system for melting and snow
496 : call picard_solver(.true., .false., &
497 : Tsf, zqsn, & ! LCOV_EXCL_LINE
498 : zqin, zSin, & ! LCOV_EXCL_LINE
499 : zTin, zTsn, & ! LCOV_EXCL_LINE
500 : phi, dt, & ! LCOV_EXCL_LINE
501 : hilyr, hslyr, & ! LCOV_EXCL_LINE
502 : km, ks, & ! LCOV_EXCL_LINE
503 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
504 : Tbot, & ! LCOV_EXCL_LINE
505 : fswint, fswsfc, & ! LCOV_EXCL_LINE
506 : rhoa, flw, & ! LCOV_EXCL_LINE
507 : potT, Qa, & ! LCOV_EXCL_LINE
508 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
509 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
510 : fadvheat, & ! LCOV_EXCL_LINE
511 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
512 : flatn, fsurfn, & ! LCOV_EXCL_LINE
513 : qpond, qocn, & ! LCOV_EXCL_LINE
514 : Spond, sss, & ! LCOV_EXCL_LINE
515 : q, dSdt, & ! LCOV_EXCL_LINE
516 110225263 : w )
517 :
518 : ! halt if solver failed
519 110225263 : if (icepack_warnings_aborted(subname)) return
520 :
521 : ! check if solution is consistent
522 : ! surface conductive heat flux should be less than
523 : ! incoming surface heat flux
524 110225263 : if (fcondtop - fsurfn < ferrmax) then
525 :
526 : ! solution is consistent - have solution so finish
527 110225263 : return
528 :
529 : else
530 :
531 : ! solution is inconsistent
532 0 : call two_stage_inconsistency(1, Tsf1, c0, fcondtop, fsurfn)
533 0 : if (icepack_warnings_aborted(subname)) return
534 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
535 0 : call icepack_warnings_add(subname//" two_stage_solver_snow: two_stage_inconsistency: cold" )
536 0 : return
537 :
538 : endif ! surface flux consistency
539 :
540 : endif ! surface temperature consistency
541 :
542 : else
543 :
544 : ! initially melting
545 53918575 : Tsf = c0
546 :
547 : ! solve the system for melting and snow
548 : call picard_solver(.true., .false., &
549 : Tsf, zqsn, & ! LCOV_EXCL_LINE
550 : zqin, zSin, & ! LCOV_EXCL_LINE
551 : zTin, zTsn, & ! LCOV_EXCL_LINE
552 : phi, dt, & ! LCOV_EXCL_LINE
553 : hilyr, hslyr, & ! LCOV_EXCL_LINE
554 : km, ks, & ! LCOV_EXCL_LINE
555 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
556 : Tbot, & ! LCOV_EXCL_LINE
557 : fswint, fswsfc, & ! LCOV_EXCL_LINE
558 : rhoa, flw, & ! LCOV_EXCL_LINE
559 : potT, Qa, & ! LCOV_EXCL_LINE
560 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
561 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
562 : fadvheat, & ! LCOV_EXCL_LINE
563 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
564 : flatn, fsurfn, & ! LCOV_EXCL_LINE
565 : qpond, qocn, & ! LCOV_EXCL_LINE
566 : Spond, sss, & ! LCOV_EXCL_LINE
567 : q, dSdt, & ! LCOV_EXCL_LINE
568 53918575 : w )
569 :
570 :
571 : ! halt if solver failed
572 53918575 : if (icepack_warnings_aborted(subname)) return
573 :
574 : ! check if solution is consistent
575 : ! surface conductive heat flux should be less than
576 : ! incoming surface heat flux
577 53918575 : if (fcondtop - fsurfn < ferrmax) then
578 :
579 : ! solution is consistent - have solution so finish
580 51254914 : return
581 :
582 : else
583 :
584 : ! solution is inconsistent - resolve assuming other surface condition
585 : ! assume surface is cold
586 2663661 : fcondtop1 = fcondtop
587 2663661 : fsurfn1 = fsurfn
588 :
589 : ! reset the solution to initial values
590 2663661 : Tsf = Tsf0
591 5342982 : zqsn = zqsn0
592 21309288 : zqin = zqin0
593 21309288 : zSin = zSin0
594 :
595 : ! solve the system for cold and snow
596 : call picard_solver(.true., .true., &
597 : Tsf, zqsn, & ! LCOV_EXCL_LINE
598 : zqin, zSin, & ! LCOV_EXCL_LINE
599 : zTin, zTsn, & ! LCOV_EXCL_LINE
600 : phi, dt, & ! LCOV_EXCL_LINE
601 : hilyr, hslyr, & ! LCOV_EXCL_LINE
602 : km, ks, & ! LCOV_EXCL_LINE
603 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
604 : Tbot, & ! LCOV_EXCL_LINE
605 : fswint, fswsfc, & ! LCOV_EXCL_LINE
606 : rhoa, flw, & ! LCOV_EXCL_LINE
607 : potT, Qa, & ! LCOV_EXCL_LINE
608 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
609 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
610 : fadvheat, & ! LCOV_EXCL_LINE
611 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
612 : flatn, fsurfn, & ! LCOV_EXCL_LINE
613 : qpond, qocn, & ! LCOV_EXCL_LINE
614 : Spond, sss, & ! LCOV_EXCL_LINE
615 : q, dSdt, & ! LCOV_EXCL_LINE
616 2663661 : w )
617 :
618 : ! halt if solver failed
619 2663661 : if (icepack_warnings_aborted(subname)) return
620 :
621 : ! check if solution is consistent - surface should be cold
622 2663661 : if (Tsf < dTemp_errmax) then
623 :
624 : ! solution is consistent - have solution so finish
625 2663661 : return
626 :
627 : else
628 :
629 : ! solution is inconsistent
630 : ! failed to find a solution so need to refine solutions until consistency found
631 0 : call two_stage_inconsistency(2, Tsf, c0, fcondtop1, fsurfn1)
632 0 : if (icepack_warnings_aborted(subname)) return
633 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
634 0 : call icepack_warnings_add(subname//" two_stage_solver_snow: two_stage_inconsistency: melting" )
635 0 : return
636 :
637 : endif ! surface temperature consistency
638 :
639 : endif ! surface flux consistency
640 :
641 : endif
642 :
643 : end subroutine two_stage_solver_snow
644 :
645 : !=======================================================================
646 :
647 196446025 : subroutine two_stage_solver_nosnow(Tsf, Tsf0, &
648 0 : zqsn, & ! LCOV_EXCL_LINE
649 392892050 : zqin, zqin0, & ! LCOV_EXCL_LINE
650 392892050 : zSin, zSin0, & ! LCOV_EXCL_LINE
651 0 : zTsn, & ! LCOV_EXCL_LINE
652 196446025 : zTin, & ! LCOV_EXCL_LINE
653 196446025 : phi, Tbot, & ! LCOV_EXCL_LINE
654 196446025 : km, ks, & ! LCOV_EXCL_LINE
655 196446025 : q, dSdt, & ! LCOV_EXCL_LINE
656 : w, dt, & ! LCOV_EXCL_LINE
657 : fswint, fswsfc, & ! LCOV_EXCL_LINE
658 : rhoa, flw, & ! LCOV_EXCL_LINE
659 : potT, Qa, & ! LCOV_EXCL_LINE
660 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
661 196446025 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
662 : qpond, qocn, & ! LCOV_EXCL_LINE
663 : Spond, sss, & ! LCOV_EXCL_LINE
664 : hilyr, hslyr, & ! LCOV_EXCL_LINE
665 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
666 : fadvheat, & ! LCOV_EXCL_LINE
667 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
668 : flatn, fsurfn )
669 :
670 : ! solve the vertical temperature and salt change for case with no snow
671 : ! 1) determine what type of surface condition existed previously - cold or melting
672 : ! 2) solve the system assuming this condition persists
673 : ! 3) check the consistency of the surface condition of the solution
674 : ! 4) If the surface condition is inconsistent resolve for the other surface condition
675 : ! 5) If neither solution is consistent the resolve the inconsistency
676 :
677 : real(kind=dbl_kind), intent(inout) :: &
678 : Tsf ! ice surface temperature (C)
679 :
680 : real(kind=dbl_kind), intent(out) :: &
681 : fcondtop , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
682 : fcondbot , & ! downward cond flux at bottom surface (W m-2) ! LCOV_EXCL_LINE
683 : flwoutn , & ! upward LW at surface (W m-2) ! LCOV_EXCL_LINE
684 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
685 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
686 : fsurfn , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
687 : fadvheat ! flow of heat to ocean due to advection (W m-2)
688 :
689 : real(kind=dbl_kind), intent(in) :: &
690 : Tsf0 ! ice surface temperature (C) at beginning of timestep
691 :
692 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
693 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
694 : zTsn ! snow layer temperature (C)
695 :
696 : real(kind=dbl_kind), dimension(:), intent(in) :: &
697 : ks , & ! snow conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
698 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
699 :
700 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
701 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
702 : zSin , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
703 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
704 : phi ! ice layer liquid fraction
705 :
706 : real(kind=dbl_kind), dimension(:), intent(in) :: &
707 : zqin0 , & ! ice layer enthalpy (J m-3) at beginning of timestep ! LCOV_EXCL_LINE
708 : zSin0 , & ! ice layer bulk salinity (ppt) at beginning of timestep ! LCOV_EXCL_LINE
709 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
710 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
711 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
712 :
713 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
714 : q ! upward interface vertical Darcy flow (m s-1)
715 :
716 : real(kind=dbl_kind), intent(in) :: &
717 : dt , & ! time step (s) ! LCOV_EXCL_LINE
718 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
719 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
720 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
721 : fswint , & ! SW absorbed in ice interior below surface (W m-2) ! LCOV_EXCL_LINE
722 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
723 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
724 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
725 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
726 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
727 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
728 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
729 : w , & ! vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
730 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
731 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
732 : Spond , & ! melt pond salinity (ppt) ! LCOV_EXCL_LINE
733 : sss ! sea surface salinity (PSU)
734 :
735 : real(kind=dbl_kind) :: &
736 : Tmlt , & ! upper ice layer melting temperature (C) ! LCOV_EXCL_LINE
737 : fcondtop1 , & ! first stage downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
738 : fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
739 : Tsf1 ! first stage ice surface temperature (C)
740 :
741 : character(len=*),parameter :: subname='(two_stage_solver_nosnow)'
742 :
743 : ! initial surface melting temperature
744 196446025 : Tmlt = liquidus_temperature_mush(zSin0(1))
745 :
746 : ! determine if surface is initially cold or melting
747 196446025 : if (Tsf < Tmlt) then
748 :
749 : ! initially cold
750 :
751 : ! solve the system for cold and no snow
752 : call picard_solver(.false., .true., &
753 : Tsf, zqsn, & ! LCOV_EXCL_LINE
754 : zqin, zSin, & ! LCOV_EXCL_LINE
755 : zTin, zTsn, & ! LCOV_EXCL_LINE
756 : phi, dt, & ! LCOV_EXCL_LINE
757 : hilyr, hslyr, & ! LCOV_EXCL_LINE
758 : km, ks, & ! LCOV_EXCL_LINE
759 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
760 : Tbot, & ! LCOV_EXCL_LINE
761 : fswint, fswsfc, & ! LCOV_EXCL_LINE
762 : rhoa, flw, & ! LCOV_EXCL_LINE
763 : potT, Qa, & ! LCOV_EXCL_LINE
764 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
765 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
766 : fadvheat, & ! LCOV_EXCL_LINE
767 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
768 : flatn, fsurfn, & ! LCOV_EXCL_LINE
769 : qpond, qocn, & ! LCOV_EXCL_LINE
770 : Spond, sss, & ! LCOV_EXCL_LINE
771 : q, dSdt, & ! LCOV_EXCL_LINE
772 147833986 : w )
773 :
774 : ! halt if solver failed
775 344280011 : if (icepack_warnings_aborted(subname)) return
776 :
777 : ! check if solution is consistent - surface should still be cold
778 147833986 : if (Tsf < Tmlt + dTemp_errmax) then
779 :
780 : ! solution is consistent - have solution so finish
781 93580737 : return
782 :
783 : else
784 : ! solution is inconsistent - surface is warmer than melting
785 : ! resolve assuming surface is melting
786 54253249 : Tsf1 = Tsf
787 :
788 : ! reset the solution to initial values
789 54253249 : Tsf = liquidus_temperature_mush(zSin0(1))
790 434025992 : zqin = zqin0
791 434025992 : zSin = zSin0
792 :
793 : ! solve the system for melt and no snow
794 : call picard_solver(.false., .false., &
795 : Tsf, zqsn, & ! LCOV_EXCL_LINE
796 : zqin, zSin, & ! LCOV_EXCL_LINE
797 : zTin, zTsn, & ! LCOV_EXCL_LINE
798 : phi, dt, & ! LCOV_EXCL_LINE
799 : hilyr, hslyr, & ! LCOV_EXCL_LINE
800 : km, ks, & ! LCOV_EXCL_LINE
801 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
802 : Tbot, & ! LCOV_EXCL_LINE
803 : fswint, fswsfc, & ! LCOV_EXCL_LINE
804 : rhoa, flw, & ! LCOV_EXCL_LINE
805 : potT, Qa, & ! LCOV_EXCL_LINE
806 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
807 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
808 : fadvheat, & ! LCOV_EXCL_LINE
809 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
810 : flatn, fsurfn, & ! LCOV_EXCL_LINE
811 : qpond, qocn, & ! LCOV_EXCL_LINE
812 : Spond, sss, & ! LCOV_EXCL_LINE
813 : q, dSdt, & ! LCOV_EXCL_LINE
814 54253249 : w )
815 54253249 : if (icepack_warnings_aborted(subname)) return
816 :
817 : ! halt if solver failed
818 54253249 : if (icepack_warnings_aborted(subname)) return
819 :
820 : ! check if solution is consistent
821 : ! surface conductive heat flux should be less than
822 : ! incoming surface heat flux
823 54253249 : if (fcondtop - fsurfn < ferrmax) then
824 :
825 : ! solution is consistent - have solution so finish
826 54253249 : return
827 :
828 : else
829 :
830 : ! solution is inconsistent
831 0 : call two_stage_inconsistency(3, Tsf1, Tmlt, fcondtop, fsurfn)
832 0 : if (icepack_warnings_aborted(subname)) return
833 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
834 0 : call icepack_warnings_add(subname//" two_stage_solver_nosnow: two_stage_inconsistency: cold" )
835 0 : return
836 :
837 : endif
838 :
839 : endif
840 :
841 : else
842 : ! initially melting
843 :
844 : ! solve the system for melt and no snow
845 48612039 : Tsf = Tmlt
846 :
847 : call picard_solver(.false., .false., &
848 : Tsf, zqsn, & ! LCOV_EXCL_LINE
849 : zqin, zSin, & ! LCOV_EXCL_LINE
850 : zTin, zTsn, & ! LCOV_EXCL_LINE
851 : phi, dt, & ! LCOV_EXCL_LINE
852 : hilyr, hslyr, & ! LCOV_EXCL_LINE
853 : km, ks, & ! LCOV_EXCL_LINE
854 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
855 : Tbot, & ! LCOV_EXCL_LINE
856 : fswint, fswsfc, & ! LCOV_EXCL_LINE
857 : rhoa, flw, & ! LCOV_EXCL_LINE
858 : potT, Qa, & ! LCOV_EXCL_LINE
859 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
860 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
861 : fadvheat, & ! LCOV_EXCL_LINE
862 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
863 : flatn, fsurfn, & ! LCOV_EXCL_LINE
864 : qpond, qocn, & ! LCOV_EXCL_LINE
865 : Spond, sss, & ! LCOV_EXCL_LINE
866 : q, dSdt, & ! LCOV_EXCL_LINE
867 48612039 : w )
868 48612039 : if (icepack_warnings_aborted(subname)) return
869 :
870 : ! halt if solver failed
871 48612039 : if (icepack_warnings_aborted(subname)) return
872 :
873 : ! check if solution is consistent
874 : ! surface conductive heat flux should be less than
875 : ! incoming surface heat flux
876 48612039 : if (fcondtop - fsurfn < ferrmax) then
877 :
878 : ! solution is consistent - have solution so finish
879 47686001 : return
880 :
881 : else
882 :
883 : ! solution is inconsistent - resolve assuming other surface condition
884 : ! assume surface is cold
885 926038 : fcondtop1 = fcondtop
886 926038 : fsurfn1 = fsurfn
887 :
888 : ! reset the solution to initial values
889 926038 : Tsf = Tsf0
890 7408304 : zqin = zqin0
891 7408304 : zSin = zSin0
892 :
893 : ! solve the system for cold and no snow
894 : call picard_solver(.false., .true., &
895 : Tsf, zqsn, & ! LCOV_EXCL_LINE
896 : zqin, zSin, & ! LCOV_EXCL_LINE
897 : zTin, zTsn, & ! LCOV_EXCL_LINE
898 : phi, dt, & ! LCOV_EXCL_LINE
899 : hilyr, hslyr, & ! LCOV_EXCL_LINE
900 : km, ks, & ! LCOV_EXCL_LINE
901 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
902 : Tbot, & ! LCOV_EXCL_LINE
903 : fswint, fswsfc, & ! LCOV_EXCL_LINE
904 : rhoa, flw, & ! LCOV_EXCL_LINE
905 : potT, Qa, & ! LCOV_EXCL_LINE
906 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
907 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
908 : fadvheat, & ! LCOV_EXCL_LINE
909 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
910 : flatn, fsurfn, & ! LCOV_EXCL_LINE
911 : qpond, qocn, & ! LCOV_EXCL_LINE
912 : Spond, sss, & ! LCOV_EXCL_LINE
913 : q, dSdt, & ! LCOV_EXCL_LINE
914 926038 : w )
915 926038 : if (icepack_warnings_aborted(subname)) return
916 :
917 : ! halt if solver failed
918 926038 : if (icepack_warnings_aborted(subname)) return
919 :
920 : ! check if solution is consistent - surface should be cold
921 926038 : if (Tsf < Tmlt + dTemp_errmax) then
922 :
923 : ! solution is consistent - have solution so finish
924 926038 : return
925 :
926 : else
927 :
928 : ! solution is inconsistent
929 0 : call two_stage_inconsistency(4, Tsf, Tmlt, fcondtop1, fsurfn1)
930 0 : if (icepack_warnings_aborted(subname)) return
931 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
932 0 : call icepack_warnings_add(subname//" two_stage_solver_nosnow: two_stage_inconsistency: melting" )
933 0 : return
934 :
935 : endif
936 :
937 : endif
938 :
939 : endif
940 :
941 : end subroutine two_stage_solver_nosnow
942 :
943 : !=======================================================================
944 :
945 0 : subroutine two_stage_inconsistency(type, Tsf, Tmlt, fcondtop, fsurfn)
946 :
947 : integer (kind=int_kind), intent(in) :: &
948 : type
949 :
950 : real(kind=dbl_kind), intent(in) :: &
951 : Tsf, & ! LCOV_EXCL_LINE
952 : Tmlt, & ! LCOV_EXCL_LINE
953 : fcondtop, & ! LCOV_EXCL_LINE
954 : fsurfn
955 :
956 : character(len=*),parameter :: subname='(two_stage_inconsistency)'
957 :
958 0 : write(warnstr,*) subname, "icepack_therm_mushy: two stage inconsistency"
959 0 : call icepack_warnings_add(warnstr)
960 0 : write(warnstr,*) subname, "type:", type
961 0 : call icepack_warnings_add(warnstr)
962 :
963 0 : if (type == 1) then
964 :
965 0 : write(warnstr,*) subname, "First stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
966 0 : call icepack_warnings_add(warnstr)
967 0 : write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
968 0 : call icepack_warnings_add(warnstr)
969 0 : write(warnstr,*) subname, "Second stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
970 0 : call icepack_warnings_add(warnstr)
971 0 : write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
972 0 : call icepack_warnings_add(warnstr)
973 :
974 0 : else if (type == 2) then
975 :
976 0 : write(warnstr,*) subname, "First stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
977 0 : call icepack_warnings_add(warnstr)
978 0 : write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
979 0 : call icepack_warnings_add(warnstr)
980 0 : write(warnstr,*) subname, "Second stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
981 0 : call icepack_warnings_add(warnstr)
982 0 : write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
983 0 : call icepack_warnings_add(warnstr)
984 :
985 0 : else if (type == 3) then
986 :
987 0 : write(warnstr,*) subname, "First stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
988 0 : call icepack_warnings_add(warnstr)
989 0 : write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
990 0 : call icepack_warnings_add(warnstr)
991 0 : write(warnstr,*) subname, "Second stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
992 0 : call icepack_warnings_add(warnstr)
993 0 : write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
994 0 : call icepack_warnings_add(warnstr)
995 :
996 0 : else if (type == 4) then
997 :
998 0 : write(warnstr,*) subname, "First stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
999 0 : call icepack_warnings_add(warnstr)
1000 0 : write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
1001 0 : call icepack_warnings_add(warnstr)
1002 0 : write(warnstr,*) subname, "Second stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
1003 0 : call icepack_warnings_add(warnstr)
1004 0 : write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
1005 0 : call icepack_warnings_add(warnstr)
1006 :
1007 : endif
1008 :
1009 0 : end subroutine two_stage_inconsistency
1010 :
1011 : !=======================================================================
1012 : ! Picard/TDMA based solver
1013 : !=======================================================================
1014 :
1015 0 : subroutine prep_picard(lsnow, zqsn, &
1016 1576096427 : zqin, zSin, & ! LCOV_EXCL_LINE
1017 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1018 1576096427 : km, ks, & ! LCOV_EXCL_LINE
1019 3152192854 : zTin, zTsn, & ! LCOV_EXCL_LINE
1020 1576096427 : Sbr, phi, & ! LCOV_EXCL_LINE
1021 3152192854 : dxp, kcstar, & ! LCOV_EXCL_LINE
1022 : einit)
1023 :
1024 : logical, intent(in) :: &
1025 : lsnow ! snow presence: T: has snow, F: no snow
1026 :
1027 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1028 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1029 : zSin , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
1030 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
1031 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1032 : ks ! snow conductivity (W m-1 K-1)
1033 :
1034 : real(kind=dbl_kind), intent(in) :: &
1035 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1036 : hslyr ! snow layer thickness (m)
1037 :
1038 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1039 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1040 : Sbr , & ! ice layer brine salinity (ppt) ! LCOV_EXCL_LINE
1041 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
1042 : zTsn , & ! snow layer temperature (C) ! LCOV_EXCL_LINE
1043 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
1044 : kcstar ! interface conductivities (W m-1 K-1)
1045 :
1046 : real(kind=dbl_kind), intent(out) :: &
1047 : einit ! initial total energy (J)
1048 :
1049 : integer(kind=int_kind) :: k
1050 :
1051 : character(len=*),parameter :: subname='(prep_picard)'
1052 :
1053 : ! calculate initial ice temperatures
1054 12608771416 : do k = 1, nilyr
1055 11032674989 : zTin(k) = icepack_mushy_temperature_mush(zqin(k), zSin(k))
1056 11032674989 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
1057 12608771416 : phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
1058 : enddo ! k
1059 :
1060 1576096427 : if (lsnow) then
1061 :
1062 2805635750 : do k = 1, nslyr
1063 2805635750 : zTsn(k) = temperature_snow(zqsn(k))
1064 : enddo ! k
1065 :
1066 : endif ! lsnow
1067 :
1068 : ! interface distances
1069 1576096427 : call calc_intercell_thickness(lsnow, hilyr, hslyr, dxp)
1070 1576096427 : if (icepack_warnings_aborted(subname)) return
1071 :
1072 : ! interface conductivities
1073 1576096427 : call calc_intercell_conductivity(lsnow, km, ks, hilyr, hslyr, kcstar)
1074 1576096427 : if (icepack_warnings_aborted(subname)) return
1075 :
1076 : ! total energy content
1077 : call total_energy_content(lsnow, &
1078 : zqin, zqsn, & ! LCOV_EXCL_LINE
1079 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1080 1576096427 : einit)
1081 1576096427 : if (icepack_warnings_aborted(subname)) return
1082 :
1083 : end subroutine prep_picard
1084 :
1085 : !=======================================================================
1086 :
1087 1576096427 : subroutine picard_solver(lsnow, lcold, &
1088 1576096427 : Tsf, zqsn, & ! LCOV_EXCL_LINE
1089 1576096427 : zqin, zSin, & ! LCOV_EXCL_LINE
1090 1576096427 : zTin, zTsn, & ! LCOV_EXCL_LINE
1091 1576096427 : phi, dt, & ! LCOV_EXCL_LINE
1092 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1093 1576096427 : km, ks, & ! LCOV_EXCL_LINE
1094 1576096427 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
1095 : Tbot, & ! LCOV_EXCL_LINE
1096 : fswint, fswsfc, & ! LCOV_EXCL_LINE
1097 : rhoa, flw, & ! LCOV_EXCL_LINE
1098 : potT, Qa, & ! LCOV_EXCL_LINE
1099 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
1100 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
1101 : fadvheat, & ! LCOV_EXCL_LINE
1102 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
1103 : flatn, fsurfn, & ! LCOV_EXCL_LINE
1104 : qpond, qocn, & ! LCOV_EXCL_LINE
1105 : Spond, sss, & ! LCOV_EXCL_LINE
1106 3152192854 : q, dSdt, & ! LCOV_EXCL_LINE
1107 : w )
1108 :
1109 : logical, intent(in) :: &
1110 : lsnow , & ! snow presence: T: has snow, F: no snow ! LCOV_EXCL_LINE
1111 : lcold ! surface cold: T: surface is cold, F: surface is melting
1112 :
1113 : real(kind=dbl_kind), intent(inout) :: &
1114 : Tsf ! snow surface temperature (C)
1115 :
1116 : real(kind=dbl_kind), intent(out) :: &
1117 : fcondtop , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
1118 : fcondbot , & ! downward cond flux at bottom surface (W m-2) ! LCOV_EXCL_LINE
1119 : fadvheat ! flow of heat to ocean due to advection (W m-2)
1120 :
1121 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1122 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1123 : zSin , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
1124 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1125 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
1126 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1127 : zTsn ! snow layer temperature (C)
1128 :
1129 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1130 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
1131 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
1132 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
1133 :
1134 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
1135 : q ! upward interface vertical Darcy flow (m s-1)
1136 :
1137 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1138 : ks , & ! snow conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
1139 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
1140 :
1141 : real(kind=dbl_kind), intent(out) :: &
1142 : flwoutn , & ! upward LW at surface (W m-2) ! LCOV_EXCL_LINE
1143 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
1144 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
1145 : fsurfn ! net flux to top surface, excluding fcondtop
1146 :
1147 : real(kind=dbl_kind), intent(in) :: &
1148 : dt , & ! time step (s) ! LCOV_EXCL_LINE
1149 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1150 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
1151 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
1152 : fswint , & ! SW absorbed in ice interior below surface (W m-2) ! LCOV_EXCL_LINE
1153 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
1154 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
1155 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
1156 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
1157 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
1158 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
1159 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
1160 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
1161 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
1162 : Spond , & ! melt pond salinity (ppt) ! LCOV_EXCL_LINE
1163 : sss , & ! sea surface salinity (ppt) ! LCOV_EXCL_LINE
1164 : w ! vertical flushing Darcy velocity (m/s)
1165 :
1166 : real(kind=dbl_kind), dimension(nilyr) :: &
1167 3152192854 : Sbr , & ! ice layer brine salinity (ppt) ! LCOV_EXCL_LINE
1168 3152192854 : qbr , & ! ice layer brine enthalpy (J m-3) ! LCOV_EXCL_LINE
1169 3152192854 : zTin0 , & ! ice layer temperature (C) at start of timestep ! LCOV_EXCL_LINE
1170 3152192854 : zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
1171 3152192854 : zSin0 , & ! ice layer bulk salinity (ppt) at start of timestep ! LCOV_EXCL_LINE
1172 3152192854 : zTin_prev ! ice layer temperature at previous iteration
1173 :
1174 : real(kind=dbl_kind), dimension(nslyr) :: &
1175 3152192854 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
1176 3152192854 : zTsn0 , & ! snow layer temperature (C) at start of timestep ! LCOV_EXCL_LINE
1177 3152192854 : zTsn_prev ! snow layer temperature at previous iteration
1178 :
1179 : real(kind=dbl_kind), dimension(nslyr+nilyr+1) :: &
1180 1576096427 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
1181 1576096427 : kcstar ! interface conductivities (W m-1 K-1)
1182 :
1183 : real(kind=dbl_kind) :: &
1184 : Tsf0 , & ! snow surface temperature (C) at start of timestep ! LCOV_EXCL_LINE
1185 : dfsurfn_dTsf , & ! derivative of net flux to top surface, excluding fcondtopn ! LCOV_EXCL_LINE
1186 : dflwoutn_dTsf , & ! derivative of longwave flux wrt surface temperature ! LCOV_EXCL_LINE
1187 : dfsensn_dTsf , & ! derivative of sensible heat flux wrt surface temperature ! LCOV_EXCL_LINE
1188 : dflatn_dTsf , & ! derivative of latent heat flux wrt surface temperature ! LCOV_EXCL_LINE
1189 : Tsf_prev , & ! snow surface temperature at previous iteration ! LCOV_EXCL_LINE
1190 : einit , & ! initial total energy (J) ! LCOV_EXCL_LINE
1191 : fadvheat_nit ! heat to ocean due to advection (W m-2) during iteration
1192 :
1193 : logical :: &
1194 : lconverged ! has Picard solver converged?
1195 :
1196 : integer :: &
1197 : nit ! Picard iteration count
1198 :
1199 : integer, parameter :: &
1200 : nit_max = 100 ! maximum number of Picard iterations
1201 :
1202 : character(len=*),parameter :: subname='(picard_solver)'
1203 :
1204 1576096427 : lconverged = .false.
1205 :
1206 : ! prepare quantities for picard iteration
1207 : call prep_picard(lsnow, zqsn, &
1208 : zqin, zSin, & ! LCOV_EXCL_LINE
1209 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1210 : km, ks, & ! LCOV_EXCL_LINE
1211 : zTin, zTsn, & ! LCOV_EXCL_LINE
1212 : Sbr, phi, & ! LCOV_EXCL_LINE
1213 : dxp, kcstar, & ! LCOV_EXCL_LINE
1214 1576096427 : einit)
1215 1576096427 : if (icepack_warnings_aborted(subname)) return
1216 :
1217 1576096427 : Tsf0 = Tsf
1218 12608771416 : zqin0 = zqin
1219 3349671340 : zqsn0 = zqsn
1220 12608771416 : zTin0 = zTin
1221 3349671340 : zTsn0 = zTsn
1222 12608771416 : zSin0 = zSin
1223 :
1224 : ! set prev variables
1225 1576096427 : Tsf_prev = Tsf
1226 3349671340 : zTsn_prev = zTsn
1227 12608771416 : zTin_prev = zTin
1228 :
1229 : ! picard iteration
1230 3618924854 : picard: do nit = 1, nit_max
1231 :
1232 : ! surface heat flux
1233 : call surface_heat_flux(Tsf, fswsfc, &
1234 : rhoa, flw, & ! LCOV_EXCL_LINE
1235 : potT, Qa, & ! LCOV_EXCL_LINE
1236 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
1237 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
1238 3618924854 : flatn, fsurfn)
1239 3618924854 : if (icepack_warnings_aborted(subname)) return
1240 :
1241 : ! derivative of heat flux with respect to surface temperature
1242 : call dsurface_heat_flux_dTsf(Tsf, rhoa, &
1243 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
1244 : dfsurfn_dTsf, dflwoutn_dTsf, & ! LCOV_EXCL_LINE
1245 3618924854 : dfsensn_dTsf, dflatn_dTsf)
1246 3618924854 : if (icepack_warnings_aborted(subname)) return
1247 :
1248 : ! tridiagonal solve of new temperatures
1249 : call solve_heat_conduction(lsnow, lcold, &
1250 : Tsf, Tbot, & ! LCOV_EXCL_LINE
1251 : zqin0, zqsn0, & ! LCOV_EXCL_LINE
1252 : phi, dt, & ! LCOV_EXCL_LINE
1253 : qpond, qocn, & ! LCOV_EXCL_LINE
1254 : q, w, & ! LCOV_EXCL_LINE
1255 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1256 : dxp, kcstar, & ! LCOV_EXCL_LINE
1257 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
1258 : fsurfn, dfsurfn_dTsf, & ! LCOV_EXCL_LINE
1259 3618924854 : zTin, zTsn)
1260 3618924854 : if (icepack_warnings_aborted(subname)) return
1261 :
1262 : ! update brine enthalpy
1263 3618924854 : call picard_updates_enthalpy(zTin, qbr)
1264 3618924854 : if (icepack_warnings_aborted(subname)) return
1265 :
1266 : ! drainage fluxes
1267 : call picard_drainage_fluxes(fadvheat_nit, q, &
1268 3618924854 : qbr, qocn)
1269 3618924854 : if (icepack_warnings_aborted(subname)) return
1270 :
1271 : ! flushing fluxes
1272 : call picard_flushing_fluxes(fadvheat_nit, w, &
1273 : qbr, & ! LCOV_EXCL_LINE
1274 3618924854 : qpond)
1275 3618924854 : if (icepack_warnings_aborted(subname)) return
1276 :
1277 : ! perform convergence check
1278 : call check_picard_convergence(lsnow, lconverged, &
1279 : Tsf, Tsf_prev, & ! LCOV_EXCL_LINE
1280 : zTin, zTin_prev,& ! LCOV_EXCL_LINE
1281 : zTsn, zTsn_prev,& ! LCOV_EXCL_LINE
1282 : phi, Tbot, & ! LCOV_EXCL_LINE
1283 : zqin, zqsn, & ! LCOV_EXCL_LINE
1284 : km, ks, & ! LCOV_EXCL_LINE
1285 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1286 : fswint, & ! LCOV_EXCL_LINE
1287 : einit, dt, & ! LCOV_EXCL_LINE
1288 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
1289 3618924854 : fadvheat_nit)
1290 3618924854 : if (icepack_warnings_aborted(subname)) return
1291 :
1292 3618924854 : if (lconverged) exit
1293 :
1294 2042828427 : Tsf_prev = Tsf
1295 4333854840 : zTsn_prev = zTsn
1296 19961552270 : zTin_prev = zTin
1297 :
1298 : enddo picard
1299 :
1300 1576096427 : fadvheat = fadvheat_nit
1301 :
1302 : ! update the picard iterants
1303 1576096427 : call picard_updates(zTin, Sbr, qbr)
1304 1576096427 : if (icepack_warnings_aborted(subname)) return
1305 :
1306 : ! solve for the salinity
1307 : call solve_salinity(zSin, Sbr, &
1308 : Spond, sss, & ! LCOV_EXCL_LINE
1309 : q, dSdt, & ! LCOV_EXCL_LINE
1310 : w, hilyr, & ! LCOV_EXCL_LINE
1311 1576096427 : dt)
1312 1576096427 : if (icepack_warnings_aborted(subname)) return
1313 :
1314 : ! final surface heat flux
1315 : call surface_heat_flux(Tsf, fswsfc, &
1316 : rhoa, flw, & ! LCOV_EXCL_LINE
1317 : potT, Qa, & ! LCOV_EXCL_LINE
1318 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
1319 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
1320 1576096427 : flatn, fsurfn)
1321 1576096427 : if (icepack_warnings_aborted(subname)) return
1322 :
1323 : ! if not converged
1324 1576096427 : if (.not. lconverged) then
1325 :
1326 : call picard_nonconvergence(Tsf0, Tsf, &
1327 : zTsn0, zTsn, & ! LCOV_EXCL_LINE
1328 : zTin0, zTin, & ! LCOV_EXCL_LINE
1329 : zSin0, zSin, & ! LCOV_EXCL_LINE
1330 : zqsn0, zqsn, & ! LCOV_EXCL_LINE
1331 : zqin0, phi, & ! LCOV_EXCL_LINE
1332 : dt, & ! LCOV_EXCL_LINE
1333 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1334 : km, ks, & ! LCOV_EXCL_LINE
1335 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
1336 : Tbot, & ! LCOV_EXCL_LINE
1337 : fswint, fswsfc, & ! LCOV_EXCL_LINE
1338 : rhoa, flw, & ! LCOV_EXCL_LINE
1339 : potT, Qa, & ! LCOV_EXCL_LINE
1340 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
1341 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
1342 : fadvheat, & ! LCOV_EXCL_LINE
1343 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
1344 : flatn, fsurfn, & ! LCOV_EXCL_LINE
1345 : qpond, qocn, & ! LCOV_EXCL_LINE
1346 : Spond, sss, & ! LCOV_EXCL_LINE
1347 : q, dSdt, & ! LCOV_EXCL_LINE
1348 0 : w)
1349 :
1350 0 : if (icepack_warnings_aborted(subname)) return
1351 0 : call icepack_warnings_add(subname//" picard_solver: Picard solver non-convergence" )
1352 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1353 0 : if (icepack_warnings_aborted(subname)) return
1354 :
1355 : endif
1356 :
1357 : end subroutine picard_solver
1358 :
1359 : !=======================================================================
1360 :
1361 0 : subroutine picard_nonconvergence(Tsf0, Tsf, &
1362 0 : zTsn0, zTsn, & ! LCOV_EXCL_LINE
1363 0 : zTin0, zTin, & ! LCOV_EXCL_LINE
1364 0 : zSin0, zSin, & ! LCOV_EXCL_LINE
1365 0 : zqsn0, zqsn, & ! LCOV_EXCL_LINE
1366 0 : zqin0, phi, & ! LCOV_EXCL_LINE
1367 : dt, & ! LCOV_EXCL_LINE
1368 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1369 0 : km, ks, & ! LCOV_EXCL_LINE
1370 0 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
1371 : Tbot, & ! LCOV_EXCL_LINE
1372 : fswint, fswsfc, & ! LCOV_EXCL_LINE
1373 : rhoa, flw, & ! LCOV_EXCL_LINE
1374 : potT, Qa, & ! LCOV_EXCL_LINE
1375 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
1376 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
1377 : fadvheat, & ! LCOV_EXCL_LINE
1378 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
1379 : flatn, fsurfn, & ! LCOV_EXCL_LINE
1380 : qpond, qocn, & ! LCOV_EXCL_LINE
1381 : Spond, sss, & ! LCOV_EXCL_LINE
1382 0 : q, dSdt, & ! LCOV_EXCL_LINE
1383 : w)
1384 :
1385 : real(kind=dbl_kind), intent(in) :: &
1386 : Tsf0 , & ! snow surface temperature (C) at beginning of timestep ! LCOV_EXCL_LINE
1387 : Tsf ! snow surface temperature (C)
1388 :
1389 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1390 : zTsn0 , & ! snow layer temperature (C) at beginning of timestep ! LCOV_EXCL_LINE
1391 : zTsn , & ! snow layer temperature (C) ! LCOV_EXCL_LINE
1392 : zqsn0 , & ! LCOV_EXCL_LINE
1393 : zqsn , & ! LCOV_EXCL_LINE
1394 : zTin0 , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1395 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1396 : zSin0 , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
1397 : zSin , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
1398 : zqin0 , & ! LCOV_EXCL_LINE
1399 : phi ! ice layer liquid fraction
1400 :
1401 : real(kind=dbl_kind), intent(in) :: &
1402 : dt , & ! time step (s) ! LCOV_EXCL_LINE
1403 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1404 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
1405 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
1406 : fswint , & ! SW absorbed in ice interior below surface (W m-2) ! LCOV_EXCL_LINE
1407 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
1408 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
1409 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
1410 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
1411 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
1412 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
1413 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
1414 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
1415 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
1416 : Spond , & ! melt pond salinity (ppt) ! LCOV_EXCL_LINE
1417 : sss , & ! sea surface salinity (ppt) ! LCOV_EXCL_LINE
1418 : w ! vertical flushing Darcy velocity (m/s)
1419 :
1420 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1421 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
1422 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
1423 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
1424 :
1425 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
1426 : q ! upward interface vertical Darcy flow (m s-1)
1427 :
1428 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1429 : ks , & ! snow conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
1430 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
1431 :
1432 : real(kind=dbl_kind), intent(in) :: &
1433 : flwoutn , & ! upward LW at surface (W m-2) ! LCOV_EXCL_LINE
1434 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
1435 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
1436 : fsurfn ! net flux to top surface, excluding fcondtop
1437 :
1438 : real(kind=dbl_kind), intent(in) :: &
1439 : fcondtop , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
1440 : fcondbot , & ! downward cond flux at bottom surface (W m-2) ! LCOV_EXCL_LINE
1441 : fadvheat ! flow of heat to ocean due to advection (W m-2)
1442 :
1443 : integer :: &
1444 : k ! vertical layer index
1445 :
1446 : character(len=char_len_long) :: &
1447 : warning ! warning message
1448 :
1449 : character(len=*),parameter :: subname='(picard_nonconvergence)'
1450 :
1451 0 : write(warning,*) "-------------------------------------"
1452 0 : call icepack_warnings_add(warning)
1453 0 : write(warning,*)
1454 0 : call icepack_warnings_add(warning)
1455 :
1456 0 : write(warning,*) trim(subname)//":picard convergence failed!"
1457 0 : call icepack_warnings_add(warning)
1458 0 : write(warning,*) "=========================="
1459 0 : call icepack_warnings_add(warning)
1460 0 : write(warning,*)
1461 0 : call icepack_warnings_add(warning)
1462 :
1463 0 : write(warning,*) "Surface: Tsf0, Tsf"
1464 0 : call icepack_warnings_add(warning)
1465 0 : write(warning,*) 0, Tsf0, Tsf
1466 0 : call icepack_warnings_add(warning)
1467 0 : write(warning,*)
1468 0 : call icepack_warnings_add(warning)
1469 :
1470 0 : write(warning,*) "Snow: zTsn0(k), zTsn(k), zqsn0(k), ks(k), Sswabs(k)"
1471 0 : call icepack_warnings_add(warning)
1472 0 : do k = 1, nslyr
1473 0 : write(warning,*) k, zTsn0(k), zTsn(k), zqsn0(k), ks(k), Sswabs(k)
1474 0 : call icepack_warnings_add(warning)
1475 : enddo ! k
1476 0 : write(warning,*)
1477 0 : call icepack_warnings_add(warning)
1478 :
1479 0 : write(warning,*) "Ice: zTin0(k), zTin(k), zSin0(k), zSin(k), phi(k), zqin0(k), km(k), Iswabs(k), dSdt(k)"
1480 0 : call icepack_warnings_add(warning)
1481 0 : do k = 1, nilyr
1482 0 : write(warning,*) k, zTin0(k), zTin(k), zSin0(k), zSin(k), phi(k), zqin0(k), km(k), Iswabs(k), dSdt(k)
1483 0 : call icepack_warnings_add(warning)
1484 : enddo ! k
1485 0 : write(warning,*)
1486 0 : call icepack_warnings_add(warning)
1487 :
1488 0 : write(warning,*) "Ice boundary: q(k)"
1489 0 : call icepack_warnings_add(warning)
1490 0 : do k = 0, nilyr
1491 0 : write(warning,*) k, q(k)
1492 0 : call icepack_warnings_add(warning)
1493 : enddo ! k
1494 0 : write(warning,*)
1495 0 : call icepack_warnings_add(warning)
1496 :
1497 0 : write(warning,*) "dt: ", dt
1498 0 : call icepack_warnings_add(warning)
1499 0 : write(warning,*) "hilyr: ", hilyr
1500 0 : call icepack_warnings_add(warning)
1501 0 : write(warning,*) "hslyr: ", hslyr
1502 0 : call icepack_warnings_add(warning)
1503 0 : write(warning,*) "Tbot: ", Tbot
1504 0 : call icepack_warnings_add(warning)
1505 0 : write(warning,*) "fswint: ", fswint
1506 0 : call icepack_warnings_add(warning)
1507 0 : write(warning,*) "fswsfc: ", fswsfc
1508 0 : call icepack_warnings_add(warning)
1509 0 : write(warning,*) "rhoa: ", rhoa
1510 0 : call icepack_warnings_add(warning)
1511 0 : write(warning,*) "flw: ", flw
1512 0 : call icepack_warnings_add(warning)
1513 0 : write(warning,*) "potT: ", potT
1514 0 : call icepack_warnings_add(warning)
1515 0 : write(warning,*) "Qa: ", Qa
1516 0 : call icepack_warnings_add(warning)
1517 0 : write(warning,*) "shcoef: ", shcoef
1518 0 : call icepack_warnings_add(warning)
1519 0 : write(warning,*) "lhcoef: ", lhcoef
1520 0 : call icepack_warnings_add(warning)
1521 0 : write(warning,*) "qpond: ", qpond
1522 0 : call icepack_warnings_add(warning)
1523 0 : write(warning,*) "qocn: ", qocn
1524 0 : call icepack_warnings_add(warning)
1525 0 : write(warning,*) "Spond: ", Spond
1526 0 : call icepack_warnings_add(warning)
1527 0 : write(warning,*) "sss: ", sss
1528 0 : call icepack_warnings_add(warning)
1529 0 : write(warning,*) "w: ", w
1530 0 : call icepack_warnings_add(warning)
1531 0 : write(warning,*) "flwoutn: ", flwoutn
1532 0 : call icepack_warnings_add(warning)
1533 0 : write(warning,*) "fsensn: ", fsensn
1534 0 : call icepack_warnings_add(warning)
1535 0 : write(warning,*) "flatn: ", flatn
1536 0 : call icepack_warnings_add(warning)
1537 0 : write(warning,*) "fsurfn: ", fsurfn
1538 0 : call icepack_warnings_add(warning)
1539 0 : write(warning,*) "fcondtop: ", fcondtop
1540 0 : call icepack_warnings_add(warning)
1541 0 : write(warning,*) "fcondbot: ", fcondbot
1542 0 : call icepack_warnings_add(warning)
1543 0 : write(warning,*) "fadvheat: ", fadvheat
1544 0 : call icepack_warnings_add(warning)
1545 0 : write(warning,*)
1546 0 : call icepack_warnings_add(warning)
1547 :
1548 0 : write(warning,*) "-------------------------------------"
1549 0 : call icepack_warnings_add(warning)
1550 :
1551 0 : end subroutine picard_nonconvergence
1552 :
1553 : !=======================================================================
1554 :
1555 3618924854 : subroutine check_picard_convergence(lsnow, &
1556 : lconverged, & ! LCOV_EXCL_LINE
1557 : Tsf, Tsf_prev, & ! LCOV_EXCL_LINE
1558 0 : zTin, zTin_prev,& ! LCOV_EXCL_LINE
1559 7237849708 : zTsn, zTsn_prev,& ! LCOV_EXCL_LINE
1560 3618924854 : phi, Tbot, & ! LCOV_EXCL_LINE
1561 3618924854 : zqin, zqsn, & ! LCOV_EXCL_LINE
1562 3618924854 : km, ks, & ! LCOV_EXCL_LINE
1563 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1564 : fswint, & ! LCOV_EXCL_LINE
1565 : einit, dt, & ! LCOV_EXCL_LINE
1566 : fcondtop, fcondbot, & ! LCOV_EXCL_LINE
1567 : fadvheat)
1568 :
1569 : logical, intent(inout) :: &
1570 : lconverged ! has Picard solver converged?
1571 :
1572 : logical, intent(in) :: &
1573 : lsnow ! snow presence: T: has snow, F: no snow
1574 :
1575 : real(kind=dbl_kind), intent(in) :: &
1576 : dt , & ! time step (s) ! LCOV_EXCL_LINE
1577 : Tsf , & ! snow surface temperature (C) ! LCOV_EXCL_LINE
1578 : Tsf_prev , & ! snow surface temperature at previous iteration ! LCOV_EXCL_LINE
1579 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1580 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
1581 : fswint , & ! SW absorbed in ice interior below surface (W m-2) ! LCOV_EXCL_LINE
1582 : einit , & ! initial total energy (J) ! LCOV_EXCL_LINE
1583 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
1584 : fadvheat ! flow of heat to ocean due to advection (W m-2)
1585 :
1586 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1587 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1588 : zTin_prev, & ! ice layer temperature at previous iteration ! LCOV_EXCL_LINE
1589 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
1590 : km ! ice conductivity (W m-1 K-1)
1591 :
1592 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1593 : zqin ! ice layer enthalpy (J m-3)
1594 :
1595 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1596 : zTsn , & ! snow layer temperature (C) ! LCOV_EXCL_LINE
1597 : zTsn_prev, & ! snow layer temperature at previous iteration ! LCOV_EXCL_LINE
1598 : ks ! snow conductivity (W m-1 K-1)
1599 :
1600 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1601 : zqsn ! snow layer enthalpy (J m-3)
1602 :
1603 : real(kind=dbl_kind), intent(out) :: &
1604 : fcondtop , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
1605 : fcondbot ! downward cond flux at bottom surface (W m-2)
1606 :
1607 : real(kind=dbl_kind) :: &
1608 : ferr , & ! energy flux error ! LCOV_EXCL_LINE
1609 : efinal , & ! initial total energy (J) at iteration ! LCOV_EXCL_LINE
1610 : dzTsn , & ! change in snow temperature (C) between iterations ! LCOV_EXCL_LINE
1611 : dzTin , & ! change in ice temperature (C) between iterations ! LCOV_EXCL_LINE
1612 : dTsf ! change in surface temperature (C) between iterations
1613 :
1614 : character(len=*),parameter :: subname='(check_picard_convergence)'
1615 :
1616 : call picard_final(lsnow, &
1617 : zqin, zqsn, & ! LCOV_EXCL_LINE
1618 : zTin, zTsn, & ! LCOV_EXCL_LINE
1619 3618924854 : phi)
1620 3618924854 : if (icepack_warnings_aborted(subname)) return
1621 :
1622 : call total_energy_content(lsnow, &
1623 : zqin, zqsn, & ! LCOV_EXCL_LINE
1624 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1625 3618924854 : efinal)
1626 3618924854 : if (icepack_warnings_aborted(subname)) return
1627 :
1628 : call maximum_variables_changes(lsnow, &
1629 : Tsf, Tsf_prev, dTsf, & ! LCOV_EXCL_LINE
1630 : zTsn, zTsn_prev, dzTsn, & ! LCOV_EXCL_LINE
1631 3618924854 : zTin, zTin_prev, dzTin)
1632 3618924854 : if (icepack_warnings_aborted(subname)) return
1633 :
1634 3618924854 : fcondbot = c2 * km(nilyr) * ((zTin(nilyr) - Tbot) / hilyr)
1635 :
1636 3618924854 : if (lsnow) then
1637 3048334638 : fcondtop = c2 * ks(1) * ((Tsf - zTsn(1)) / hslyr)
1638 : else
1639 570590216 : fcondtop = c2 * km(1) * ((Tsf - zTin(1)) / hilyr)
1640 : endif
1641 :
1642 3618924854 : ferr = (efinal - einit) / dt - (fcondtop - fcondbot + fswint - fadvheat)
1643 :
1644 : lconverged = (dTsf < dTemp_errmax .and. &
1645 : dzTsn < dTemp_errmax .and. & ! LCOV_EXCL_LINE
1646 : dzTin < dTemp_errmax .and. & ! LCOV_EXCL_LINE
1647 3618924854 : abs(ferr) < 0.9_dbl_kind*ferrmax)
1648 :
1649 : end subroutine check_picard_convergence
1650 :
1651 : !=======================================================================
1652 :
1653 3618924854 : subroutine picard_drainage_fluxes(fadvheat, q, &
1654 3618924854 : qbr, qocn)
1655 :
1656 : real(kind=dbl_kind), intent(out) :: &
1657 : fadvheat ! flow of heat to ocean due to advection (W m-2)
1658 :
1659 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
1660 : q ! upward interface vertical Darcy flow (m s-1)
1661 :
1662 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1663 : qbr ! ice layer brine enthalpy (J m-3)
1664 :
1665 : real(kind=dbl_kind), intent(in) :: &
1666 : qocn ! ocean brine enthalpy (J m-3)
1667 :
1668 : integer :: &
1669 : k ! vertical layer index
1670 :
1671 : character(len=*),parameter :: subname='(picard_drainage_fluxes)'
1672 :
1673 3618924854 : fadvheat = c0
1674 :
1675 : ! calculate fluxes from base upwards
1676 25332473978 : do k = 1, nilyr-1
1677 :
1678 25332473978 : fadvheat = fadvheat - q(k) * (qbr(k+1) - qbr(k))
1679 :
1680 : enddo ! k
1681 :
1682 3618924854 : k = nilyr
1683 :
1684 3618924854 : fadvheat = fadvheat - q(k) * (qocn - qbr(k))
1685 :
1686 3618924854 : end subroutine picard_drainage_fluxes
1687 :
1688 : !=======================================================================
1689 :
1690 3618924854 : subroutine picard_flushing_fluxes(fadvheat, w, &
1691 3618924854 : qbr, & ! LCOV_EXCL_LINE
1692 : qpond)
1693 :
1694 : real(kind=dbl_kind), intent(inout) :: &
1695 : fadvheat ! flow of heat to ocean due to advection (W m-2)
1696 :
1697 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1698 : qbr ! ice layer brine enthalpy (J m-3)
1699 :
1700 : real(kind=dbl_kind), intent(in) :: &
1701 : w , & ! vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
1702 : qpond ! melt pond brine enthalpy (J m-3)
1703 :
1704 : character(len=*),parameter :: subname='(picard_flushing_fluxes)'
1705 :
1706 3618924854 : fadvheat = fadvheat + w * (qbr(nilyr) - qpond)
1707 :
1708 3618924854 : end subroutine picard_flushing_fluxes
1709 :
1710 : !=======================================================================
1711 :
1712 3618924854 : subroutine maximum_variables_changes(lsnow, &
1713 : Tsf, Tsf_prev, dTsf, & ! LCOV_EXCL_LINE
1714 3618924854 : zTsn, zTsn_prev, dzTsn, & ! LCOV_EXCL_LINE
1715 3618924854 : zTin, zTin_prev, dzTin)
1716 :
1717 : logical, intent(in) :: &
1718 : lsnow ! snow presence: T: has snow, F: no snow
1719 :
1720 : real(kind=dbl_kind), intent(in) :: &
1721 : Tsf , & ! snow surface temperature (C) ! LCOV_EXCL_LINE
1722 : Tsf_prev ! snow surface temperature at previous iteration
1723 :
1724 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1725 : zTsn , & ! snow layer temperature (C) ! LCOV_EXCL_LINE
1726 : zTsn_prev, & ! snow layer temperature at previous iteration ! LCOV_EXCL_LINE
1727 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1728 : zTin_prev ! ice layer temperature at previous iteration
1729 :
1730 : real(kind=dbl_kind), intent(out) :: &
1731 : dTsf , & ! change in surface temperature (C) between iterations ! LCOV_EXCL_LINE
1732 : dzTsn , & ! change in snow temperature (C) between iterations ! LCOV_EXCL_LINE
1733 : dzTin ! change in surface temperature (C) between iterations
1734 :
1735 : character(len=*),parameter :: subname='(maximum_variables_changes)'
1736 :
1737 3618924854 : dTsf = abs(Tsf - Tsf_prev)
1738 :
1739 3618924854 : if (lsnow) then
1740 6452515890 : dzTsn = maxval(abs(zTsn - zTsn_prev))
1741 : else ! lsnow
1742 570590216 : dzTsn = c0
1743 : endif ! lsnow
1744 :
1745 28951398832 : dzTin = maxval(abs(zTin - zTin_prev))
1746 :
1747 3618924854 : end subroutine maximum_variables_changes
1748 :
1749 : !=======================================================================
1750 :
1751 5195021281 : subroutine total_energy_content(lsnow, &
1752 5195021281 : zqin, zqsn, & ! LCOV_EXCL_LINE
1753 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1754 : energy)
1755 :
1756 : logical, intent(in) :: &
1757 : lsnow ! snow presence: T: has snow, F: no snow
1758 :
1759 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1760 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1761 : zqsn ! snow layer enthalpy (J m-3)
1762 :
1763 : real(kind=dbl_kind), intent(in) :: &
1764 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1765 : hslyr ! snow layer thickness (m)
1766 :
1767 : real(kind=dbl_kind), intent(out) :: &
1768 : energy ! total energy of ice and snow
1769 :
1770 : integer :: &
1771 : k ! vertical layer index
1772 :
1773 : character(len=*),parameter :: subname='(total_energy_content)'
1774 :
1775 5195021281 : energy = c0
1776 :
1777 5195021281 : if (lsnow) then
1778 :
1779 9258151640 : do k = 1, nslyr
1780 :
1781 9258151640 : energy = energy + hslyr * zqsn(k)
1782 :
1783 : enddo ! k
1784 :
1785 : endif ! lsnow
1786 :
1787 41560170248 : do k = 1, nilyr
1788 :
1789 41560170248 : energy = energy + hilyr * zqin(k)
1790 :
1791 : enddo ! k
1792 :
1793 5195021281 : end subroutine total_energy_content
1794 :
1795 : !=======================================================================
1796 :
1797 0 : subroutine picard_updates(zTin, &
1798 3152192854 : Sbr, qbr)
1799 :
1800 : ! update brine salinity and liquid fraction based on new temperatures
1801 :
1802 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1803 : zTin ! ice layer temperature (C)
1804 :
1805 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1806 : Sbr , & ! ice layer brine salinity (ppt) ! LCOV_EXCL_LINE
1807 : qbr ! ice layer brine enthalpy (J m-3)
1808 :
1809 : integer :: &
1810 : k ! vertical layer index
1811 :
1812 : character(len=*),parameter :: subname='(picard_updates)'
1813 :
1814 12608771416 : do k = 1, nilyr
1815 :
1816 11032674989 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
1817 12608771416 : qbr(k) = enthalpy_brine(zTin(k))
1818 :
1819 : enddo ! k
1820 :
1821 1576096427 : end subroutine picard_updates
1822 :
1823 : !=======================================================================
1824 :
1825 3618924854 : subroutine picard_updates_enthalpy(zTin, qbr)
1826 :
1827 : ! update brine salinity and liquid fraction based on new temperatures
1828 :
1829 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1830 : zTin ! ice layer temperature (C)
1831 :
1832 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1833 : qbr ! ice layer brine enthalpy (J m-3)
1834 :
1835 : integer :: &
1836 : k ! vertical layer index
1837 :
1838 : character(len=*),parameter :: subname='(picard_updates_enthalpy)'
1839 :
1840 28951398832 : do k = 1, nilyr
1841 :
1842 28951398832 : qbr(k) = enthalpy_brine(zTin(k))
1843 :
1844 : enddo ! k
1845 :
1846 3618924854 : end subroutine picard_updates_enthalpy
1847 :
1848 : !=======================================================================
1849 :
1850 3618924854 : subroutine picard_final(lsnow, &
1851 7237849708 : zqin, zqsn, & ! LCOV_EXCL_LINE
1852 7237849708 : zTin, zTsn, & ! LCOV_EXCL_LINE
1853 3618924854 : phi)
1854 :
1855 : logical, intent(in) :: &
1856 : lsnow ! snow presence: T: has snow, F: no snow
1857 :
1858 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
1859 : zqin, & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1860 : zqsn ! snow layer enthalpy (J m-3)
1861 :
1862 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1863 : zTin, & ! ice layer temperature (C) ! LCOV_EXCL_LINE
1864 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
1865 : zTsn ! snow layer temperature (C)
1866 :
1867 : integer :: &
1868 : k ! vertical layer index
1869 :
1870 : character(len=*),parameter :: subname='(picard_final)'
1871 :
1872 28951398832 : do k = 1, nilyr
1873 28951398832 : zqin(k) = enthalpy_mush_liquid_fraction(zTin(k), phi(k))
1874 : enddo ! k
1875 :
1876 3618924854 : if (lsnow) then
1877 :
1878 6452515890 : do k = 1, nslyr
1879 6452515890 : zqsn(k) = icepack_enthalpy_snow(zTsn(k))
1880 : enddo ! k
1881 :
1882 : endif ! lsnow
1883 :
1884 3618924854 : end subroutine picard_final
1885 :
1886 : !=======================================================================
1887 :
1888 1576096427 : subroutine calc_intercell_thickness(lsnow, hilyr, hslyr, dxp)
1889 :
1890 : logical, intent(in) :: &
1891 : lsnow ! snow presence: T: has snow, F: no snow
1892 :
1893 : real(kind=dbl_kind), intent(in) :: &
1894 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1895 : hslyr ! snow layer thickness (m)
1896 :
1897 : real(kind=dbl_kind), dimension(:), intent(out) :: &
1898 : dxp ! distances between grid points (m)
1899 :
1900 : integer :: &
1901 : l ! vertical index
1902 :
1903 : character(len=*),parameter :: subname='(calc_intercell_thickness)'
1904 :
1905 1576096427 : if (lsnow) then
1906 :
1907 1324471115 : dxp(1) = hslyr / c2
1908 :
1909 1481164635 : do l = 2, nslyr
1910 :
1911 207631330 : dxp(l) = hslyr
1912 :
1913 : enddo ! l
1914 :
1915 1324471115 : dxp(nslyr+1) = (hilyr + hslyr) / c2
1916 :
1917 9271297805 : do l = nslyr+2, nilyr+nslyr
1918 :
1919 9271297805 : dxp(l) = hilyr
1920 :
1921 : enddo ! l
1922 :
1923 1324471115 : dxp(nilyr+nslyr+1) = hilyr / c2
1924 :
1925 : else ! lsnow
1926 :
1927 251625312 : dxp(1) = hilyr / c2
1928 :
1929 1761377184 : do l = 2, nilyr
1930 :
1931 1761377184 : dxp(l) = hilyr
1932 :
1933 : enddo ! l
1934 :
1935 251625312 : dxp(nilyr+1) = hilyr / c2
1936 :
1937 544035590 : do l = nilyr+2, nilyr+nslyr+1
1938 :
1939 544035590 : dxp(l) = c0
1940 :
1941 : enddo ! l
1942 :
1943 : endif ! lsnow
1944 :
1945 1576096427 : end subroutine calc_intercell_thickness
1946 :
1947 : !=======================================================================
1948 :
1949 1576096427 : subroutine calc_intercell_conductivity(lsnow, &
1950 1576096427 : km, ks, & ! LCOV_EXCL_LINE
1951 : hilyr, hslyr, & ! LCOV_EXCL_LINE
1952 1576096427 : kcstar)
1953 :
1954 : logical, intent(in) :: &
1955 : lsnow ! snow presence: T: has snow, F: no snow
1956 :
1957 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1958 : km , & ! ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
1959 : ks ! snow conductivity (W m-1 K-1)
1960 :
1961 : real(kind=dbl_kind), intent(in) :: &
1962 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
1963 : hslyr ! snow layer thickness (m)
1964 :
1965 : real(kind=dbl_kind), dimension(:), intent(out) :: &
1966 : kcstar ! interface conductivities (W m-1 K-1)
1967 :
1968 : real(kind=dbl_kind) :: &
1969 : fe ! distance fraction at interface
1970 :
1971 : integer :: &
1972 : k, & ! vertical layer index ! LCOV_EXCL_LINE
1973 : l ! vertical index
1974 :
1975 : character(len=*),parameter :: subname='(calc_intercell_conductivity)'
1976 :
1977 1576096427 : if (lsnow) then
1978 :
1979 1324471115 : kcstar(1) = ks(1)
1980 :
1981 1481164635 : do l = 2, nslyr
1982 :
1983 156693520 : k = l
1984 207631330 : kcstar(l) = (c2 * ks(k) * ks(k-1)) / (ks(k) + ks(k-1))
1985 :
1986 : enddo ! l
1987 :
1988 1324471115 : fe = hilyr / (hilyr + hslyr)
1989 1324471115 : kcstar(nslyr+1) = c1 / ((c1 - fe) / ks(nslyr) + fe / km(1))
1990 :
1991 9271297805 : do k = 2, nilyr
1992 :
1993 7946826690 : l = k + nslyr
1994 9271297805 : kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
1995 :
1996 : enddo ! k
1997 :
1998 1324471115 : kcstar(nilyr+nslyr+1) = km(nilyr)
1999 :
2000 : else ! lsnow
2001 :
2002 251625312 : kcstar(1) = km(1)
2003 :
2004 1761377184 : do k = 2, nilyr
2005 :
2006 1509751872 : l = k
2007 1761377184 : kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
2008 :
2009 : enddo ! k
2010 :
2011 251625312 : kcstar(nilyr+1) = km(nilyr)
2012 :
2013 544035590 : do l = nilyr+2, nilyr+nslyr+1
2014 :
2015 544035590 : kcstar(l) = c0
2016 :
2017 : enddo ! l
2018 :
2019 : endif ! lsnow
2020 :
2021 1576096427 : end subroutine calc_intercell_conductivity
2022 :
2023 : !=======================================================================
2024 :
2025 3618924854 : subroutine solve_heat_conduction(lsnow, lcold, &
2026 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2027 3618924854 : zqin0, zqsn0, & ! LCOV_EXCL_LINE
2028 3618924854 : phi, dt, & ! LCOV_EXCL_LINE
2029 : qpond, qocn, & ! LCOV_EXCL_LINE
2030 3618924854 : q, w, & ! LCOV_EXCL_LINE
2031 : hilyr, hslyr, & ! LCOV_EXCL_LINE
2032 3618924854 : dxp, kcstar, & ! LCOV_EXCL_LINE
2033 3618924854 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
2034 : fsurfn, dfsurfn_dTsf, & ! LCOV_EXCL_LINE
2035 3618924854 : zTin, zTsn)
2036 :
2037 : logical, intent(in) :: &
2038 : lsnow , & ! snow presence: T: has snow, F: no snow ! LCOV_EXCL_LINE
2039 : lcold ! surface cold: T: surface is cold, F: surface is melting
2040 :
2041 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2042 : zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep ! LCOV_EXCL_LINE
2043 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
2044 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
2045 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
2046 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
2047 :
2048 : real(kind=dbl_kind), intent(inout) :: &
2049 : Tsf ! snow surface temperature (C)
2050 :
2051 : real(kind=dbl_kind), intent(in) :: &
2052 : dt , & ! timestep (s) ! LCOV_EXCL_LINE
2053 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2054 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
2055 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
2056 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2057 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2058 : w , & ! vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
2059 : fsurfn , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
2060 : dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
2061 :
2062 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
2063 : q ! upward interface vertical Darcy flow (m s-1)
2064 :
2065 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2066 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
2067 : kcstar ! interface conductivities (W m-1 K-1)
2068 :
2069 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
2070 : zTin ! ice layer temperature (C)
2071 :
2072 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
2073 : zTsn ! snow layer temperature (C)
2074 :
2075 : real(kind=dbl_kind), dimension(nilyr+nslyr+1) :: &
2076 7237849708 : Ap , & ! diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2077 7237849708 : As , & ! lower off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2078 3618924854 : An , & ! upper off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2079 7237849708 : b , & ! right hand side of matrix solve ! LCOV_EXCL_LINE
2080 3618924854 : T ! ice and snow temperatures
2081 :
2082 : integer :: &
2083 : nyn ! matrix size
2084 :
2085 : character(len=*),parameter :: subname='(solve_heat_conduction)'
2086 :
2087 : ! set up matrix and right hand side - snow
2088 3618924854 : if (lsnow) then
2089 :
2090 3048334638 : if (lcold) then
2091 :
2092 : call matrix_elements_snow_cold(Ap, As, An, b, nyn, &
2093 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2094 : zqin0, zqsn0, & ! LCOV_EXCL_LINE
2095 : qpond, qocn, & ! LCOV_EXCL_LINE
2096 : phi, q, & ! LCOV_EXCL_LINE
2097 : w, & ! LCOV_EXCL_LINE
2098 : hilyr, hslyr, & ! LCOV_EXCL_LINE
2099 : dxp, kcstar, & ! LCOV_EXCL_LINE
2100 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
2101 : fsurfn, dfsurfn_dTsf, & ! LCOV_EXCL_LINE
2102 2720047707 : dt)
2103 2720047707 : if (icepack_warnings_aborted(subname)) return
2104 :
2105 : else ! lcold
2106 :
2107 : call matrix_elements_snow_melt(Ap, As, An, b, nyn, &
2108 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2109 : zqin0, zqsn0, & ! LCOV_EXCL_LINE
2110 : qpond, qocn, & ! LCOV_EXCL_LINE
2111 : phi, q, & ! LCOV_EXCL_LINE
2112 : w, & ! LCOV_EXCL_LINE
2113 : hilyr, hslyr, & ! LCOV_EXCL_LINE
2114 : dxp, kcstar, & ! LCOV_EXCL_LINE
2115 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
2116 328286931 : dt)
2117 328286931 : if (icepack_warnings_aborted(subname)) return
2118 :
2119 : endif ! lcold
2120 :
2121 : else ! lsnow
2122 :
2123 570590216 : if (lcold) then
2124 :
2125 : call matrix_elements_nosnow_cold(Ap, As, An, b, nyn, &
2126 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2127 : zqin0, & ! LCOV_EXCL_LINE
2128 : qpond, qocn, & ! LCOV_EXCL_LINE
2129 : phi, q, & ! LCOV_EXCL_LINE
2130 : w, & ! LCOV_EXCL_LINE
2131 : hilyr, & ! LCOV_EXCL_LINE
2132 : dxp, kcstar, & ! LCOV_EXCL_LINE
2133 : Iswabs, & ! LCOV_EXCL_LINE
2134 : fsurfn, dfsurfn_dTsf, & ! LCOV_EXCL_LINE
2135 364863343 : dt)
2136 364863343 : if (icepack_warnings_aborted(subname)) return
2137 :
2138 : else ! lcold
2139 :
2140 : call matrix_elements_nosnow_melt(Ap, As, An, b, nyn, &
2141 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2142 : zqin0, & ! LCOV_EXCL_LINE
2143 : qpond, qocn, & ! LCOV_EXCL_LINE
2144 : phi, q, & ! LCOV_EXCL_LINE
2145 : w, & ! LCOV_EXCL_LINE
2146 : hilyr, & ! LCOV_EXCL_LINE
2147 : dxp, kcstar, & ! LCOV_EXCL_LINE
2148 : Iswabs, & ! LCOV_EXCL_LINE
2149 205726873 : dt)
2150 205726873 : if (icepack_warnings_aborted(subname)) return
2151 :
2152 : endif ! lcold
2153 :
2154 : endif ! lsnow
2155 :
2156 : ! tridiag to get new temperatures
2157 : call tdma_solve_sparse( &
2158 3618924854 : An(1:nyn), Ap(1:nyn), As(1:nyn), b(1:nyn), T(1:nyn), nyn)
2159 3618924854 : if (icepack_warnings_aborted(subname)) return
2160 :
2161 : call update_temperatures(lsnow, lcold, &
2162 : T, Tsf, & ! LCOV_EXCL_LINE
2163 3618924854 : zTin, zTsn)
2164 3618924854 : if (icepack_warnings_aborted(subname)) return
2165 :
2166 : end subroutine solve_heat_conduction
2167 :
2168 : !=======================================================================
2169 :
2170 3618924854 : subroutine update_temperatures(lsnow, lcold, &
2171 3618924854 : T, Tsf, & ! LCOV_EXCL_LINE
2172 3618924854 : zTin, zTsn)
2173 :
2174 : logical, intent(in) :: &
2175 : lsnow , & ! snow presence: T: has snow, F: no snow ! LCOV_EXCL_LINE
2176 : lcold ! surface cold: T: surface is cold, F: surface is melting
2177 :
2178 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2179 : T ! matrix solution vector
2180 :
2181 : real(kind=dbl_kind), intent(inout) :: &
2182 : Tsf ! snow surface temperature (C)
2183 :
2184 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
2185 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
2186 : zTsn ! snow layer temperature (C)
2187 :
2188 : integer :: &
2189 : l , & ! vertical index ! LCOV_EXCL_LINE
2190 : k ! vertical layer index
2191 :
2192 : character(len=*),parameter :: subname='(update_temperatures)'
2193 :
2194 3618924854 : if (lsnow) then
2195 :
2196 3048334638 : if (lcold) then
2197 :
2198 2720047707 : Tsf = T(1)
2199 :
2200 5772644940 : do k = 1, nslyr
2201 3052597233 : l = k + 1
2202 5772644940 : zTsn(k) = T(l)
2203 : enddo ! k
2204 :
2205 21760381656 : do k = 1, nilyr
2206 19040333949 : l = k + nslyr + 1
2207 21760381656 : zTin(k) = T(l)
2208 : enddo ! k
2209 :
2210 : else ! lcold
2211 :
2212 679870950 : do k = 1, nslyr
2213 351584019 : l = k
2214 679870950 : zTsn(k) = T(l)
2215 : enddo ! k
2216 :
2217 2626295448 : do k = 1, nilyr
2218 2298008517 : l = k + nslyr
2219 2626295448 : zTin(k) = T(l)
2220 : enddo ! k
2221 :
2222 : endif ! lcold
2223 :
2224 : else ! lsnow
2225 :
2226 570590216 : if (lcold) then
2227 :
2228 364863343 : Tsf = T(1)
2229 :
2230 2918906744 : do k = 1, nilyr
2231 2554043401 : l = k + 1
2232 2918906744 : zTin(k) = T(l)
2233 : enddo ! k
2234 :
2235 : else ! lcold
2236 :
2237 1645814984 : do k = 1, nilyr
2238 1440088111 : l = k
2239 1645814984 : zTin(k) = T(l)
2240 : enddo ! k
2241 :
2242 : endif ! lcold
2243 :
2244 : endif ! lsnow
2245 :
2246 3618924854 : end subroutine update_temperatures
2247 :
2248 : !=======================================================================
2249 :
2250 411453746 : subroutine matrix_elements_nosnow_melt(Ap, As, An, b, nyn, &
2251 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2252 205726873 : zqin0, & ! LCOV_EXCL_LINE
2253 : qpond, qocn, & ! LCOV_EXCL_LINE
2254 205726873 : phi, q, & ! LCOV_EXCL_LINE
2255 : w, & ! LCOV_EXCL_LINE
2256 : hilyr, & ! LCOV_EXCL_LINE
2257 205726873 : dxp, kcstar, & ! LCOV_EXCL_LINE
2258 205726873 : Iswabs, & ! LCOV_EXCL_LINE
2259 : dt)
2260 :
2261 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2262 : Ap , & ! diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2263 : As , & ! lower off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2264 : An , & ! upper off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2265 : b ! right hand side of matrix solve
2266 :
2267 : integer, intent(out) :: &
2268 : nyn ! matrix size
2269 :
2270 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2271 : zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep ! LCOV_EXCL_LINE
2272 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
2273 : phi ! ice layer liquid fraction
2274 :
2275 : real(kind=dbl_kind), intent(in) :: &
2276 : Tsf , & ! snow surface temperature (C) ! LCOV_EXCL_LINE
2277 : dt , & ! timestep (s) ! LCOV_EXCL_LINE
2278 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2279 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
2280 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2281 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2282 : w ! downwards vertical flushing Darcy velocity (m/s)
2283 :
2284 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
2285 : q ! upward interface vertical Darcy flow (m s-1)
2286 :
2287 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2288 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
2289 : kcstar ! interface conductivities (W m-1 K-1)
2290 :
2291 : integer :: &
2292 : k , & ! vertical layer index ! LCOV_EXCL_LINE
2293 : l ! vertical index
2294 :
2295 : character(len=*),parameter :: subname='(matrix_elements_nosnow_melt)'
2296 :
2297 : ! surface layer
2298 205726873 : k = 1
2299 205726873 : l = k
2300 :
2301 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2302 : kcstar(k+1) / dxp(k+1) + & ! LCOV_EXCL_LINE
2303 : kcstar(k) / dxp(k) + & ! LCOV_EXCL_LINE
2304 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2305 205726873 : w * cp_ocn * rhow
2306 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2307 205726873 : q(k) * cp_ocn * rhow
2308 205726873 : An(l) = c0
2309 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2310 : (kcstar(k) / dxp(k)) * Tsf + & ! LCOV_EXCL_LINE
2311 205726873 : w * qpond
2312 :
2313 : ! interior ice layers
2314 1234361238 : do k = 2, nilyr-1
2315 :
2316 1028634365 : l = k
2317 :
2318 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2319 : kcstar(k+1) / dxp(k+1) + & ! LCOV_EXCL_LINE
2320 : kcstar(k) / dxp(k) + & ! LCOV_EXCL_LINE
2321 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2322 1028634365 : w * cp_ocn * rhow
2323 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2324 1028634365 : q(k) * cp_ocn * rhow
2325 : An(l) = -kcstar(k) / dxp(k) - &
2326 1028634365 : w * cp_ocn * rhow
2327 1234361238 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2328 :
2329 : enddo ! k
2330 :
2331 : ! bottom layer
2332 205726873 : k = nilyr
2333 205726873 : l = k
2334 :
2335 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2336 : kcstar(k+1) / dxp(k+1) + & ! LCOV_EXCL_LINE
2337 : kcstar(k) / dxp(k) + & ! LCOV_EXCL_LINE
2338 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2339 205726873 : w * cp_ocn * rhow
2340 205726873 : As(l) = c0
2341 : An(l) = -kcstar(k) / dxp(k) - &
2342 205726873 : w * cp_ocn * rhow
2343 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2344 : (kcstar(k+1) * Tbot) / dxp(k+1) + & ! LCOV_EXCL_LINE
2345 205726873 : q(k) * qocn
2346 :
2347 205726873 : nyn = nilyr
2348 :
2349 205726873 : end subroutine matrix_elements_nosnow_melt
2350 :
2351 : !=======================================================================
2352 :
2353 729726686 : subroutine matrix_elements_nosnow_cold(Ap, As, An, b, nyn, &
2354 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2355 0 : zqin0, & ! LCOV_EXCL_LINE
2356 : qpond, qocn, & ! LCOV_EXCL_LINE
2357 364863343 : phi, q, & ! LCOV_EXCL_LINE
2358 : w, & ! LCOV_EXCL_LINE
2359 : hilyr, & ! LCOV_EXCL_LINE
2360 364863343 : dxp, kcstar, & ! LCOV_EXCL_LINE
2361 364863343 : Iswabs, & ! LCOV_EXCL_LINE
2362 : fsurfn, dfsurfn_dTsf, & ! LCOV_EXCL_LINE
2363 : dt)
2364 :
2365 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2366 : Ap , & ! diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2367 : As , & ! lower off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2368 : An , & ! upper off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2369 : b ! right hand side of matrix solve
2370 :
2371 : integer, intent(out) :: &
2372 : nyn ! matrix size
2373 :
2374 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2375 : zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep ! LCOV_EXCL_LINE
2376 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
2377 : phi ! ice layer liquid fraction
2378 :
2379 : real(kind=dbl_kind), intent(in) :: &
2380 : Tsf , & ! snow surface temperature (C) ! LCOV_EXCL_LINE
2381 : dt , & ! timestep (s) ! LCOV_EXCL_LINE
2382 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2383 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
2384 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2385 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2386 : w , & ! downwards vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
2387 : fsurfn , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
2388 : dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
2389 :
2390 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
2391 : q ! upward interface vertical Darcy flow (m s-1)
2392 :
2393 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2394 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
2395 : kcstar ! interface conductivities (W m-1 K-1)
2396 :
2397 : integer :: &
2398 : k , & ! vertical layer index ! LCOV_EXCL_LINE
2399 : l ! vertical index
2400 :
2401 : character(len=*),parameter :: subname='(matrix_elements_nosnow_cold)'
2402 :
2403 : ! surface temperature
2404 364863343 : l = 1
2405 364863343 : Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
2406 364863343 : As(l) = kcstar(1) / dxp(1)
2407 364863343 : An(l) = c0
2408 364863343 : b (l) = dfsurfn_dTsf * Tsf - fsurfn
2409 :
2410 : ! surface layer
2411 364863343 : k = 1
2412 364863343 : l = k + 1
2413 :
2414 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2415 : kcstar(k+1) / dxp(k+1) + & ! LCOV_EXCL_LINE
2416 : kcstar(k) / dxp(k) + & ! LCOV_EXCL_LINE
2417 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2418 364863343 : w * cp_ocn * rhow
2419 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2420 364863343 : q(k) * cp_ocn * rhow
2421 364863343 : An(l) = -kcstar(k) / dxp(k)
2422 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2423 364863343 : w * qpond
2424 :
2425 : ! interior ice layers
2426 2189180058 : do k = 2, nilyr-1
2427 :
2428 1824316715 : l = k + 1
2429 :
2430 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2431 : kcstar(k+1) / dxp(k+1) + & ! LCOV_EXCL_LINE
2432 : kcstar(k) / dxp(k) + & ! LCOV_EXCL_LINE
2433 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2434 1824316715 : w * cp_ocn * rhow
2435 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2436 1824316715 : q(k) * cp_ocn * rhow
2437 : An(l) = -kcstar(k) / dxp(k) - &
2438 1824316715 : w * cp_ocn * rhow
2439 2189180058 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2440 :
2441 : enddo ! k
2442 :
2443 : ! bottom layer
2444 364863343 : k = nilyr
2445 364863343 : l = k + 1
2446 :
2447 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2448 : kcstar(k+1) / dxp(k+1) + & ! LCOV_EXCL_LINE
2449 : kcstar(k) / dxp(k) + & ! LCOV_EXCL_LINE
2450 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2451 364863343 : w * cp_ocn * rhow
2452 364863343 : As(l) = c0
2453 : An(l) = -kcstar(k) / dxp(k) - &
2454 364863343 : w * cp_ocn * rhow
2455 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2456 : (kcstar(k+1) * Tbot) / dxp(k+1) + & ! LCOV_EXCL_LINE
2457 364863343 : q(k) * qocn
2458 :
2459 364863343 : nyn = nilyr + 1
2460 :
2461 364863343 : end subroutine matrix_elements_nosnow_cold
2462 :
2463 : !=======================================================================
2464 :
2465 656573862 : subroutine matrix_elements_snow_melt(Ap, As, An, b, nyn, &
2466 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2467 656573862 : zqin0, zqsn0, & ! LCOV_EXCL_LINE
2468 : qpond, qocn, & ! LCOV_EXCL_LINE
2469 328286931 : phi, q, & ! LCOV_EXCL_LINE
2470 : w, & ! LCOV_EXCL_LINE
2471 : hilyr, hslyr, & ! LCOV_EXCL_LINE
2472 328286931 : dxp, kcstar, & ! LCOV_EXCL_LINE
2473 656573862 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
2474 : dt)
2475 :
2476 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2477 : Ap , & ! diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2478 : As , & ! lower off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2479 : An , & ! upper off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2480 : b ! right hand side of matrix solve
2481 :
2482 : integer, intent(out) :: &
2483 : nyn ! matrix size
2484 :
2485 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2486 : zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep ! LCOV_EXCL_LINE
2487 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
2488 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
2489 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
2490 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
2491 :
2492 : real(kind=dbl_kind), intent(in) :: &
2493 : Tsf , & ! snow surface temperature (C) ! LCOV_EXCL_LINE
2494 : dt , & ! timestep (s) ! LCOV_EXCL_LINE
2495 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2496 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
2497 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
2498 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2499 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2500 : w ! downwards vertical flushing Darcy velocity (m/s)
2501 :
2502 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
2503 : q ! upward interface vertical Darcy flow (m s-1)
2504 :
2505 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2506 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
2507 : kcstar ! interface conductivities (W m-1 K-1)
2508 :
2509 : integer :: &
2510 : k , & ! vertical layer index ! LCOV_EXCL_LINE
2511 : l ! vertical index
2512 :
2513 : character(len=*),parameter :: subname='(matrix_elements_snow_melt)'
2514 :
2515 : ! surface layer
2516 328286931 : k = 1
2517 328286931 : l = k
2518 :
2519 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2520 : kcstar(l+1) / dxp(l+1) + & ! LCOV_EXCL_LINE
2521 328286931 : kcstar(l) / dxp(l)
2522 328286931 : As(l) = -kcstar(l+1) / dxp(l+1)
2523 328286931 : An(l) = c0
2524 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k) + &
2525 328286931 : (kcstar(l) * Tsf) / dxp(l)
2526 :
2527 : ! interior snow layers
2528 351584019 : do k = 2, nslyr
2529 :
2530 23297088 : l = k
2531 :
2532 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2533 : kcstar(l+1) / dxp(l+1) + & ! LCOV_EXCL_LINE
2534 23297088 : kcstar(l) / dxp(l)
2535 23297088 : As(l) = -kcstar(l+1) / dxp(l+1)
2536 23297088 : An(l) = -kcstar(l) / dxp(l)
2537 32580284 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2538 :
2539 : enddo ! k
2540 :
2541 : ! top ice layer
2542 328286931 : k = 1
2543 328286931 : l = nslyr + k
2544 :
2545 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2546 : kcstar(l+1) / dxp(l+1) + & ! LCOV_EXCL_LINE
2547 : kcstar(l) / dxp(l) + & ! LCOV_EXCL_LINE
2548 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2549 328286931 : w * cp_ocn * rhow
2550 : As(l) = -kcstar(l+1) / dxp(l+1) - &
2551 328286931 : q(k) * cp_ocn * rhow
2552 328286931 : An(l) = -kcstar(l) / dxp(l)
2553 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2554 328286931 : w * qpond
2555 :
2556 : ! interior ice layers
2557 1969721586 : do k = 2, nilyr-1
2558 :
2559 1641434655 : l = nslyr + k
2560 :
2561 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2562 : kcstar(l+1) / dxp(l+1) + & ! LCOV_EXCL_LINE
2563 : kcstar(l) / dxp(l) + & ! LCOV_EXCL_LINE
2564 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2565 1641434655 : w * cp_ocn * rhow
2566 : As(l) = -kcstar(l+1) / dxp(l+1) - &
2567 1641434655 : q(k) * cp_ocn * rhow
2568 : An(l) = -kcstar(l) / dxp(l) - &
2569 1641434655 : w * cp_ocn * rhow
2570 1969721586 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2571 :
2572 : enddo ! k
2573 :
2574 : ! bottom layer
2575 328286931 : k = nilyr
2576 328286931 : l = nilyr + nslyr
2577 :
2578 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2579 : kcstar(l+1) / dxp(l+1) + & ! LCOV_EXCL_LINE
2580 : kcstar(l) / dxp(l) + & ! LCOV_EXCL_LINE
2581 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2582 328286931 : w * cp_ocn * rhow
2583 328286931 : As(l) = c0
2584 : An(l) = -kcstar(l) / dxp(l) - &
2585 328286931 : w * cp_ocn * rhow
2586 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2587 : (kcstar(l+1) * Tbot) / dxp(l+1) + & ! LCOV_EXCL_LINE
2588 328286931 : q(k) * qocn
2589 :
2590 328286931 : nyn = nilyr + nslyr
2591 :
2592 328286931 : end subroutine matrix_elements_snow_melt
2593 :
2594 : !=======================================================================
2595 :
2596 5440095414 : subroutine matrix_elements_snow_cold(Ap, As, An, b, nyn, &
2597 : Tsf, Tbot, & ! LCOV_EXCL_LINE
2598 5440095414 : zqin0, zqsn0, & ! LCOV_EXCL_LINE
2599 : qpond, qocn, & ! LCOV_EXCL_LINE
2600 2720047707 : phi, q, & ! LCOV_EXCL_LINE
2601 : w, & ! LCOV_EXCL_LINE
2602 : hilyr, hslyr, & ! LCOV_EXCL_LINE
2603 2720047707 : dxp, kcstar, & ! LCOV_EXCL_LINE
2604 5440095414 : Iswabs, Sswabs, & ! LCOV_EXCL_LINE
2605 : fsurfn, dfsurfn_dTsf, & ! LCOV_EXCL_LINE
2606 : dt)
2607 :
2608 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2609 : Ap , & ! diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2610 : As , & ! lower off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2611 : An , & ! upper off-diagonal of tridiagonal matrix ! LCOV_EXCL_LINE
2612 : b ! right hand side of matrix solve
2613 :
2614 : integer, intent(out) :: &
2615 : nyn ! matrix size
2616 :
2617 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2618 : zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep ! LCOV_EXCL_LINE
2619 : Iswabs , & ! SW radiation absorbed in ice layers (W m-2) ! LCOV_EXCL_LINE
2620 : phi , & ! ice layer liquid fraction ! LCOV_EXCL_LINE
2621 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep ! LCOV_EXCL_LINE
2622 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
2623 :
2624 : real(kind=dbl_kind), intent(in) :: &
2625 : Tsf , & ! snow surface temperature (C) ! LCOV_EXCL_LINE
2626 : dt , & ! timestep (s) ! LCOV_EXCL_LINE
2627 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2628 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
2629 : Tbot , & ! ice bottom surfce temperature (deg C) ! LCOV_EXCL_LINE
2630 : qpond , & ! melt pond brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2631 : qocn , & ! ocean brine enthalpy (J m-3) ! LCOV_EXCL_LINE
2632 : w , & ! downwards vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
2633 : fsurfn , & ! net flux to top surface, excluding fcondtop ! LCOV_EXCL_LINE
2634 : dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
2635 :
2636 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
2637 : q ! upward interface vertical Darcy flow (m s-1)
2638 :
2639 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2640 : dxp , & ! distances between grid points (m) ! LCOV_EXCL_LINE
2641 : kcstar ! interface conductivities (W m-1 K-1)
2642 :
2643 : integer :: &
2644 : k , & ! vertical layer index ! LCOV_EXCL_LINE
2645 : l , & ! matrix index ! LCOV_EXCL_LINE
2646 : m ! vertical index
2647 :
2648 : character(len=*),parameter :: subname='(matrix_elements_snow_cold)'
2649 :
2650 : ! surface temperature
2651 2720047707 : l = 1
2652 2720047707 : Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
2653 2720047707 : As(l) = kcstar(1) / dxp(1)
2654 2720047707 : An(l) = c0
2655 2720047707 : b (l) = dfsurfn_dTsf * Tsf - fsurfn
2656 :
2657 : ! surface layer
2658 2720047707 : k = 1
2659 2720047707 : l = k + 1
2660 2720047707 : m = 1
2661 :
2662 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2663 : kcstar(m+1) / dxp(m+1) + & ! LCOV_EXCL_LINE
2664 2720047707 : kcstar(m) / dxp(m)
2665 2720047707 : As(l) = -kcstar(m+1) / dxp(m+1)
2666 2720047707 : An(l) = -kcstar(m) / dxp(m)
2667 2720047707 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2668 :
2669 : ! interior snow layers
2670 3052597233 : do k = 2, nslyr
2671 :
2672 332549526 : l = k + 1
2673 332549526 : m = k
2674 :
2675 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2676 : kcstar(m+1) / dxp(m+1) + & ! LCOV_EXCL_LINE
2677 332549526 : kcstar(m) / dxp(m)
2678 332549526 : As(l) = -kcstar(m+1) / dxp(m+1)
2679 332549526 : An(l) = -kcstar(m) / dxp(m)
2680 439430205 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2681 :
2682 : enddo ! k
2683 :
2684 : ! top ice layer
2685 2720047707 : k = 1
2686 2720047707 : l = nslyr + k + 1
2687 2720047707 : m = k + nslyr
2688 :
2689 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2690 : kcstar(m+1) / dxp(m+1) + & ! LCOV_EXCL_LINE
2691 : kcstar(m) / dxp(m) + & ! LCOV_EXCL_LINE
2692 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2693 2720047707 : w * cp_ocn * rhow
2694 : As(l) = -kcstar(m+1) / dxp(m+1) - &
2695 2720047707 : q(k) * cp_ocn * rhow
2696 2720047707 : An(l) = -kcstar(m) / dxp(m)
2697 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2698 2720047707 : w * qpond
2699 :
2700 : ! interior ice layers
2701 16320286242 : do k = 2, nilyr-1
2702 :
2703 13600238535 : l = nslyr + k + 1
2704 13600238535 : m = k + nslyr
2705 :
2706 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2707 : kcstar(m+1) / dxp(m+1) + & ! LCOV_EXCL_LINE
2708 : kcstar(m) / dxp(m) + & ! LCOV_EXCL_LINE
2709 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2710 13600238535 : w * cp_ocn * rhow
2711 : As(l) = -kcstar(m+1) / dxp(m+1) - &
2712 13600238535 : q(k) * cp_ocn * rhow
2713 : An(l) = -kcstar(m) / dxp(m) - &
2714 13600238535 : w * cp_ocn * rhow
2715 16320286242 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2716 :
2717 : enddo ! k
2718 :
2719 : ! bottom layer
2720 2720047707 : k = nilyr
2721 2720047707 : l = nilyr + nslyr + 1
2722 2720047707 : m = k + nslyr
2723 :
2724 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2725 : kcstar(m+1) / dxp(m+1) + & ! LCOV_EXCL_LINE
2726 : kcstar(m) / dxp(m) + & ! LCOV_EXCL_LINE
2727 : q(k) * cp_ocn * rhow + & ! LCOV_EXCL_LINE
2728 2720047707 : w * cp_ocn * rhow
2729 2720047707 : As(l) = c0
2730 : An(l) = -kcstar(m) / dxp(m) - &
2731 2720047707 : w * cp_ocn * rhow
2732 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2733 : (kcstar(m+1) * Tbot) / dxp(m+1) + & ! LCOV_EXCL_LINE
2734 2720047707 : q(k) * qocn
2735 :
2736 2720047707 : nyn = nilyr + nslyr + 1
2737 :
2738 2720047707 : end subroutine matrix_elements_snow_cold
2739 :
2740 : !=======================================================================
2741 :
2742 3152192854 : subroutine solve_salinity(zSin, Sbr, &
2743 : Spond, sss, & ! LCOV_EXCL_LINE
2744 3152192854 : q, dSdt, & ! LCOV_EXCL_LINE
2745 : w, hilyr, & ! LCOV_EXCL_LINE
2746 : dt)
2747 :
2748 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
2749 : zSin ! ice layer bulk salinity (ppt)
2750 :
2751 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2752 : Sbr , & ! ice layer brine salinity (ppt) ! LCOV_EXCL_LINE
2753 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
2754 :
2755 : real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
2756 : q ! upward interface vertical Darcy flow (m s-1)
2757 :
2758 : real(kind=dbl_kind), intent(in) :: &
2759 : Spond , & ! melt pond salinity (ppt) ! LCOV_EXCL_LINE
2760 : sss , & ! sea surface salinity (ppt) ! LCOV_EXCL_LINE
2761 : w , & ! vertical flushing Darcy velocity (m/s) ! LCOV_EXCL_LINE
2762 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2763 : dt ! timestep (s)
2764 :
2765 : integer :: &
2766 : k ! vertical layer index
2767 :
2768 : real(kind=dbl_kind), parameter :: &
2769 : S_min = p01
2770 :
2771 : real(kind=dbl_kind), dimension(nilyr) :: &
2772 1576096427 : zSin0
2773 :
2774 : character(len=*),parameter :: subname='(solve_salinity)'
2775 :
2776 12608771416 : zSin0 = zSin
2777 :
2778 1576096427 : k = 1
2779 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2780 : ((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + & ! LCOV_EXCL_LINE
2781 : dSdt(k) + & ! LCOV_EXCL_LINE
2782 1576096427 : (w * (Spond - Sbr(k))) / hilyr) * dt)
2783 :
2784 9456578562 : do k = 2, nilyr-1
2785 :
2786 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2787 : ((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + & ! LCOV_EXCL_LINE
2788 : dSdt(k) + & ! LCOV_EXCL_LINE
2789 9456578562 : (w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
2790 :
2791 : enddo ! k
2792 :
2793 1576096427 : k = nilyr
2794 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2795 : ((q(k) * (sss - Sbr(k))) / hilyr + & ! LCOV_EXCL_LINE
2796 : dSdt(k) + & ! LCOV_EXCL_LINE
2797 1576096427 : (w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
2798 :
2799 :
2800 12608771416 : if (minval(zSin) < c0) then
2801 :
2802 :
2803 0 : write(warnstr,*) subname, (q(k) * (Sbr(k+1) - Sbr(k))) / hilyr, &
2804 0 : dSdt(k) , & ! LCOV_EXCL_LINE
2805 0 : (w * (Spond - Sbr(k))) / hilyr
2806 0 : call icepack_warnings_add(warnstr)
2807 :
2808 0 : do k = 1, nilyr
2809 :
2810 0 : write(warnstr,*) subname, k, zSin(k), zSin0(k)
2811 0 : call icepack_warnings_add(warnstr)
2812 :
2813 : enddo
2814 :
2815 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
2816 0 : return
2817 :
2818 : endif
2819 :
2820 : end subroutine solve_salinity
2821 :
2822 : !=======================================================================
2823 :
2824 3618924854 : subroutine tdma_solve_sparse(a, b, c, d, x, n)
2825 :
2826 : ! perform a tri-diagonal solve with TDMA using a sparse tridiagoinal matrix
2827 :
2828 : integer(kind=int_kind), intent(in) :: &
2829 : n ! matrix size
2830 :
2831 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2832 : a , & ! matrix lower off-diagonal ! LCOV_EXCL_LINE
2833 : b , & ! matrix diagonal ! LCOV_EXCL_LINE
2834 : c , & ! matrix upper off-diagonal ! LCOV_EXCL_LINE
2835 : d ! right hand side vector
2836 :
2837 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2838 : x ! solution vector
2839 :
2840 : real(kind=dbl_kind), dimension(nilyr+nslyr+1) :: &
2841 7237849708 : cp , & ! modified upper off-diagonal vector ! LCOV_EXCL_LINE
2842 7237849708 : dp ! modified right hand side vector
2843 :
2844 : integer(kind=int_kind) :: &
2845 : i ! vector index
2846 :
2847 : character(len=*),parameter :: subname='(tdma_solve_sparse)'
2848 :
2849 : ! forward sweep
2850 3618924854 : cp(1) = c(1) / b(1)
2851 28202641426 : do i = 2, n-1
2852 28202641426 : cp(i) = c(i) / (b(i) - cp(i-1)*a(i))
2853 : enddo
2854 :
2855 3618924854 : dp(1) = d(1) / b(1)
2856 31821566280 : do i = 2, n
2857 31821566280 : dp(i) = (d(i) - dp(i-1)*a(i)) / (b(i) - cp(i-1)*a(i))
2858 : enddo
2859 :
2860 : ! back substitution
2861 3618924854 : x(n) = dp(n)
2862 31821566280 : do i = n-1,1,-1
2863 31821566280 : x(i) = dp(i) - cp(i)*x(i+1)
2864 : enddo
2865 :
2866 3618924854 : end subroutine tdma_solve_sparse
2867 :
2868 : !=======================================================================
2869 : ! Effect of salinity
2870 : !=======================================================================
2871 :
2872 19567540244 : function permeability(phi) result(perm)
2873 :
2874 : ! given the liquid fraction calculate the permeability
2875 : ! See Golden et al. 2007
2876 :
2877 : real(kind=dbl_kind), intent(in) :: &
2878 : phi ! liquid fraction
2879 :
2880 : real(kind=dbl_kind) :: &
2881 : perm ! permeability (m2)
2882 :
2883 : real(kind=dbl_kind), parameter :: &
2884 : phic = p05 ! critical liquid fraction for impermeability
2885 :
2886 : character(len=*),parameter :: subname='(permeability)'
2887 :
2888 19567540244 : perm = 3.0e-8_dbl_kind * max(phi - phic, c0)**3
2889 :
2890 19567540244 : end function permeability
2891 :
2892 : !=======================================================================
2893 :
2894 0 : subroutine explicit_flow_velocities(zSin, &
2895 1408028216 : zTin, Tsf, & ! LCOV_EXCL_LINE
2896 1408028216 : Tbot, q, & ! LCOV_EXCL_LINE
2897 2816056432 : dSdt, Sbr, & ! LCOV_EXCL_LINE
2898 1408028216 : qbr, dt, & ! LCOV_EXCL_LINE
2899 : sss, qocn, & ! LCOV_EXCL_LINE
2900 : hilyr, hin)
2901 :
2902 : ! calculate the rapid gravity drainage mode Darcy velocity and the
2903 : ! slow mode drainage rate
2904 :
2905 : real(kind=dbl_kind), dimension(:), intent(in) :: &
2906 : zSin, & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
2907 : zTin ! ice layer temperature (C)
2908 :
2909 : real(kind=dbl_kind), intent(in) :: &
2910 : Tsf , & ! ice/snow surface temperature (C) ! LCOV_EXCL_LINE
2911 : Tbot , & ! ice bottom temperature (C) ! LCOV_EXCL_LINE
2912 : dt , & ! time step (s) ! LCOV_EXCL_LINE
2913 : sss , & ! sea surface salinty (ppt) ! LCOV_EXCL_LINE
2914 : qocn , & ! ocean enthalpy (J m-3) ! LCOV_EXCL_LINE
2915 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
2916 : hin ! ice thickness (m)
2917 :
2918 : real(kind=dbl_kind), dimension(0:nilyr), intent(out) :: &
2919 : q ! rapid mode upward interface vertical Darcy flow (m s-1)
2920 :
2921 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2922 : dSdt ! slow mode drainage rate (ppt s-1)
2923 :
2924 : real(kind=dbl_kind), dimension(:), intent(out) :: &
2925 : Sbr , & ! ice layer brine salinity (ppt) ! LCOV_EXCL_LINE
2926 : qbr ! ice layer brine enthalpy (J m-3)
2927 :
2928 : real(kind=dbl_kind), parameter :: &
2929 : kappal = 8.824e-8_dbl_kind, & ! heat diffusivity of liquid ! LCOV_EXCL_LINE
2930 : fracmax = p2 , & ! limiting advective layer fraction ! LCOV_EXCL_LINE
2931 : zSin_min = p1 , & ! minimum bulk salinity (ppt) ! LCOV_EXCL_LINE
2932 : safety_factor = c10 ! to prevent negative salinities
2933 :
2934 : real(kind=dbl_kind), dimension(1:nilyr) :: &
2935 2816056432 : phi ! ice layer liquid fraction
2936 :
2937 : real(kind=dbl_kind), dimension(0:nilyr) :: &
2938 2816056432 : rho ! ice layer brine density (kg m-3)
2939 :
2940 : real(kind=dbl_kind) :: &
2941 : rho_ocn , & ! ocean density (kg m-3) ! LCOV_EXCL_LINE
2942 : perm_min , & ! minimum permeability from layer to ocean (m2) ! LCOV_EXCL_LINE
2943 : perm_harm , & ! harmonic mean of permeability from layer to ocean (m2) ! LCOV_EXCL_LINE
2944 : rho_sum , & ! sum of the brine densities from layer to ocean (kg m-3) ! LCOV_EXCL_LINE
2945 : rho_pipe , & ! density of the brine in the channel (kg m-3) ! LCOV_EXCL_LINE
2946 : z , & ! distance to layer from top surface (m) ! LCOV_EXCL_LINE
2947 : perm , & ! ice layer permeability (m2) ! LCOV_EXCL_LINE
2948 : drho , & ! brine density difference between layer and ocean (kg m-3) ! LCOV_EXCL_LINE
2949 : Ra , & ! local mush Rayleigh number ! LCOV_EXCL_LINE
2950 : rn , & ! real value of number of layers considered ! LCOV_EXCL_LINE
2951 : L , & ! thickness of the layers considered (m) ! LCOV_EXCL_LINE
2952 : dx , & ! horizontal size of convective flow (m) ! LCOV_EXCL_LINE
2953 : dx2 , & ! square of the horizontal size of convective flow (m2) ! LCOV_EXCL_LINE
2954 : Am , & ! A parameter for mush ! LCOV_EXCL_LINE
2955 : Bm , & ! B parameter for mush ! LCOV_EXCL_LINE
2956 : Ap , & ! A parameter for channel ! LCOV_EXCL_LINE
2957 : Bp , & ! B parameter for channel ! LCOV_EXCL_LINE
2958 : qlimit , & ! limit to vertical Darcy flow for numerical stability ! LCOV_EXCL_LINE
2959 : dS_guess , & ! expected bulk salinity without limits ! LCOV_EXCL_LINE
2960 : alpha , & ! desalination limiting factor ! LCOV_EXCL_LINE
2961 : ra_constants ! for Rayleigh number
2962 :
2963 : integer(kind=int_kind) :: &
2964 : k ! ice layer index
2965 :
2966 : character(len=*),parameter :: subname='(explicit_flow_velocities)'
2967 :
2968 : ! initial downward sweep - determine derived physical quantities
2969 11264225728 : do k = 1, nilyr
2970 :
2971 9856197512 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
2972 9856197512 : phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
2973 9856197512 : qbr(k) = enthalpy_brine(zTin(k))
2974 11264225728 : rho(k) = icepack_mushy_density_brine(Sbr(k))
2975 :
2976 : enddo ! k
2977 :
2978 1408028216 : rho(0) = rho(1)
2979 :
2980 : ! ocean conditions
2981 1408028216 : Sbr(nilyr+1) = sss
2982 1408028216 : qbr(nilyr+1) = qocn
2983 1408028216 : rho_ocn = icepack_mushy_density_brine(sss)
2984 :
2985 : ! initialize accumulated quantities
2986 1408028216 : perm_min = bignum
2987 1408028216 : perm_harm = c0
2988 1408028216 : rho_sum = c0
2989 :
2990 : ! limit to q for numerical stability
2991 1408028216 : qlimit = (fracmax * hilyr) / dt
2992 :
2993 : ! no flow through ice top surface
2994 1408028216 : q(0) = c0
2995 :
2996 1408028216 : ra_constants = gravit / (viscosity_dyn * kappal)
2997 : ! first iterate over layers going up
2998 11264225728 : do k = nilyr, 1, -1
2999 :
3000 : ! vertical position from ice top surface
3001 9856197512 : z = ((real(k, dbl_kind) - p5) / real(nilyr, dbl_kind)) * hin
3002 :
3003 : ! permeabilities
3004 9856197512 : perm = permeability(phi(k))
3005 9856197512 : perm_min = min(perm_min,perm)
3006 9856197512 : perm_harm = perm_harm + (c1 / max(perm,1.0e-30_dbl_kind))
3007 :
3008 : ! densities
3009 9856197512 : rho_sum = rho_sum + rho(k)
3010 : !rho_pipe = rho(k)
3011 9856197512 : rho_pipe = p5 * (rho(k) + rho(k-1))
3012 9856197512 : drho = max(rho(k) - rho_ocn, c0)
3013 :
3014 : ! mush Rayleigh number
3015 9856197512 : Ra = drho * (hin-z) * perm_min * ra_constants
3016 :
3017 : ! height of mush layer to layer k
3018 9856197512 : rn = real(nilyr-k+1,dbl_kind)
3019 9856197512 : L = rn * hilyr
3020 :
3021 : ! horizontal size of convection
3022 9856197512 : dx = L * c2 * aspect_rapid_mode
3023 9856197512 : dx2 = dx**2
3024 :
3025 : ! determine vertical Darcy flow
3026 9856197512 : Am = (dx2 * rn) / (viscosity_dyn * perm_harm)
3027 9856197512 : Bm = (-gravit * rho_sum) / rn
3028 :
3029 9856197512 : Ap = (pi * a_rapid_mode**4) / (c8 * viscosity_dyn)
3030 9856197512 : Bp = -rho_pipe * gravit
3031 :
3032 9856197512 : q(k) = max((Am / dx2) * ((-Ap*Bp - Am*Bm) / (Am + Ap) + Bm), 1.0e-30_dbl_kind)
3033 :
3034 : ! modify by Rayleigh number and advection limit
3035 9856197512 : q(k) = min(q(k) * (max(Ra - Rac_rapid_mode, c0) / (Ra+puny)), qlimit)
3036 :
3037 : ! late stage drainage
3038 : dSdt(k) = dSdt_slow_mode * (max((zSin(k) - phi_c_slow_mode*Sbr(k)), c0) &
3039 9856197512 : * max((Tbot - Tsf), c0)) / (hin + 0.001_dbl_kind)
3040 :
3041 9856197512 : dSdt(k) = max(dSdt(k), (-zSin(k) * 0.5_dbl_kind) / dt)
3042 :
3043 : ! restrict flows to prevent too much salt loss
3044 9856197512 : dS_guess = (((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + dSdt(k)) * dt) * safety_factor
3045 :
3046 9856197512 : if (abs(dS_guess) < puny) then
3047 4623469265 : alpha = c1
3048 : else
3049 5232728247 : alpha = (zSin_min - zSin(k)) / dS_guess
3050 : endif
3051 :
3052 9856197512 : if (alpha < c0 .or. alpha > c1) alpha = c1
3053 :
3054 9856197512 : q(k) = q(k) * alpha
3055 11264225728 : dSdt(k) = dSdt(k) * alpha
3056 :
3057 : enddo ! k
3058 :
3059 1408028216 : end subroutine explicit_flow_velocities
3060 :
3061 : !=======================================================================
3062 : ! Flushing
3063 : !=======================================================================
3064 :
3065 1408028216 : subroutine flushing_velocity(zTin, phi, &
3066 : hin, hsn, & ! LCOV_EXCL_LINE
3067 : hilyr, & ! LCOV_EXCL_LINE
3068 : hpond, apond, & ! LCOV_EXCL_LINE
3069 : dt, w)
3070 :
3071 : ! calculate the vertical flushing Darcy velocity (positive downward)
3072 :
3073 : real(kind=dbl_kind), dimension(:), intent(in) :: &
3074 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
3075 : phi ! ice layer liquid fraction
3076 :
3077 : real(kind=dbl_kind), intent(in) :: &
3078 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
3079 : hpond , & ! melt pond thickness (m) ! LCOV_EXCL_LINE
3080 : apond , & ! melt pond area fraction of category (-) ! LCOV_EXCL_LINE
3081 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
3082 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
3083 : dt ! time step (s)
3084 :
3085 : real(kind=dbl_kind), intent(out) :: &
3086 : w ! vertical flushing Darcy flow rate (m s-1)
3087 :
3088 : real(kind=dbl_kind), parameter :: &
3089 : advection_limit = 0.005_dbl_kind ! limit to fraction of brine in
3090 : ! any layer that can be advected
3091 :
3092 : real(kind=dbl_kind) :: &
3093 : perm , & ! ice layer permeability (m2) ! LCOV_EXCL_LINE
3094 : ice_mass , & ! mass of ice (kg m-2) ! LCOV_EXCL_LINE
3095 : perm_harm , & ! harmonic mean of ice permeability (m2) ! LCOV_EXCL_LINE
3096 : hocn , & ! ocean surface height above ice base (m) ! LCOV_EXCL_LINE
3097 : hbrine , & ! brine surface height above ice base (m) ! LCOV_EXCL_LINE
3098 : w_down_max , & ! maximum downward flushing Darcy flow rate (m s-1) ! LCOV_EXCL_LINE
3099 : phi_min , & ! minimum porosity in the mush ! LCOV_EXCL_LINE
3100 : wlimit , & ! limit to w to avoid advecting all brine in layer ! LCOV_EXCL_LINE
3101 : dhhead ! hydraulic head (m)
3102 :
3103 : integer(kind=int_kind) :: &
3104 : k ! ice layer index
3105 :
3106 : character(len=*),parameter :: subname='(flushing_velocity)'
3107 :
3108 : ! initialize
3109 1408028216 : w = c0
3110 :
3111 : ! only flush if ponds are active
3112 1408028216 : if (tr_pond) then
3113 :
3114 1387334676 : ice_mass = c0
3115 1387334676 : perm_harm = c0
3116 1387334676 : phi_min = c1
3117 :
3118 11098677408 : do k = 1, nilyr
3119 :
3120 : ! liquid fraction
3121 : !phi = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
3122 9711342732 : phi_min = min(phi_min,phi(k))
3123 :
3124 : ! permeability
3125 9711342732 : perm = permeability(phi(k))
3126 :
3127 : ! ice mass
3128 : ice_mass = ice_mass + phi(k) * icepack_mushy_density_brine(liquidus_brine_salinity_mush(zTin(k))) + &
3129 9711342732 : (c1 - phi(k)) * rhoi
3130 :
3131 : ! permeability harmonic mean
3132 11098677408 : perm_harm = perm_harm + c1 / (perm + 1e-30_dbl_kind)
3133 :
3134 : enddo ! k
3135 :
3136 1387334676 : ice_mass = ice_mass * hilyr
3137 :
3138 1387334676 : perm_harm = real(nilyr,dbl_kind) / perm_harm
3139 :
3140 : ! calculate ocean surface height above bottom of ice
3141 1387334676 : hocn = (ice_mass + hpond * apond * rhofresh + hsn * rhos) / rhow
3142 :
3143 : ! calculate brine height above bottom of ice
3144 1387334676 : hbrine = hin + hpond
3145 :
3146 : ! pressure head
3147 1387334676 : dhhead = max(hbrine - hocn,c0)
3148 :
3149 : ! darcy flow through ice
3150 1387334676 : w = (perm_harm * rhow * gravit * (dhhead / hin)) / viscosity_dyn
3151 :
3152 : ! maximum down flow to drain pond
3153 1387334676 : w_down_max = (hpond * apond) / dt
3154 :
3155 : ! limit flow
3156 1387334676 : w = min(w,w_down_max)
3157 :
3158 : ! limit amount of brine that can be advected out of any particular layer
3159 1387334676 : wlimit = (advection_limit * phi_min * hilyr) / dt
3160 :
3161 1387334676 : if (abs(w) > puny) then
3162 347258867 : w = w * max(min(abs(wlimit/w),c1),c0)
3163 : else
3164 1040075809 : w = c0
3165 : endif
3166 :
3167 1387334676 : w = max(w, c0)
3168 :
3169 : endif
3170 :
3171 1408028216 : end subroutine flushing_velocity
3172 :
3173 : !=======================================================================
3174 :
3175 1408028216 : subroutine flush_pond(w, hpond, apond, dt)
3176 :
3177 : ! given a flushing velocity drain the meltponds
3178 :
3179 : real(kind=dbl_kind), intent(in) :: &
3180 : w , & ! vertical flushing Darcy flow rate (m s-1) ! LCOV_EXCL_LINE
3181 : apond , & ! melt pond area fraction of category (-) ! LCOV_EXCL_LINE
3182 : dt ! time step (s)
3183 :
3184 : real(kind=dbl_kind), intent(inout) :: &
3185 : hpond ! melt pond thickness (m)
3186 :
3187 : real(kind=dbl_kind), parameter :: &
3188 : hpond0 = 0.01_dbl_kind
3189 :
3190 : real(kind=dbl_kind) :: &
3191 : lambda_pond ! 1 / macroscopic drainage time scale (s)
3192 :
3193 : character(len=*),parameter :: subname='(flush_pond)'
3194 :
3195 1408028216 : if (tr_pond) then
3196 1387334676 : if (apond > c0 .and. hpond > c0) then
3197 :
3198 : ! flush pond through mush
3199 898572779 : hpond = hpond - w * dt / apond
3200 :
3201 898572779 : hpond = max(hpond, c0)
3202 :
3203 : ! exponential decay of pond
3204 898572779 : lambda_pond = c1 / (tscale_pnd_drain * 24.0_dbl_kind * 3600.0_dbl_kind)
3205 898572779 : hpond = hpond - lambda_pond * dt * (hpond + hpond0)
3206 :
3207 898572779 : hpond = max(hpond, c0)
3208 :
3209 : endif
3210 : endif
3211 :
3212 1408028216 : end subroutine flush_pond
3213 :
3214 : !=======================================================================
3215 :
3216 1408028216 : subroutine flood_ice(hsn, hin, &
3217 : hslyr, hilyr, & ! LCOV_EXCL_LINE
3218 0 : zqsn, zqin, & ! LCOV_EXCL_LINE
3219 1408028216 : phi, dt, & ! LCOV_EXCL_LINE
3220 2816056432 : zSin, Sbr, & ! LCOV_EXCL_LINE
3221 : sss, qocn, & ! LCOV_EXCL_LINE
3222 1408028216 : smice, smliq, & ! LCOV_EXCL_LINE
3223 : snoice, fadvheat)
3224 :
3225 : ! given upwards flushing brine flow calculate amount of snow ice and
3226 : ! convert snow to ice with appropriate properties
3227 :
3228 : real(kind=dbl_kind), intent(in) :: &
3229 : dt , & ! time step (s) ! LCOV_EXCL_LINE
3230 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
3231 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
3232 : sss , & ! sea surface salinity (ppt) ! LCOV_EXCL_LINE
3233 : qocn ! ocean brine enthalpy (J m-3)
3234 :
3235 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
3236 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
3237 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
3238 : zSin , & ! ice layer bulk salinity (ppt) ! LCOV_EXCL_LINE
3239 : phi , & ! ice liquid fraction ! LCOV_EXCL_LINE
3240 : smice , & ! ice mass tracer in snow (kg/m^3) ! LCOV_EXCL_LINE
3241 : smliq ! liquid water mass tracer in snow (kg/m^3)
3242 :
3243 : real(kind=dbl_kind), dimension(:), intent(in) :: &
3244 : Sbr ! ice layer brine salinity (ppt)
3245 :
3246 : real(kind=dbl_kind), intent(inout) :: &
3247 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
3248 : hilyr ! ice layer thickness (m)
3249 :
3250 : real(kind=dbl_kind), intent(out) :: &
3251 : snoice ! snow ice formation
3252 :
3253 : real(kind=dbl_kind), intent(inout) :: &
3254 : fadvheat ! advection heat flux to ocean
3255 :
3256 : real(kind=dbl_kind) :: &
3257 : hin2 , & ! new ice thickness (m) ! LCOV_EXCL_LINE
3258 : hsn2 , & ! new snow thickness (m) ! LCOV_EXCL_LINE
3259 : hilyr2 , & ! new ice layer thickness (m) ! LCOV_EXCL_LINE
3260 : hslyr2 , & ! new snow layer thickness (m) ! LCOV_EXCL_LINE
3261 : dh , & ! thickness of snowice formation (m) ! LCOV_EXCL_LINE
3262 : phi_snowice , & ! liquid fraction of new snow ice ! LCOV_EXCL_LINE
3263 : rho_snowice , & ! density of snowice (kg m-3) ! LCOV_EXCL_LINE
3264 : zSin_snowice , & ! bulk salinity of new snowice (ppt) ! LCOV_EXCL_LINE
3265 : zqin_snowice , & ! ice enthalpy of new snowice (J m-3) ! LCOV_EXCL_LINE
3266 : zqsn_snowice , & ! snow enthalpy of snow that is becoming snowice (J m-3) ! LCOV_EXCL_LINE
3267 : freeboard_density , & ! negative of ice surface freeboard times the ocean density (kg m-2) ! LCOV_EXCL_LINE
3268 : ice_mass , & ! mass of the ice (kg m-2) ! LCOV_EXCL_LINE
3269 : ! snow_mass , & ! mass of the ice (kg m-2) ! LCOV_EXCL_LINE
3270 : rho_ocn , & ! density of the ocean (kg m-3) ! LCOV_EXCL_LINE
3271 : ice_density , & ! density of ice layer (kg m-3) ! LCOV_EXCL_LINE
3272 : hadded , & ! thickness rate of water used from ocean (m/s) ! LCOV_EXCL_LINE
3273 : wadded , & ! mass rate of water used from ocean (kg/m^2/s) ! LCOV_EXCL_LINE
3274 : eadded ! energy rate of water used from ocean (W/m^2)
3275 :
3276 : ! real(kind=dbl_kind) :: &
3277 : ! sadded ! salt rate of water used from ocean (kg/m^2/s)
3278 :
3279 : integer :: &
3280 : k ! vertical index
3281 :
3282 : character(len=*),parameter :: subname='(flood_ice)'
3283 :
3284 1408028216 : snoice = c0
3285 :
3286 : ! check we have snow
3287 1408028216 : if (hsn > puny) then
3288 :
3289 1211582191 : rho_ocn = icepack_mushy_density_brine(sss)
3290 :
3291 : ! ice mass
3292 1211582191 : ice_mass = c0
3293 9692657528 : do k = 1, nilyr
3294 8481075337 : ice_density = min(phi(k) * icepack_mushy_density_brine(Sbr(k)) + (c1 - phi(k)) * rhoi,rho_ocn)
3295 9692657528 : ice_mass = ice_mass + ice_density
3296 : enddo ! k
3297 1211582191 : ice_mass = ice_mass * hilyr
3298 :
3299 : ! for now, do not use variable snow density
3300 : ! snow_mass = c0
3301 : ! if (snwgrain) then
3302 : ! do k = 1,nslyr
3303 : ! snow_mass = snow_mass + (smice(k) + smliq(k)) * hslyr
3304 : ! enddo
3305 :
3306 : ! ! negative freeboard times ocean density
3307 : ! freeboard_density = max(ice_mass + snow_mass - hin * rho_ocn, c0) ! may not be BFB
3308 :
3309 : ! if (freeboard_density > c0) then ! ice is flooded
3310 : ! phi_snowice = (c1 - snow_mass / hsn / rhoi) ! sea ice fraction of newly formed snow-ice
3311 : ! ! density of newly formed snowice
3312 : ! ! use rhos instead of (c1-phi_snowice)*rhoi to conserve ice and liquid
3313 : ! ! snow tracers when rhos = smice + smliq
3314 : ! rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
3315 : ! endif
3316 : ! else
3317 : ! snow_mass = rhos * hsn
3318 : ! negative freeboard times ocean density
3319 1211582191 : freeboard_density = max(ice_mass + hsn * rhos - hin * rho_ocn, c0)
3320 :
3321 1211582191 : if (freeboard_density > c0) then ! ice is flooded
3322 59440891 : phi_snowice = (c1 - rhos / rhoi) ! sea ice fraction of newly formed snow-ice
3323 : ! density of newly formed snow-ice
3324 59440891 : rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
3325 : endif ! freeboard_density > c0
3326 : ! endif ! snwgrain
3327 :
3328 1211582191 : if (freeboard_density > c0) then ! ice is flooded
3329 :
3330 : ! calculate thickness of new ice added
3331 59440891 : dh = freeboard_density / (rho_ocn - rho_snowice + rhos)
3332 59440891 : dh = max(min(dh,hsn),c0)
3333 :
3334 : ! enthalpy of snow that becomes snow-ice
3335 59440891 : call enthalpy_snow_snowice(dh, hsn, zqsn, zqsn_snowice)
3336 59440891 : if (icepack_warnings_aborted(subname)) return
3337 :
3338 : ! change thicknesses
3339 59440891 : hin2 = hin + dh
3340 59440891 : hsn2 = hsn - dh
3341 :
3342 59440891 : hilyr2 = hin2 / real(nilyr,dbl_kind)
3343 59440891 : hslyr2 = hsn2 / real(nslyr,dbl_kind)
3344 :
3345 : ! properties of new snow ice
3346 59440891 : zSin_snowice = phi_snowice * sss
3347 59440891 : zqin_snowice = phi_snowice * qocn + zqsn_snowice
3348 :
3349 : ! change snow properties
3350 59440891 : call update_vertical_tracers_snow(zqsn, hslyr, hslyr2)
3351 59440891 : if (icepack_warnings_aborted(subname)) return
3352 :
3353 59440891 : if (snwgrain .and. hslyr2 > puny) then
3354 768203 : call update_vertical_tracers_snow(smice, hslyr, hslyr2)
3355 768203 : call update_vertical_tracers_snow(smliq, hslyr, hslyr2)
3356 768203 : if (icepack_warnings_aborted(subname)) return
3357 : endif
3358 :
3359 : ! change ice properties
3360 : call update_vertical_tracers_ice(zqin, hilyr, hilyr2, &
3361 59440891 : hin, hin2, zqin_snowice)
3362 59440891 : if (icepack_warnings_aborted(subname)) return
3363 : call update_vertical_tracers_ice(zSin, hilyr, hilyr2, &
3364 59440891 : hin, hin2, zSin_snowice)
3365 59440891 : if (icepack_warnings_aborted(subname)) return
3366 : call update_vertical_tracers_ice(phi, hilyr, hilyr2, &
3367 59440891 : hin, hin2, phi_snowice)
3368 59440891 : if (icepack_warnings_aborted(subname)) return
3369 :
3370 : ! change thicknesses
3371 59440891 : hilyr = hilyr2
3372 59440891 : hslyr = hslyr2
3373 59440891 : snoice = dh
3374 :
3375 59440891 : hadded = (dh * phi_snowice) / dt
3376 59440891 : wadded = hadded * rhoi
3377 59440891 : eadded = hadded * qocn
3378 : ! sadded = wadded * ice_ref_salinity * p001
3379 :
3380 : ! conservation
3381 59440891 : fadvheat = fadvheat - eadded
3382 :
3383 : endif
3384 :
3385 : endif
3386 :
3387 : end subroutine flood_ice
3388 :
3389 : !=======================================================================
3390 :
3391 59440891 : subroutine enthalpy_snow_snowice(dh, hsn, zqsn, zqsn_snowice)
3392 :
3393 : ! determine enthalpy of the snow being converted to snow ice
3394 :
3395 : real(kind=dbl_kind), intent(in) :: &
3396 : dh , & ! thickness of new snowice formation (m) ! LCOV_EXCL_LINE
3397 : hsn ! initial snow thickness
3398 :
3399 : real(kind=dbl_kind), dimension(:), intent(in) :: &
3400 : zqsn ! snow layer enthalpy (J m-3)
3401 :
3402 : real(kind=dbl_kind), intent(out) :: &
3403 : zqsn_snowice ! enthalpy of snow becoming snowice (J m-3)
3404 :
3405 : real(kind=dbl_kind) :: &
3406 : rnlyr ! real value of number of snow layers turning to snowice
3407 :
3408 : integer(kind=int_kind) :: &
3409 : nlyr , & ! no of snow layers completely converted to snowice ! LCOV_EXCL_LINE
3410 : k ! snow layer index
3411 :
3412 : character(len=*),parameter :: subname='(enthalpy_snow_snowice)'
3413 :
3414 59440891 : zqsn_snowice = c0
3415 :
3416 : ! snow depth and snow layers affected by snowice formation
3417 59440891 : if (hsn > puny) then
3418 59440891 : rnlyr = (dh / hsn) * nslyr
3419 59440891 : nlyr = min(floor(rnlyr),nslyr-1) ! nlyr=0 if nslyr=1
3420 :
3421 : ! loop over full snow layers affected
3422 : ! not executed if nlyr=0
3423 59736443 : do k = nslyr, nslyr-nlyr+1, -1
3424 1494662 : zqsn_snowice = zqsn_snowice + zqsn(k) / rnlyr
3425 : enddo ! k
3426 :
3427 : ! partially converted snow layer
3428 : zqsn_snowice = zqsn_snowice + &
3429 59440891 : ((rnlyr - real(nlyr,dbl_kind)) / rnlyr) * zqsn(nslyr-nlyr)
3430 : endif
3431 :
3432 59440891 : end subroutine enthalpy_snow_snowice
3433 :
3434 : !=======================================================================
3435 :
3436 60977297 : subroutine update_vertical_tracers_snow(trc, hlyr1, hlyr2)
3437 :
3438 : ! given some snow ice formation regrid snow layers
3439 :
3440 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
3441 : trc ! vertical tracer
3442 :
3443 : real(kind=dbl_kind), intent(in) :: &
3444 : hlyr1 , & ! old cell thickness ! LCOV_EXCL_LINE
3445 : hlyr2 ! new cell thickness
3446 :
3447 : real(kind=dbl_kind), dimension(1:nslyr) :: &
3448 121954594 : trc2 ! temporary array for updated tracer
3449 :
3450 : ! vertical indexes for old and new grid
3451 : integer(kind=int_kind) :: &
3452 : k1 , & ! vertical index for old grid ! LCOV_EXCL_LINE
3453 : k2 ! vertical index for new grid
3454 :
3455 : real(kind=dbl_kind) :: &
3456 : z1a , & ! lower boundary of old cell ! LCOV_EXCL_LINE
3457 : z1b , & ! upper boundary of old cell ! LCOV_EXCL_LINE
3458 : z2a , & ! lower boundary of new cell ! LCOV_EXCL_LINE
3459 : z2b , & ! upper boundary of new cell ! LCOV_EXCL_LINE
3460 : overlap ! overlap between old and new cell
3461 :
3462 : character(len=*),parameter :: subname='(update_vertical_tracers_snow)'
3463 :
3464 : ! loop over new grid cells
3465 132035412 : do k2 = 1, nslyr
3466 :
3467 : ! initialize new tracer
3468 71058115 : trc2(k2) = c0
3469 :
3470 : ! calculate upper and lower boundary of new cell
3471 71058115 : z2a = (k2 - 1) * hlyr2
3472 71058115 : z2b = k2 * hlyr2
3473 :
3474 : ! loop over old grid cells
3475 190797828 : do k1 = 1, nslyr
3476 :
3477 : ! calculate upper and lower boundary of old cell
3478 119739713 : z1a = (k1 - 1) * hlyr1
3479 119739713 : z1b = k1 * hlyr1
3480 :
3481 : ! calculate overlap between old and new cell
3482 119739713 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3483 :
3484 : ! aggregate old grid cell contribution to new cell
3485 190797828 : trc2(k2) = trc2(k2) + overlap * trc(k1)
3486 :
3487 : enddo ! k1
3488 :
3489 : ! renormalize new grid cell
3490 132035412 : trc2(k2) = trc2(k2) / hlyr2
3491 :
3492 : enddo ! k2
3493 :
3494 : ! update vertical tracer array with the adjusted tracer
3495 132035412 : trc = trc2
3496 :
3497 60977297 : end subroutine update_vertical_tracers_snow
3498 :
3499 : !=======================================================================
3500 :
3501 178322673 : subroutine update_vertical_tracers_ice(trc, hlyr1, hlyr2, &
3502 : h1, h2, trc0)
3503 :
3504 : ! given some snow ice formation regrid ice layers
3505 :
3506 : real(kind=dbl_kind), dimension(:), intent(inout) :: &
3507 : trc ! vertical tracer
3508 :
3509 : real(kind=dbl_kind), intent(in) :: &
3510 : hlyr1 , & ! old cell thickness ! LCOV_EXCL_LINE
3511 : hlyr2 , & ! new cell thickness ! LCOV_EXCL_LINE
3512 : h1 , & ! old total thickness ! LCOV_EXCL_LINE
3513 : h2 , & ! new total thickness ! LCOV_EXCL_LINE
3514 : trc0 ! tracer value of added snow ice on ice top
3515 :
3516 : real(kind=dbl_kind), dimension(1:nilyr) :: &
3517 356645346 : trc2 ! temporary array for updated tracer
3518 :
3519 : integer(kind=int_kind) :: &
3520 : k1 , & ! vertical indexes for old grid ! LCOV_EXCL_LINE
3521 : k2 ! vertical indexes for new grid
3522 :
3523 : real(kind=dbl_kind) :: &
3524 : z1a , & ! lower boundary of old cell ! LCOV_EXCL_LINE
3525 : z1b , & ! upper boundary of old cell ! LCOV_EXCL_LINE
3526 : z2a , & ! lower boundary of new cell ! LCOV_EXCL_LINE
3527 : z2b , & ! upper boundary of new cell ! LCOV_EXCL_LINE
3528 : overlap ! overlap between old and new cell
3529 :
3530 : character(len=*),parameter :: subname='(update_vertical_tracers_ice)'
3531 :
3532 : ! loop over new grid cells
3533 1426581384 : do k2 = 1, nilyr
3534 :
3535 : ! initialize new tracer
3536 1248258711 : trc2(k2) = c0
3537 :
3538 : ! calculate upper and lower boundary of new cell
3539 1248258711 : z2a = (k2 - 1) * hlyr2
3540 1248258711 : z2b = k2 * hlyr2
3541 :
3542 : ! calculate upper and lower boundary of added snow ice at top
3543 1248258711 : z1a = c0
3544 1248258711 : z1b = h2 - h1
3545 :
3546 : ! calculate overlap between added ice and new cell
3547 1248258711 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3548 :
3549 : ! aggregate added ice contribution to new cell
3550 1248258711 : trc2(k2) = trc2(k2) + overlap * trc0
3551 :
3552 : ! loop over old grid cells
3553 9986069688 : do k1 = 1, nilyr
3554 :
3555 : ! calculate upper and lower boundary of old cell
3556 8737810977 : z1a = (k1 - 1) * hlyr1 + h2 - h1
3557 8737810977 : z1b = k1 * hlyr1 + h2 - h1
3558 :
3559 : ! calculate overlap between old and new cell
3560 8737810977 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3561 :
3562 : ! aggregate old grid cell contribution to new cell
3563 9986069688 : trc2(k2) = trc2(k2) + overlap * trc(k1)
3564 :
3565 : enddo ! k1
3566 :
3567 : ! renormalize new grid cell
3568 1426581384 : trc2(k2) = trc2(k2) / hlyr2
3569 :
3570 : enddo ! k2
3571 :
3572 : ! update vertical tracer array with the adjusted tracer
3573 1426581384 : trc = trc2
3574 :
3575 178322673 : end subroutine update_vertical_tracers_ice
3576 :
3577 : !=======================================================================
3578 :
3579 : end module icepack_therm_mushy
3580 :
3581 : !=======================================================================
|
