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