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