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 609320240 : subroutine temperature_changes_salinity(dt, &
43 : nilyr, nslyr, &
44 : rhoa, flw, &
45 : potT, Qa, &
46 : shcoef, lhcoef, &
47 : fswsfc, fswint, &
48 1218640480 : Sswabs, Iswabs, &
49 : hilyr, hslyr, &
50 : apond, hpond, &
51 1218640480 : zqin, zTin, &
52 1218640480 : zqsn, zTsn, &
53 609320240 : 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 1413349498 : zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep
119 1413349498 : zSin0 , & ! internal ice layer salinities (ppt) at start of timestep
120 1413349498 : phi , & ! liquid fraction
121 1413349498 : km , & ! ice conductivity (W m-1 K-1)
122 1478252504 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
123 :
124 : real(kind=dbl_kind), dimension(1:nilyr+1) :: &
125 1445801001 : Sbr , & ! brine salinity (ppt)
126 1445801001 : qbr ! brine enthalpy (J m-3)
127 :
128 : real(kind=dbl_kind), dimension(0:nilyr) :: &
129 1445801001 : q ! upward interface vertical Darcy flow (m s-1)
130 :
131 : real(kind=dbl_kind), dimension(1:nslyr) :: &
132 1218640480 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
133 706674749 : ks ! snow conductivity (W m-1 K-1)
134 :
135 : real(kind=dbl_kind) :: &
136 32451503 : Tsf0 , & ! ice/snow surface temperature (C) at start of timestep
137 32451503 : hin , & ! ice thickness (m)
138 32451503 : hsn , & ! snow thickness (m)
139 32451503 : hslyr_min , & ! minimum snow layer thickness (m)
140 32451503 : w , & ! vertical flushing Darcy velocity (m/s)
141 32451503 : qocn , & ! ocean brine enthalpy (J m-3)
142 32451503 : qpond , & ! melt pond brine enthalpy (J m-3)
143 32451503 : Spond ! melt pond salinity (ppt)
144 :
145 64903006 : 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 609320240 : fadvheat = c0
156 609320240 : snoice = c0
157 :
158 609320240 : Tsf0 = Tsf
159 1218640480 : zqsn0 = zqsn
160 4874561920 : zqin0 = zqin
161 4874561920 : zSin0 = zSin
162 :
163 609320240 : Spond = c0
164 609320240 : qpond = enthalpy_brine(c0)
165 :
166 609320240 : hslyr_min = hs_min / real(nslyr, dbl_kind)
167 :
168 609320240 : lsnow = (hslyr > hslyr_min)
169 :
170 609320240 : hin = hilyr * real(nilyr,dbl_kind)
171 :
172 609320240 : qocn = enthalpy_brine(Tbot)
173 :
174 609320240 : if (lsnow) then
175 531603793 : hsn = hslyr * real(nslyr,dbl_kind)
176 : else
177 77716447 : hsn = c0
178 : endif
179 :
180 4874561920 : do k = 1, nilyr
181 4874561920 : 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 609320240 : dt, w)
191 609320240 : 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 609320240 : hilyr, hin)
201 609320240 : if (icepack_warnings_aborted(subname)) return
202 :
203 : ! calculate the conductivities
204 609320240 : call conductivity_mush_array(nilyr, zqin0, zSin0, km)
205 609320240 : 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 609320240 : if (sw_redist) then
221 :
222 0 : dt_rhoi_hlyr = dt / (rhoi*hilyr)
223 :
224 0 : do k = 1, nilyr
225 :
226 0 : Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc
227 :
228 0 : Tmlt = liquidus_temperature_mush(zSin(k))
229 :
230 0 : if (zTin(k) <= Tmlt - sw_dtemp) then
231 0 : ci = cp_ice - Lfresh * Tmlt / (zTin(k)**2)
232 0 : Iswabs_tmp = min(Iswabs(k), &
233 0 : sw_frac*(Tmlt-zTin(k))*ci/dt_rhoi_hlyr)
234 : endif
235 0 : if (Iswabs_tmp < puny) Iswabs_tmp = c0
236 :
237 0 : dswabs = min(Iswabs(k) - Iswabs_tmp, fswint)
238 :
239 0 : fswsfc = fswsfc + dswabs
240 0 : fswint = fswint - dswabs
241 0 : Iswabs(k) = Iswabs_tmp
242 :
243 : enddo
244 :
245 : endif
246 :
247 609320240 : if (lsnow) then
248 : ! case with snow
249 :
250 : ! calculate the snow conductivities
251 531603793 : call conductivity_snow_array(ks)
252 531603793 : 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 531603793 : flatn, fsurfn )
278 :
279 531603793 : 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 1063207586 : do k = 1, nslyr
287 1063207586 : zTsn(k) = temperature_snow(zqsn(k))
288 : enddo ! k
289 :
290 4252830344 : do k = 1, nilyr
291 3721226551 : zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
292 3721226551 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
293 4252830344 : 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 77716447 : flatn, fsurfn )
323 :
324 77716447 : 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 621731576 : do k = 1, nilyr
332 544015129 : zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
333 544015129 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
334 621731576 : qbr(k) = enthalpy_brine(zTin(k))
335 : enddo ! k
336 :
337 : endif
338 :
339 : ! drain ponds from flushing
340 609320240 : call flush_pond(w, hpond, apond, dt)
341 609320240 : 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 609320240 : snoice, fadvheat)
352 609320240 : if (icepack_warnings_aborted(subname)) return
353 :
354 : end subroutine temperature_changes_salinity
355 :
356 : !=======================================================================
357 :
358 531603793 : subroutine two_stage_solver_snow(nilyr, nslyr, &
359 : Tsf, Tsf0, &
360 1063207586 : zqsn, zqsn0, &
361 1063207586 : zqin, zqin0, &
362 1063207586 : zSin, zSin0, &
363 531603793 : zTsn, &
364 531603793 : zTin, &
365 531603793 : phi, Tbot, &
366 531603793 : km, ks, &
367 531603793 : q, dSdt, &
368 : w, dt, &
369 : fswint, fswsfc, &
370 : rhoa, flw, &
371 : potT, Qa, &
372 : shcoef, lhcoef, &
373 531603793 : 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 30172660 : fcondtop1 , & ! first stage downward cond flux at top surface (W m-2)
454 30172660 : fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop
455 30172660 : 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 531603793 : 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 511870466 : w )
488 1043474259 : if (icepack_warnings_aborted(subname)) return
489 :
490 : ! halt if solver failed
491 511870466 : if (icepack_warnings_aborted(subname)) return
492 :
493 : ! check if solution is consistent - surface should still be cold
494 511870466 : if (Tsf < dTemp_errmax) then
495 :
496 : ! solution is consistent - have solution so finish
497 477649810 : return
498 :
499 : else
500 :
501 : ! solution is inconsistent - surface is warmer than melting
502 : ! resolve assuming surface is melting
503 34220656 : Tsf1 = Tsf
504 :
505 : ! reset the solution to initial values
506 34220656 : Tsf = c0
507 68441312 : zqsn = zqsn0
508 273765248 : zqin = zqin0
509 273765248 : 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 34220656 : w )
534 34220656 : if (icepack_warnings_aborted(subname)) return
535 :
536 : ! halt if solver failed
537 34220656 : 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 34220656 : if (fcondtop - fsurfn < ferrmax) then
543 :
544 : ! solution is consistent - have solution so finish
545 34220656 : 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 19733327 : 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 19733327 : w )
588 :
589 19733327 : if (icepack_warnings_aborted(subname)) return
590 :
591 : ! halt if solver failed
592 19733327 : 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 19733327 : if (fcondtop - fsurfn < ferrmax) then
598 :
599 : ! solution is consistent - have solution so finish
600 19082927 : return
601 :
602 : else
603 :
604 : ! solution is inconsistent - resolve assuming other surface condition
605 : ! assume surface is cold
606 650400 : fcondtop1 = fcondtop
607 650400 : fsurfn1 = fsurfn
608 :
609 : ! reset the solution to initial values
610 650400 : Tsf = Tsf0
611 1300800 : zqsn = zqsn0
612 5203200 : zqin = zqin0
613 5203200 : 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 650400 : w )
638 650400 : if (icepack_warnings_aborted(subname)) return
639 :
640 : ! halt if solver failed
641 650400 : if (icepack_warnings_aborted(subname)) return
642 :
643 : ! check if solution is consistent - surface should be cold
644 650400 : if (Tsf < dTemp_errmax) then
645 :
646 : ! solution is consistent - have solution so finish
647 650400 : 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 77716447 : subroutine two_stage_solver_nosnow(nilyr, nslyr, &
670 : Tsf, Tsf0, &
671 77716447 : zqsn, &
672 155432894 : zqin, zqin0, &
673 155432894 : zSin, zSin0, &
674 77716447 : zTsn, &
675 77716447 : zTin, &
676 77716447 : phi, Tbot, &
677 77716447 : km, ks, &
678 77716447 : q, dSdt, &
679 : w, dt, &
680 : fswint, fswsfc, &
681 : rhoa, flw, &
682 : potT, Qa, &
683 : shcoef, lhcoef, &
684 77716447 : 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 2278843 : Tmlt , & ! upper ice layer melting temperature (C)
764 2278843 : fcondtop1 , & ! first stage downward cond flux at top surface (W m-2)
765 2278843 : fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop
766 2278843 : Tsf1 ! first stage ice surface temperature (C)
767 :
768 : character(len=*),parameter :: subname='(two_stage_solver_nosnow)'
769 :
770 : ! initial surface melting temperature
771 77716447 : Tmlt = liquidus_temperature_mush(zSin0(1))
772 :
773 : ! determine if surface is initially cold or melting
774 77716447 : 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 62332236 : w )
801 140048683 : if (icepack_warnings_aborted(subname)) return
802 :
803 : ! halt if solver failed
804 62332236 : if (icepack_warnings_aborted(subname)) return
805 :
806 : ! check if solution is consistent - surface should still be cold
807 62332236 : if (Tsf < Tmlt + dTemp_errmax) then
808 :
809 : ! solution is consistent - have solution so finish
810 41136293 : return
811 :
812 : else
813 : ! solution is inconsistent - surface is warmer than melting
814 : ! resolve assuming surface is melting
815 21195943 : Tsf1 = Tsf
816 :
817 : ! reset the solution to initial values
818 21195943 : Tsf = liquidus_temperature_mush(zSin0(1))
819 169567544 : zqin = zqin0
820 169567544 : 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 21195943 : w )
845 21195943 : if (icepack_warnings_aborted(subname)) return
846 :
847 : ! halt if solver failed
848 21195943 : 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 21195943 : if (fcondtop - fsurfn < ferrmax) then
854 :
855 : ! solution is consistent - have solution so finish
856 21195943 : 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 15384211 : 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 15384211 : w )
899 15384211 : if (icepack_warnings_aborted(subname)) return
900 :
901 : ! halt if solver failed
902 15384211 : 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 15384211 : if (fcondtop - fsurfn < ferrmax) then
908 :
909 : ! solution is consistent - have solution so finish
910 15221204 : return
911 :
912 : else
913 :
914 : ! solution is inconsistent - resolve assuming other surface condition
915 : ! assume surface is cold
916 163007 : fcondtop1 = fcondtop
917 163007 : fsurfn1 = fsurfn
918 :
919 : ! reset the solution to initial values
920 163007 : Tsf = Tsf0
921 1304056 : zqin = zqin0
922 1304056 : 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 163007 : w )
947 163007 : if (icepack_warnings_aborted(subname)) return
948 :
949 : ! halt if solver failed
950 163007 : if (icepack_warnings_aborted(subname)) return
951 :
952 : ! check if solution is consistent - surface should be cold
953 163007 : if (Tsf < Tmlt + dTemp_errmax) then
954 :
955 : ! solution is consistent - have solution so finish
956 163007 : 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 665550246 : subroutine prep_picard(nilyr, nslyr, &
1048 665550246 : lsnow, zqsn, &
1049 665550246 : zqin, zSin, &
1050 : hilyr, hslyr, &
1051 665550246 : km, ks, &
1052 1331100492 : zTin, zTsn, &
1053 665550246 : Sbr, phi, &
1054 1331100492 : 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 5324401968 : do k = 1, nilyr
1092 4658851722 : zTin(k) = icepack_mushy_temperature_mush(zqin(k), zSin(k))
1093 4658851722 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
1094 5324401968 : phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
1095 : enddo ! k
1096 :
1097 665550246 : if (lsnow) then
1098 :
1099 1132949698 : do k = 1, nslyr
1100 1132949698 : zTsn(k) = temperature_snow(zqsn(k))
1101 : enddo ! k
1102 :
1103 : endif ! lsnow
1104 :
1105 : ! interface distances
1106 665550246 : call calc_intercell_thickness(nilyr, nslyr, lsnow, hilyr, hslyr, dxp)
1107 665550246 : if (icepack_warnings_aborted(subname)) return
1108 :
1109 : ! interface conductivities
1110 : call calc_intercell_conductivity(lsnow, nilyr, nslyr, &
1111 665550246 : km, ks, hilyr, hslyr, kcstar)
1112 665550246 : 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 665550246 : einit)
1120 665550246 : if (icepack_warnings_aborted(subname)) return
1121 :
1122 : end subroutine prep_picard
1123 :
1124 : !=======================================================================
1125 :
1126 665550246 : subroutine picard_solver(nilyr, nslyr, &
1127 : lsnow, lcold, &
1128 665550246 : Tsf, zqsn, &
1129 665550246 : zqin, zSin, &
1130 665550246 : zTin, zTsn, &
1131 665550246 : phi, dt, &
1132 : hilyr, hslyr, &
1133 665550246 : km, ks, &
1134 665550246 : 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 1331100492 : 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 1551941388 : Sbr , & ! ice layer brine salinity (ppt)
1212 1551941388 : qbr , & ! ice layer brine enthalpy (J m-3)
1213 1551941388 : zTin0 , & ! ice layer temperature (C) at start of timestep
1214 1551941388 : zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep
1215 1551941388 : zSin0 , & ! ice layer bulk salinity (ppt) at start of timestep
1216 1551941388 : zTin_prev ! ice layer temperature at previous iteration
1217 :
1218 : real(kind=dbl_kind), dimension(nslyr) :: &
1219 1331100492 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
1220 1331100492 : zTsn0 , & ! snow layer temperature (C) at start of timestep
1221 1331100492 : zTsn_prev ! snow layer temperature at previous iteration
1222 :
1223 : real(kind=dbl_kind), dimension(nslyr+nilyr+1) :: &
1224 1699168652 : dxp , & ! distances between grid points (m)
1225 1070425222 : kcstar ! interface conductivities (W m-1 K-1)
1226 :
1227 : real(kind=dbl_kind) :: &
1228 36806816 : Tsf0 , & ! snow surface temperature (C) at start of timestep
1229 36806816 : dfsurfn_dTsf , & ! derivative of net flux to top surface, excluding fcondtopn
1230 36806816 : dflwoutn_dTsf , & ! derivative of longwave flux wrt surface temperature
1231 36806816 : dfsensn_dTsf , & ! derivative of sensible heat flux wrt surface temperature
1232 36806816 : dflatn_dTsf , & ! derivative of latent heat flux wrt surface temperature
1233 36806816 : Tsf_prev , & ! snow surface temperature at previous iteration
1234 36806816 : einit , & ! initial total energy (J)
1235 36806816 : 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 665550246 : 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 665550246 : einit)
1260 665550246 : if (icepack_warnings_aborted(subname)) return
1261 :
1262 665550246 : Tsf0 = Tsf
1263 5324401968 : zqin0 = zqin
1264 1331100492 : zqsn0 = zqsn
1265 5324401968 : zTin0 = zTin
1266 1331100492 : zTsn0 = zTsn
1267 5324401968 : zSin0 = zSin
1268 :
1269 : ! set prev variables
1270 665550246 : Tsf_prev = Tsf
1271 1331100492 : zTsn_prev = zTsn
1272 5324401968 : zTin_prev = zTin
1273 :
1274 : ! picard iteration
1275 1490258537 : 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 1490258537 : flatn, fsurfn)
1284 1490258537 : 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 1490258537 : dfsensn_dTsf, dflatn_dTsf)
1291 1490258537 : 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 1490258537 : zTin, zTsn)
1306 1490258537 : if (icepack_warnings_aborted(subname)) return
1307 :
1308 : ! update brine enthalpy
1309 1490258537 : call picard_updates_enthalpy(nilyr, zTin, qbr)
1310 1490258537 : if (icepack_warnings_aborted(subname)) return
1311 :
1312 : ! drainage fluxes
1313 : call picard_drainage_fluxes(fadvheat_nit, q, &
1314 0 : qbr, qocn, &
1315 1490258537 : nilyr)
1316 1490258537 : if (icepack_warnings_aborted(subname)) return
1317 :
1318 : ! flushing fluxes
1319 : call picard_flushing_fluxes(nilyr, &
1320 : fadvheat_nit, w, &
1321 0 : qbr, &
1322 1490258537 : qpond)
1323 1490258537 : 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 1490258537 : fadvheat_nit)
1340 1490258537 : if (icepack_warnings_aborted(subname)) return
1341 :
1342 1490258537 : if (lconverged) exit
1343 :
1344 824708291 : Tsf_prev = Tsf
1345 1649416582 : zTsn_prev = zTsn
1346 7263216574 : zTin_prev = zTin
1347 :
1348 : enddo picard
1349 :
1350 665550246 : fadvheat = fadvheat_nit
1351 :
1352 : ! update the picard iterants
1353 0 : call picard_updates(nilyr, zTin, &
1354 665550246 : Sbr, qbr)
1355 665550246 : 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 665550246 : dt, nilyr)
1363 665550246 : 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 665550246 : flatn, fsurfn)
1372 665550246 : if (icepack_warnings_aborted(subname)) return
1373 :
1374 : ! if not converged
1375 665550246 : 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 1490258537 : subroutine check_picard_convergence(nilyr, nslyr, &
1456 : lsnow, &
1457 : lconverged, &
1458 : Tsf, Tsf_prev, &
1459 2980517074 : zTin, zTin_prev,&
1460 2980517074 : zTsn, zTsn_prev,&
1461 1490258537 : phi, Tbot, &
1462 1490258537 : zqin, zqsn, &
1463 1490258537 : 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 82004691 : ferr , & ! energy flux error
1514 82004691 : efinal , & ! initial total energy (J) at iteration
1515 82004691 : dzTsn , & ! change in snow temperature (C) between iterations
1516 82004691 : dzTin , & ! change in ice temperature (C) between iterations
1517 82004691 : 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 1490258537 : phi)
1526 1490258537 : if (icepack_warnings_aborted(subname)) return
1527 :
1528 : call total_energy_content(lsnow, &
1529 : nilyr, nslyr, &
1530 0 : zqin, zqsn, &
1531 : hilyr, hslyr, &
1532 1490258537 : efinal)
1533 1490258537 : 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 1490258537 : zTin, zTin_prev, dzTin)
1539 1490258537 : if (icepack_warnings_aborted(subname)) return
1540 :
1541 1490258537 : fcondbot = c2 * km(nilyr) * ((zTin(nilyr) - Tbot) / hilyr)
1542 :
1543 1490258537 : if (lsnow) then
1544 1267651203 : fcondtop = c2 * ks(1) * ((Tsf - zTsn(1)) / hslyr)
1545 : else
1546 222607334 : fcondtop = c2 * km(1) * ((Tsf - zTin(1)) / hilyr)
1547 : endif
1548 :
1549 1490258537 : 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 1490258537 : abs(ferr) < 0.9_dbl_kind*ferrmax)
1555 :
1556 : end subroutine check_picard_convergence
1557 :
1558 : !=======================================================================
1559 :
1560 1490258537 : subroutine picard_drainage_fluxes(fadvheat, q, &
1561 1490258537 : 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 1490258537 : fadvheat = c0
1585 :
1586 : ! calculate fluxes from base upwards
1587 10431809759 : do k = 1, nilyr-1
1588 :
1589 10431809759 : fadvheat = fadvheat - q(k) * (qbr(k+1) - qbr(k))
1590 :
1591 : enddo ! k
1592 :
1593 1490258537 : k = nilyr
1594 :
1595 1490258537 : fadvheat = fadvheat - q(k) * (qocn - qbr(k))
1596 :
1597 1490258537 : end subroutine picard_drainage_fluxes
1598 :
1599 : !=======================================================================
1600 :
1601 1490258537 : subroutine picard_flushing_fluxes(nilyr, &
1602 : fadvheat, w, &
1603 1490258537 : 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 1490258537 : fadvheat = fadvheat + w * (qbr(nilyr) - qpond)
1622 :
1623 1490258537 : end subroutine picard_flushing_fluxes
1624 :
1625 : !=======================================================================
1626 :
1627 1490258537 : subroutine maximum_variables_changes(lsnow, &
1628 : Tsf, Tsf_prev, dTsf, &
1629 1490258537 : zTsn, zTsn_prev, dzTsn, &
1630 1490258537 : 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 1490258537 : dTsf = abs(Tsf - Tsf_prev)
1653 :
1654 1490258537 : if (lsnow) then
1655 2535302406 : dzTsn = maxval(abs(zTsn - zTsn_prev))
1656 : else ! lsnow
1657 222607334 : dzTsn = c0
1658 : endif ! lsnow
1659 :
1660 11922068296 : dzTin = maxval(abs(zTin - zTin_prev))
1661 :
1662 1490258537 : end subroutine maximum_variables_changes
1663 :
1664 : !=======================================================================
1665 :
1666 2155808783 : subroutine total_energy_content(lsnow, &
1667 : nilyr, nslyr, &
1668 2155808783 : 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 2155808783 : energy = c0
1696 :
1697 2155808783 : if (lsnow) then
1698 :
1699 3668252104 : do k = 1, nslyr
1700 :
1701 3668252104 : energy = energy + hslyr * zqsn(k)
1702 :
1703 : enddo ! k
1704 :
1705 : endif ! lsnow
1706 :
1707 17246470264 : do k = 1, nilyr
1708 :
1709 17246470264 : energy = energy + hilyr * zqin(k)
1710 :
1711 : enddo ! k
1712 :
1713 2155808783 : end subroutine total_energy_content
1714 :
1715 : !=======================================================================
1716 :
1717 665550246 : subroutine picard_updates(nilyr, zTin, &
1718 1331100492 : 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 5324401968 : do k = 1, nilyr
1738 :
1739 4658851722 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
1740 5324401968 : qbr(k) = enthalpy_brine(zTin(k))
1741 :
1742 : enddo ! k
1743 :
1744 665550246 : end subroutine picard_updates
1745 :
1746 : !=======================================================================
1747 :
1748 1490258537 : 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 11922068296 : do k = 1, nilyr
1767 :
1768 11922068296 : qbr(k) = enthalpy_brine(zTin(k))
1769 :
1770 : enddo ! k
1771 :
1772 1490258537 : end subroutine picard_updates_enthalpy
1773 :
1774 : !=======================================================================
1775 :
1776 1490258537 : subroutine picard_final(lsnow, &
1777 : nilyr, nslyr, &
1778 2980517074 : zqin, zqsn, &
1779 2980517074 : zTin, zTsn, &
1780 1490258537 : 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 11922068296 : do k = 1, nilyr
1804 11922068296 : zqin(k) = enthalpy_mush_liquid_fraction(zTin(k), phi(k))
1805 : enddo ! k
1806 :
1807 1490258537 : if (lsnow) then
1808 :
1809 2535302406 : do k = 1, nslyr
1810 2535302406 : zqsn(k) = enthalpy_snow(zTsn(k))
1811 : enddo ! k
1812 :
1813 : endif ! lsnow
1814 :
1815 1490258537 : end subroutine picard_final
1816 :
1817 : !=======================================================================
1818 :
1819 665550246 : 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 665550246 : if (lsnow) then
1841 :
1842 566474849 : dxp(1) = hslyr / c2
1843 :
1844 566474849 : do l = 2, nslyr
1845 :
1846 566474849 : dxp(l) = hslyr
1847 :
1848 : enddo ! l
1849 :
1850 566474849 : dxp(nslyr+1) = (hilyr + hslyr) / c2
1851 :
1852 3965323943 : do l = nslyr+2, nilyr+nslyr
1853 :
1854 3965323943 : dxp(l) = hilyr
1855 :
1856 : enddo ! l
1857 :
1858 566474849 : dxp(nilyr+nslyr+1) = hilyr / c2
1859 :
1860 : else ! lsnow
1861 :
1862 99075397 : dxp(1) = hilyr / c2
1863 :
1864 693527779 : do l = 2, nilyr
1865 :
1866 693527779 : dxp(l) = hilyr
1867 :
1868 : enddo ! l
1869 :
1870 99075397 : dxp(nilyr+1) = hilyr / c2
1871 :
1872 198150794 : do l = nilyr+2, nilyr+nslyr+1
1873 :
1874 198150794 : dxp(l) = c0
1875 :
1876 : enddo ! l
1877 :
1878 : endif ! lsnow
1879 :
1880 665550246 : end subroutine calc_intercell_thickness
1881 :
1882 : !=======================================================================
1883 :
1884 665550246 : subroutine calc_intercell_conductivity(lsnow, &
1885 : nilyr, nslyr, &
1886 665550246 : km, ks, &
1887 : hilyr, hslyr, &
1888 665550246 : 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 36806816 : 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 665550246 : if (lsnow) then
1918 :
1919 566474849 : kcstar(1) = ks(1)
1920 :
1921 566474849 : do l = 2, nslyr
1922 :
1923 0 : k = l
1924 566474849 : kcstar(l) = (c2 * ks(k) * ks(k-1)) / (ks(k) + ks(k-1))
1925 :
1926 : enddo ! l
1927 :
1928 566474849 : fe = hilyr / (hilyr + hslyr)
1929 566474849 : kcstar(nslyr+1) = c1 / ((c1 - fe) / ks(nslyr) + fe / km(1))
1930 :
1931 3965323943 : do k = 2, nilyr
1932 :
1933 3398849094 : l = k + nslyr
1934 3965323943 : kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
1935 :
1936 : enddo ! k
1937 :
1938 566474849 : kcstar(nilyr+nslyr+1) = km(nilyr)
1939 :
1940 : else ! lsnow
1941 :
1942 99075397 : kcstar(1) = km(1)
1943 :
1944 693527779 : do k = 2, nilyr
1945 :
1946 594452382 : l = k
1947 693527779 : kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
1948 :
1949 : enddo ! k
1950 :
1951 99075397 : kcstar(nilyr+1) = km(nilyr)
1952 :
1953 198150794 : do l = nilyr+2, nilyr+nslyr+1
1954 :
1955 198150794 : kcstar(l) = c0
1956 :
1957 : enddo ! l
1958 :
1959 : endif ! lsnow
1960 :
1961 665550246 : end subroutine calc_intercell_conductivity
1962 :
1963 : !=======================================================================
1964 :
1965 1490258537 : subroutine solve_heat_conduction(lsnow, lcold, &
1966 : nilyr, nslyr, &
1967 : Tsf, Tbot, &
1968 1490258537 : zqin0, zqsn0, &
1969 1490258537 : phi, dt, &
1970 : qpond, qocn, &
1971 1490258537 : q, w, &
1972 : hilyr, hslyr, &
1973 1490258537 : dxp, kcstar, &
1974 1490258537 : Iswabs, Sswabs, &
1975 : fsurfn, dfsurfn_dTsf, &
1976 1490258537 : 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 3636554602 : Ap , & ! diagonal of tridiagonal matrix
2022 3636554602 : As , & ! lower off-diagonal of tridiagonal matrix
2023 3800563984 : An , & ! upper off-diagonal of tridiagonal matrix
2024 3636554602 : b , & ! right hand side of matrix solve
2025 2310305447 : 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 1490258537 : if (lsnow) then
2034 :
2035 1267651203 : 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 1159743396 : dt)
2049 1159743396 : 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 107907807 : dt)
2064 107907807 : if (icepack_warnings_aborted(subname)) return
2065 :
2066 : endif ! lcold
2067 :
2068 : else ! lsnow
2069 :
2070 222607334 : 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 149447785 : dt)
2084 149447785 : 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 73159549 : dt)
2099 73159549 : 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 1490258537 : An(1:nyn), Ap(1:nyn), As(1:nyn), b(1:nyn), T(1:nyn), nyn)
2108 1490258537 : if (icepack_warnings_aborted(subname)) return
2109 :
2110 : call update_temperatures(lsnow, lcold, &
2111 : nilyr, nslyr, &
2112 0 : T, Tsf, &
2113 1490258537 : zTin, zTsn)
2114 1490258537 : if (icepack_warnings_aborted(subname)) return
2115 :
2116 : end subroutine solve_heat_conduction
2117 :
2118 : !=======================================================================
2119 :
2120 1490258537 : subroutine update_temperatures(lsnow, lcold, &
2121 : nilyr, nslyr, &
2122 1490258537 : T, Tsf, &
2123 1490258537 : 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 1490258537 : if (lsnow) then
2150 :
2151 1267651203 : if (lcold) then
2152 :
2153 1159743396 : Tsf = T(1)
2154 :
2155 2319486792 : do k = 1, nslyr
2156 1159743396 : l = k + 1
2157 2319486792 : zTsn(k) = T(l)
2158 : enddo ! k
2159 :
2160 9277947168 : do k = 1, nilyr
2161 8118203772 : l = k + nslyr + 1
2162 9277947168 : zTin(k) = T(l)
2163 : enddo ! k
2164 :
2165 : else ! lcold
2166 :
2167 215815614 : do k = 1, nslyr
2168 107907807 : l = k
2169 215815614 : zTsn(k) = T(l)
2170 : enddo ! k
2171 :
2172 863262456 : do k = 1, nilyr
2173 755354649 : l = k + nslyr
2174 863262456 : zTin(k) = T(l)
2175 : enddo ! k
2176 :
2177 : endif ! lcold
2178 :
2179 : else ! lsnow
2180 :
2181 222607334 : if (lcold) then
2182 :
2183 149447785 : Tsf = T(1)
2184 :
2185 1195582280 : do k = 1, nilyr
2186 1046134495 : l = k + 1
2187 1195582280 : zTin(k) = T(l)
2188 : enddo ! k
2189 :
2190 : else ! lcold
2191 :
2192 585276392 : do k = 1, nilyr
2193 512116843 : l = k
2194 585276392 : zTin(k) = T(l)
2195 : enddo ! k
2196 :
2197 : endif ! lcold
2198 :
2199 : endif ! lsnow
2200 :
2201 1490258537 : end subroutine update_temperatures
2202 :
2203 : !=======================================================================
2204 :
2205 146319098 : subroutine matrix_elements_nosnow_melt(Ap, As, An, b, nyn, &
2206 : nilyr, &
2207 : Tsf, Tbot, &
2208 73159549 : zqin0, &
2209 : qpond, qocn, &
2210 73159549 : phi, q, &
2211 : w, &
2212 : hilyr, &
2213 73159549 : dxp, kcstar, &
2214 73159549 : 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 73159549 : k = 1
2258 73159549 : l = k
2259 :
2260 4403370 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2261 8806740 : kcstar(k+1) / dxp(k+1) + &
2262 4403370 : kcstar(k) / dxp(k) + &
2263 4403370 : q(k) * cp_ocn * rhow + &
2264 77562919 : w * cp_ocn * rhow
2265 13210110 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2266 81966289 : q(k) * cp_ocn * rhow
2267 73159549 : An(l) = c0
2268 4403370 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2269 4403370 : (kcstar(k) / dxp(k)) * Tsf + &
2270 73159549 : w * qpond
2271 :
2272 : ! interior ice layers
2273 438957294 : do k = 2, nilyr-1
2274 :
2275 365797745 : l = k
2276 :
2277 22016850 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2278 44033700 : kcstar(k+1) / dxp(k+1) + &
2279 22016850 : kcstar(k) / dxp(k) + &
2280 22016850 : q(k) * cp_ocn * rhow + &
2281 387814595 : w * cp_ocn * rhow
2282 66050550 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2283 409831445 : q(k) * cp_ocn * rhow
2284 22016850 : An(l) = -kcstar(k) / dxp(k) - &
2285 365797745 : w * cp_ocn * rhow
2286 438957294 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2287 :
2288 : enddo ! k
2289 :
2290 : ! bottom layer
2291 73159549 : k = nilyr
2292 73159549 : l = k
2293 :
2294 4403370 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2295 8806740 : kcstar(k+1) / dxp(k+1) + &
2296 4403370 : kcstar(k) / dxp(k) + &
2297 4403370 : q(k) * cp_ocn * rhow + &
2298 77562919 : w * cp_ocn * rhow
2299 73159549 : As(l) = c0
2300 4403370 : An(l) = -kcstar(k) / dxp(k) - &
2301 73159549 : w * cp_ocn * rhow
2302 4403370 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2303 8806740 : (kcstar(k+1) * Tbot) / dxp(k+1) + &
2304 81966289 : q(k) * qocn
2305 :
2306 73159549 : nyn = nilyr
2307 :
2308 73159549 : end subroutine matrix_elements_nosnow_melt
2309 :
2310 : !=======================================================================
2311 :
2312 298895570 : subroutine matrix_elements_nosnow_cold(Ap, As, An, b, nyn, &
2313 : nilyr, &
2314 : Tsf, Tbot, &
2315 149447785 : zqin0, &
2316 : qpond, qocn, &
2317 149447785 : phi, q, &
2318 : w, &
2319 : hilyr, &
2320 149447785 : dxp, kcstar, &
2321 149447785 : 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 149447785 : l = 1
2368 149447785 : Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
2369 149447785 : As(l) = kcstar(1) / dxp(1)
2370 149447785 : An(l) = c0
2371 149447785 : b (l) = dfsurfn_dTsf * Tsf - fsurfn
2372 :
2373 : ! surface layer
2374 149447785 : k = 1
2375 149447785 : l = k + 1
2376 :
2377 3722439 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2378 7444878 : kcstar(k+1) / dxp(k+1) + &
2379 3722439 : kcstar(k) / dxp(k) + &
2380 3722439 : q(k) * cp_ocn * rhow + &
2381 153170224 : w * cp_ocn * rhow
2382 11167317 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2383 156892663 : q(k) * cp_ocn * rhow
2384 149447785 : An(l) = -kcstar(k) / dxp(k)
2385 3722439 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2386 149447785 : w * qpond
2387 :
2388 : ! interior ice layers
2389 896686710 : do k = 2, nilyr-1
2390 :
2391 747238925 : l = k + 1
2392 :
2393 18612195 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2394 37224390 : kcstar(k+1) / dxp(k+1) + &
2395 18612195 : kcstar(k) / dxp(k) + &
2396 18612195 : q(k) * cp_ocn * rhow + &
2397 765851120 : w * cp_ocn * rhow
2398 55836585 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2399 784463315 : q(k) * cp_ocn * rhow
2400 18612195 : An(l) = -kcstar(k) / dxp(k) - &
2401 747238925 : w * cp_ocn * rhow
2402 896686710 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2403 :
2404 : enddo ! k
2405 :
2406 : ! bottom layer
2407 149447785 : k = nilyr
2408 149447785 : l = k + 1
2409 :
2410 3722439 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2411 7444878 : kcstar(k+1) / dxp(k+1) + &
2412 3722439 : kcstar(k) / dxp(k) + &
2413 3722439 : q(k) * cp_ocn * rhow + &
2414 153170224 : w * cp_ocn * rhow
2415 149447785 : As(l) = c0
2416 3722439 : An(l) = -kcstar(k) / dxp(k) - &
2417 149447785 : w * cp_ocn * rhow
2418 3722439 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2419 7444878 : (kcstar(k+1) * Tbot) / dxp(k+1) + &
2420 156892663 : q(k) * qocn
2421 :
2422 149447785 : nyn = nilyr + 1
2423 :
2424 149447785 : end subroutine matrix_elements_nosnow_cold
2425 :
2426 : !=======================================================================
2427 :
2428 215815614 : subroutine matrix_elements_snow_melt(Ap, As, An, b, nyn, &
2429 : nilyr, nslyr, &
2430 : Tsf, Tbot, &
2431 215815614 : zqin0, zqsn0, &
2432 : qpond, qocn, &
2433 107907807 : phi, q, &
2434 : w, &
2435 : hilyr, hslyr, &
2436 107907807 : dxp, kcstar, &
2437 215815614 : 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 107907807 : k = 1
2485 107907807 : l = k
2486 :
2487 9898576 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2488 19797152 : kcstar(l+1) / dxp(l+1) + &
2489 127704959 : kcstar(l) / dxp(l)
2490 107907807 : As(l) = -kcstar(l+1) / dxp(l+1)
2491 107907807 : An(l) = c0
2492 9898576 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k) + &
2493 107907807 : (kcstar(l) * Tsf) / dxp(l)
2494 :
2495 : ! interior snow layers
2496 107907807 : 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 107907807 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2506 :
2507 : enddo ! k
2508 :
2509 : ! top ice layer
2510 107907807 : k = 1
2511 107907807 : l = nslyr + k
2512 :
2513 9898576 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2514 19797152 : kcstar(l+1) / dxp(l+1) + &
2515 9898576 : kcstar(l) / dxp(l) + &
2516 9898576 : q(k) * cp_ocn * rhow + &
2517 117806383 : w * cp_ocn * rhow
2518 29695728 : As(l) = -kcstar(l+1) / dxp(l+1) - &
2519 127704959 : q(k) * cp_ocn * rhow
2520 107907807 : An(l) = -kcstar(l) / dxp(l)
2521 9898576 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2522 107907807 : w * qpond
2523 :
2524 : ! interior ice layers
2525 647446842 : do k = 2, nilyr-1
2526 :
2527 539539035 : l = nslyr + k
2528 :
2529 49492880 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2530 98985760 : kcstar(l+1) / dxp(l+1) + &
2531 49492880 : kcstar(l) / dxp(l) + &
2532 49492880 : q(k) * cp_ocn * rhow + &
2533 589031915 : w * cp_ocn * rhow
2534 148478640 : As(l) = -kcstar(l+1) / dxp(l+1) - &
2535 638524795 : q(k) * cp_ocn * rhow
2536 49492880 : An(l) = -kcstar(l) / dxp(l) - &
2537 539539035 : w * cp_ocn * rhow
2538 647446842 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2539 :
2540 : enddo ! k
2541 :
2542 : ! bottom layer
2543 107907807 : k = nilyr
2544 107907807 : l = nilyr + nslyr
2545 :
2546 9898576 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2547 19797152 : kcstar(l+1) / dxp(l+1) + &
2548 9898576 : kcstar(l) / dxp(l) + &
2549 9898576 : q(k) * cp_ocn * rhow + &
2550 117806383 : w * cp_ocn * rhow
2551 107907807 : As(l) = c0
2552 9898576 : An(l) = -kcstar(l) / dxp(l) - &
2553 107907807 : w * cp_ocn * rhow
2554 9898576 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2555 19797152 : (kcstar(l+1) * Tbot) / dxp(l+1) + &
2556 127704959 : q(k) * qocn
2557 :
2558 107907807 : nyn = nilyr + nslyr
2559 :
2560 107907807 : end subroutine matrix_elements_snow_melt
2561 :
2562 : !=======================================================================
2563 :
2564 2319486792 : subroutine matrix_elements_snow_cold(Ap, As, An, b, nyn, &
2565 : nilyr, nslyr, &
2566 : Tsf, Tbot, &
2567 2319486792 : zqin0, zqsn0, &
2568 : qpond, qocn, &
2569 1159743396 : phi, q, &
2570 : w, &
2571 : hilyr, hslyr, &
2572 1159743396 : dxp, kcstar, &
2573 2319486792 : 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 1159743396 : l = 1
2625 1159743396 : Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
2626 1159743396 : As(l) = kcstar(1) / dxp(1)
2627 1159743396 : An(l) = c0
2628 1159743396 : b (l) = dfsurfn_dTsf * Tsf - fsurfn
2629 :
2630 : ! surface layer
2631 1159743396 : k = 1
2632 1159743396 : l = k + 1
2633 1159743396 : m = 1
2634 :
2635 63980306 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2636 127960612 : kcstar(m+1) / dxp(m+1) + &
2637 1287704008 : kcstar(m) / dxp(m)
2638 1159743396 : As(l) = -kcstar(m+1) / dxp(m+1)
2639 1159743396 : An(l) = -kcstar(m) / dxp(m)
2640 1159743396 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2641 :
2642 : ! interior snow layers
2643 1159743396 : 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 1159743396 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2654 :
2655 : enddo ! k
2656 :
2657 : ! top ice layer
2658 1159743396 : k = 1
2659 1159743396 : l = nslyr + k + 1
2660 1159743396 : m = k + nslyr
2661 :
2662 63980306 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2663 127960612 : kcstar(m+1) / dxp(m+1) + &
2664 63980306 : kcstar(m) / dxp(m) + &
2665 63980306 : q(k) * cp_ocn * rhow + &
2666 1223723702 : w * cp_ocn * rhow
2667 191940918 : As(l) = -kcstar(m+1) / dxp(m+1) - &
2668 1287704008 : q(k) * cp_ocn * rhow
2669 1159743396 : An(l) = -kcstar(m) / dxp(m)
2670 63980306 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2671 1159743396 : w * qpond
2672 :
2673 : ! interior ice layers
2674 6958460376 : do k = 2, nilyr-1
2675 :
2676 5798716980 : l = nslyr + k + 1
2677 5798716980 : m = k + nslyr
2678 :
2679 319901530 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2680 639803060 : kcstar(m+1) / dxp(m+1) + &
2681 319901530 : kcstar(m) / dxp(m) + &
2682 319901530 : q(k) * cp_ocn * rhow + &
2683 6118618510 : w * cp_ocn * rhow
2684 959704590 : As(l) = -kcstar(m+1) / dxp(m+1) - &
2685 6438520040 : q(k) * cp_ocn * rhow
2686 319901530 : An(l) = -kcstar(m) / dxp(m) - &
2687 5798716980 : w * cp_ocn * rhow
2688 6958460376 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2689 :
2690 : enddo ! k
2691 :
2692 : ! bottom layer
2693 1159743396 : k = nilyr
2694 1159743396 : l = nilyr + nslyr + 1
2695 1159743396 : m = k + nslyr
2696 :
2697 63980306 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2698 127960612 : kcstar(m+1) / dxp(m+1) + &
2699 63980306 : kcstar(m) / dxp(m) + &
2700 63980306 : q(k) * cp_ocn * rhow + &
2701 1223723702 : w * cp_ocn * rhow
2702 1159743396 : As(l) = c0
2703 63980306 : An(l) = -kcstar(m) / dxp(m) - &
2704 1159743396 : w * cp_ocn * rhow
2705 63980306 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2706 127960612 : (kcstar(m+1) * Tbot) / dxp(m+1) + &
2707 1287704008 : q(k) * qocn
2708 :
2709 1159743396 : nyn = nilyr + nslyr + 1
2710 :
2711 1159743396 : end subroutine matrix_elements_snow_cold
2712 :
2713 : !=======================================================================
2714 :
2715 1331100492 : subroutine solve_salinity(zSin, Sbr, &
2716 : Spond, sss, &
2717 1331100492 : 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 960004774 : zSin0
2749 :
2750 : character(len=*),parameter :: subname='(solve_salinity)'
2751 :
2752 5324401968 : zSin0 = zSin
2753 :
2754 665550246 : k = 1
2755 36806816 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2756 73613632 : ((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + &
2757 36806816 : dSdt(k) + &
2758 702357062 : (w * (Spond - Sbr(k))) / hilyr) * dt)
2759 :
2760 3993301476 : do k = 2, nilyr-1
2761 :
2762 184034080 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2763 368068160 : ((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + &
2764 184034080 : dSdt(k) + &
2765 4177335556 : (w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
2766 :
2767 : enddo ! k
2768 :
2769 665550246 : k = nilyr
2770 36806816 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2771 36806816 : ((q(k) * (sss - Sbr(k))) / hilyr + &
2772 36806816 : dSdt(k) + &
2773 665550246 : (w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
2774 :
2775 :
2776 5324401968 : 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 1490258537 : 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 3636554602 : cp , & ! modified upper off-diagonal vector
2822 3636554602 : 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 1490258537 : cp(1) = c(1) / b(1)
2831 11518393606 : do i = 2, n-1
2832 11518393606 : cp(i) = c(i) / (b(i) - cp(i-1)*a(i))
2833 : enddo
2834 :
2835 1490258537 : dp(1) = d(1) / b(1)
2836 13008652143 : do i = 2, n
2837 13008652143 : dp(i) = (d(i) - dp(i-1)*a(i)) / (b(i) - cp(i-1)*a(i))
2838 : enddo
2839 :
2840 : ! back substitution
2841 1490258537 : x(n) = dp(n)
2842 13008652143 : do i = n-1,1,-1
2843 13008652143 : x(i) = dp(i) - cp(i)*x(i+1)
2844 : enddo
2845 :
2846 1490258537 : end subroutine tdma_solve_sparse
2847 :
2848 : !=======================================================================
2849 : ! Effect of salinity
2850 : !=======================================================================
2851 :
2852 8494022880 : 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 8494022880 : perm = 3.0e-8_dbl_kind * max(phi - phic, c0)**3
2869 :
2870 8494022880 : end function permeability
2871 :
2872 : !=======================================================================
2873 :
2874 609320240 : subroutine explicit_flow_velocities(nilyr, zSin, &
2875 609320240 : zTin, Tsf, &
2876 609320240 : Tbot, q, &
2877 1218640480 : dSdt, Sbr, &
2878 609320240 : 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 1413349498 : phi ! ice layer liquid fraction
2919 :
2920 : real(kind=dbl_kind), dimension(0:nilyr) :: &
2921 1445801001 : rho ! ice layer brine density (kg m-3)
2922 :
2923 : real(kind=dbl_kind) :: &
2924 32451503 : rho_ocn , & ! ocean density (kg m-3)
2925 32451503 : perm_min , & ! minimum permeability from layer to ocean (m2)
2926 32451503 : perm_harm , & ! harmonic mean of permeability from layer to ocean (m2)
2927 32451503 : rho_sum , & ! sum of the brine densities from layer to ocean (kg m-3)
2928 32451503 : rho_pipe , & ! density of the brine in the channel (kg m-3)
2929 32451503 : z , & ! distance to layer from top surface (m)
2930 32451503 : perm , & ! ice layer permeability (m2)
2931 32451503 : drho , & ! brine density difference between layer and ocean (kg m-3)
2932 32451503 : Ra , & ! local mush Rayleigh number
2933 32451503 : rn , & ! real value of number of layers considered
2934 32451503 : L , & ! thickness of the layers considered (m)
2935 32451503 : dx , & ! horizontal size of convective flow (m)
2936 32451503 : dx2 , & ! square of the horizontal size of convective flow (m2)
2937 32451503 : Am , & ! A parameter for mush
2938 32451503 : Bm , & ! B parameter for mush
2939 32451503 : Ap , & ! A parameter for channel
2940 32451503 : Bp , & ! B parameter for channel
2941 32451503 : qlimit , & ! limit to vertical Darcy flow for numerical stability
2942 32451503 : dS_guess , & ! expected bulk salinity without limits
2943 32451503 : alpha , & ! desalination limiting factor
2944 32451503 : 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 4874561920 : do k = 1, nilyr
2953 :
2954 4265241680 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
2955 4265241680 : phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
2956 4265241680 : qbr(k) = enthalpy_brine(zTin(k))
2957 4874561920 : rho(k) = icepack_mushy_density_brine(Sbr(k))
2958 :
2959 : enddo ! k
2960 :
2961 609320240 : rho(0) = rho(1)
2962 :
2963 : ! ocean conditions
2964 609320240 : Sbr(nilyr+1) = sss
2965 609320240 : qbr(nilyr+1) = qocn
2966 609320240 : rho_ocn = icepack_mushy_density_brine(sss)
2967 :
2968 : ! initialize accumulated quantities
2969 609320240 : perm_min = bignum
2970 609320240 : perm_harm = c0
2971 609320240 : rho_sum = c0
2972 :
2973 : ! limit to q for numerical stability
2974 609320240 : qlimit = (fracmax * hilyr) / dt
2975 :
2976 : ! no flow through ice top surface
2977 609320240 : q(0) = c0
2978 :
2979 609320240 : ra_constants = gravit / (viscosity_dyn * kappal)
2980 : ! first iterate over layers going up
2981 4874561920 : do k = nilyr, 1, -1
2982 :
2983 : ! vertical position from ice top surface
2984 4265241680 : z = ((real(k, dbl_kind) - p5) / real(nilyr, dbl_kind)) * hin
2985 :
2986 : ! permeabilities
2987 4265241680 : perm = permeability(phi(k))
2988 4265241680 : perm_min = min(perm_min,perm)
2989 4265241680 : perm_harm = perm_harm + (c1 / max(perm,1.0e-30_dbl_kind))
2990 :
2991 : ! densities
2992 4265241680 : rho_sum = rho_sum + rho(k)
2993 : !rho_pipe = rho(k)
2994 4265241680 : rho_pipe = p5 * (rho(k) + rho(k-1))
2995 4265241680 : drho = max(rho(k) - rho_ocn, c0)
2996 :
2997 : ! mush Rayleigh number
2998 4265241680 : Ra = drho * (hin-z) * perm_min * ra_constants
2999 :
3000 : ! height of mush layer to layer k
3001 4265241680 : rn = real(nilyr-k+1,dbl_kind)
3002 4265241680 : L = rn * hilyr
3003 :
3004 : ! horizontal size of convection
3005 4265241680 : dx = L * c2 * aspect_rapid_mode
3006 4265241680 : dx2 = dx**2
3007 :
3008 : ! determine vertical Darcy flow
3009 4265241680 : Am = (dx2 * rn) / (viscosity_dyn * perm_harm)
3010 4265241680 : Bm = (-gravit * rho_sum) / rn
3011 :
3012 4265241680 : Ap = (pi * a_rapid_mode**4) / (c8 * viscosity_dyn)
3013 4265241680 : Bp = -rho_pipe * gravit
3014 :
3015 4265241680 : 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 4265241680 : q(k) = min(q(k) * (max(Ra - Rac_rapid_mode, c0) / (Ra+puny)), qlimit)
3019 :
3020 : ! late stage drainage
3021 227160521 : dSdt(k) = dSdt_slow_mode * (max((zSin(k) - phi_c_slow_mode*Sbr(k)), c0) &
3022 4265241680 : * max((Tbot - Tsf), c0)) / (hin + 0.001_dbl_kind)
3023 :
3024 4265241680 : dSdt(k) = max(dSdt(k), (-zSin(k) * 0.5_dbl_kind) / dt)
3025 :
3026 : ! restrict flows to prevent too much salt loss
3027 4265241680 : dS_guess = (((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + dSdt(k)) * dt) * safety_factor
3028 :
3029 4265241680 : if (abs(dS_guess) < puny) then
3030 1863571806 : alpha = c1
3031 : else
3032 2401669874 : alpha = (zSin_min - zSin(k)) / dS_guess
3033 : endif
3034 :
3035 4265241680 : if (alpha < c0 .or. alpha > c1) alpha = c1
3036 :
3037 4265241680 : q(k) = q(k) * alpha
3038 4874561920 : dSdt(k) = dSdt(k) * alpha
3039 :
3040 : enddo ! k
3041 :
3042 609320240 : end subroutine explicit_flow_velocities
3043 :
3044 : !=======================================================================
3045 : ! Flushing
3046 : !=======================================================================
3047 :
3048 609320240 : subroutine flushing_velocity(zTin, &
3049 609320240 : 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 32451503 : perm , & ! ice layer permeability (m2)
3081 32451503 : ice_mass , & ! mass of ice (kg m-2)
3082 32451503 : perm_harm , & ! harmonic mean of ice permeability (m2)
3083 32451503 : hocn , & ! ocean surface height above ice base (m)
3084 32451503 : hbrine , & ! brine surface height above ice base (m)
3085 32451503 : w_down_max , & ! maximum downward flushing Darcy flow rate (m s-1)
3086 32451503 : phi_min , & ! minimum porosity in the mush
3087 32451503 : wlimit , & ! limit to w to avoid advecting all brine in layer
3088 32451503 : 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 609320240 : w = c0
3097 :
3098 : ! only flush if ponds are active
3099 609320240 : if (tr_pond) then
3100 :
3101 604111600 : ice_mass = c0
3102 604111600 : perm_harm = c0
3103 604111600 : phi_min = c1
3104 :
3105 4832892800 : do k = 1, nilyr
3106 :
3107 : ! liquid fraction
3108 : !phi = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
3109 4228781200 : phi_min = min(phi_min,phi(k))
3110 :
3111 : ! permeability
3112 4228781200 : perm = permeability(phi(k))
3113 :
3114 : ! ice mass
3115 222602961 : ice_mass = ice_mass + phi(k) * icepack_mushy_density_brine(liquidus_brine_salinity_mush(zTin(k))) + &
3116 4228781200 : (c1 - phi(k)) * rhoi
3117 :
3118 : ! permeability harmonic mean
3119 4832892800 : perm_harm = perm_harm + c1 / (perm + 1e-30_dbl_kind)
3120 :
3121 : enddo ! k
3122 :
3123 604111600 : ice_mass = ice_mass * hilyr
3124 :
3125 604111600 : perm_harm = real(nilyr,dbl_kind) / perm_harm
3126 :
3127 : ! calculate ocean surface height above bottom of ice
3128 604111600 : hocn = (ice_mass + hpond * apond * rhow + hsn * rhos) / rhow
3129 :
3130 : ! calculate brine height above bottom of ice
3131 604111600 : hbrine = hin + hpond
3132 :
3133 : ! pressure head
3134 604111600 : dhhead = max(hbrine - hocn,c0)
3135 :
3136 : ! darcy flow through ice
3137 604111600 : w = (perm_harm * rhow * gravit * (dhhead / hin)) / viscosity_dyn
3138 :
3139 : ! maximum down flow to drain pond
3140 604111600 : w_down_max = (hpond * apond) / dt
3141 :
3142 : ! limit flow
3143 604111600 : w = min(w,w_down_max)
3144 :
3145 : ! limit amount of brine that can be advected out of any particular layer
3146 604111600 : wlimit = (advection_limit * phi_min * hilyr) / dt
3147 :
3148 604111600 : if (abs(w) > puny) then
3149 127188275 : w = w * max(min(abs(wlimit/w),c1),c0)
3150 : else
3151 476923325 : w = c0
3152 : endif
3153 :
3154 604111600 : w = max(w, c0)
3155 :
3156 : endif
3157 :
3158 609320240 : end subroutine flushing_velocity
3159 :
3160 : !=======================================================================
3161 :
3162 609320240 : 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 609320240 : if (tr_pond) then
3181 604111600 : if (apond > c0 .and. hpond > c0) then
3182 :
3183 : ! flush pond through mush
3184 386497542 : hpond = hpond - w * dt / apond
3185 :
3186 386497542 : hpond = max(hpond, c0)
3187 :
3188 : ! exponential decay of pond
3189 386497542 : hpond = hpond - lambda_pond * dt * (hpond + hpond0)
3190 :
3191 386497542 : hpond = max(hpond, c0)
3192 :
3193 : endif
3194 : endif
3195 :
3196 609320240 : end subroutine flush_pond
3197 :
3198 : !=======================================================================
3199 :
3200 609320240 : subroutine flood_ice(hsn, hin, &
3201 : nslyr, nilyr, &
3202 : hslyr, hilyr, &
3203 609320240 : zqsn, zqin, &
3204 609320240 : phi, dt, &
3205 1218640480 : 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 32451503 : hin2 , & ! new ice thickness (m)
3244 32451503 : hsn2 , & ! new snow thickness (m)
3245 32451503 : hilyr2 , & ! new ice layer thickness (m)
3246 32451503 : hslyr2 , & ! new snow layer thickness (m)
3247 32451503 : dh , & ! thickness of snowice formation (m)
3248 32451503 : phi_snowice , & ! liquid fraction of new snow ice
3249 32451503 : rho_snowice , & ! density of snowice (kg m-3)
3250 32451503 : zSin_snowice , & ! bulk salinity of new snowice (ppt)
3251 32451503 : zqin_snowice , & ! ice enthalpy of new snowice (J m-2)
3252 32451503 : zqsn_snowice , & ! snow enthalpy of snow thats becoming snowice (J m-2)
3253 32451503 : freeboard_density , & ! negative of ice surface freeboard times the ocean density (kg m-2)
3254 32451503 : ice_mass , & ! mass of the ice (kg m-2)
3255 32451503 : rho_ocn , & ! density of the ocean (kg m-3)
3256 32451503 : ice_density , & ! density of ice layer (kg m-3)
3257 32451503 : hadded , & ! thickness rate of water used from ocean (m/s)
3258 32451503 : wadded , & ! mass rate of water used from ocean (kg/m^2/s)
3259 32451503 : 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 609320240 : snoice = c0
3270 :
3271 : ! check we have snow
3272 609320240 : if (hsn > puny) then
3273 :
3274 531603793 : rho_ocn = icepack_mushy_density_brine(sss)
3275 :
3276 : ! ice mass
3277 531603793 : ice_mass = c0
3278 4252830344 : do k = 1, nilyr
3279 3721226551 : ice_density = min(phi(k) * icepack_mushy_density_brine(Sbr(k)) + (c1 - phi(k)) * rhoi,rho_ocn)
3280 4252830344 : ice_mass = ice_mass + ice_density
3281 : enddo ! k
3282 531603793 : ice_mass = ice_mass * hilyr
3283 :
3284 : ! negative freeboard times ocean density
3285 531603793 : freeboard_density = max(ice_mass + hsn * rhos - hin * rho_ocn, c0)
3286 :
3287 : ! check if have flooded ice
3288 531603793 : if (freeboard_density > c0) then
3289 :
3290 : ! sea ice fraction of newly formed snow ice
3291 51964778 : phi_snowice = (c1 - rhos / rhoi)
3292 :
3293 : ! density of newly formed snowice
3294 51964778 : rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
3295 :
3296 : ! calculate thickness of new ice added
3297 51964778 : dh = freeboard_density / (rho_ocn - rho_snowice + rhos)
3298 51964778 : dh = max(min(dh,hsn),c0)
3299 :
3300 : ! enthalpy of snow that becomes snowice
3301 51964778 : call enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
3302 51964778 : if (icepack_warnings_aborted(subname)) return
3303 :
3304 : ! change thicknesses
3305 51964778 : hin2 = hin + dh
3306 51964778 : hsn2 = hsn - dh
3307 :
3308 51964778 : hilyr2 = hin2 / real(nilyr,dbl_kind)
3309 51964778 : hslyr2 = hsn2 / real(nslyr,dbl_kind)
3310 :
3311 : ! properties of new snow ice
3312 51964778 : zSin_snowice = phi_snowice * sss
3313 51964778 : zqin_snowice = phi_snowice * qocn + zqsn_snowice
3314 :
3315 : ! change snow properties
3316 51964778 : call update_vertical_tracers_snow(nslyr, zqsn, hslyr, hslyr2)
3317 51964778 : if (icepack_warnings_aborted(subname)) return
3318 :
3319 : ! change ice properties
3320 0 : call update_vertical_tracers_ice(nilyr, zqin, hilyr, hilyr2, &
3321 51964778 : hin, hin2, zqin_snowice)
3322 51964778 : if (icepack_warnings_aborted(subname)) return
3323 0 : call update_vertical_tracers_ice(nilyr, zSin, hilyr, hilyr2, &
3324 51964778 : hin, hin2, zSin_snowice)
3325 51964778 : if (icepack_warnings_aborted(subname)) return
3326 0 : call update_vertical_tracers_ice(nilyr, phi, hilyr, hilyr2, &
3327 51964778 : hin, hin2, phi_snowice)
3328 51964778 : if (icepack_warnings_aborted(subname)) return
3329 :
3330 : ! change thicknesses
3331 51964778 : hilyr = hilyr2
3332 51964778 : hslyr = hslyr2
3333 51964778 : snoice = dh
3334 :
3335 51964778 : hadded = (dh * phi_snowice) / dt
3336 51964778 : wadded = hadded * rhoi
3337 51964778 : eadded = hadded * qocn
3338 : ! sadded = wadded * ice_ref_salinity * p001
3339 :
3340 : ! conservation
3341 51964778 : fadvheat = fadvheat - eadded
3342 :
3343 : endif
3344 :
3345 : endif
3346 :
3347 : end subroutine flood_ice
3348 :
3349 : !=======================================================================
3350 :
3351 51964778 : 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 4452112 : 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 51964778 : zqsn_snowice = c0
3378 :
3379 : ! snow depth and snow layers affected by snowice formation
3380 51964778 : if (hsn > puny) then
3381 51964778 : rnlyr = (dh / hsn) * nslyr
3382 51964778 : 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 51964778 : do k = nslyr, nslyr-nlyr+1, -1
3387 51964778 : zqsn_snowice = zqsn_snowice + zqsn(k) / rnlyr
3388 : enddo ! k
3389 :
3390 : ! partially converted snow layer
3391 : zqsn_snowice = zqsn_snowice + &
3392 51964778 : ((rnlyr - real(nlyr,dbl_kind)) / rnlyr) * zqsn(nslyr-nlyr)
3393 : endif
3394 :
3395 51964778 : end subroutine enthalpy_snow_snowice
3396 :
3397 : !=======================================================================
3398 :
3399 51964778 : 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 103929556 : 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 4452112 : z1a , & ! lower boundary of old cell
3423 4452112 : z1b , & ! upper boundary of old cell
3424 4452112 : z2a , & ! lower boundary of new cell
3425 4452112 : z2b , & ! upper boundary of new cell
3426 4452112 : 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 103929556 : do k2 = 1, nslyr
3432 :
3433 : ! initialize new tracer
3434 51964778 : trc2(k2) = c0
3435 :
3436 : ! calculate upper and lower boundary of new cell
3437 51964778 : z2a = (k2 - 1) * hlyr2
3438 51964778 : z2b = k2 * hlyr2
3439 :
3440 : ! loop over old grid cells
3441 103929556 : do k1 = 1, nslyr
3442 :
3443 : ! calculate upper and lower boundary of old cell
3444 51964778 : z1a = (k1 - 1) * hlyr1
3445 51964778 : z1b = k1 * hlyr1
3446 :
3447 : ! calculate overlap between old and new cell
3448 51964778 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3449 :
3450 : ! aggregate old grid cell contribution to new cell
3451 103929556 : trc2(k2) = trc2(k2) + overlap * trc(k1)
3452 :
3453 : enddo ! k1
3454 :
3455 : ! renormalize new grid cell
3456 103929556 : trc2(k2) = trc2(k2) / hlyr2
3457 :
3458 : enddo ! k2
3459 :
3460 : ! update vertical tracer array with the adjusted tracer
3461 103929556 : trc = trc2
3462 :
3463 51964778 : end subroutine update_vertical_tracers_snow
3464 :
3465 : !=======================================================================
3466 :
3467 155894334 : 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 391926684 : 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 13356336 : z1a , & ! lower boundary of old cell
3494 13356336 : z1b , & ! upper boundary of old cell
3495 13356336 : z2a , & ! lower boundary of new cell
3496 13356336 : z2b , & ! upper boundary of new cell
3497 13356336 : 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 1247154672 : do k2 = 1, nilyr
3503 :
3504 : ! initialize new tracer
3505 1091260338 : trc2(k2) = c0
3506 :
3507 : ! calculate upper and lower boundary of new cell
3508 1091260338 : z2a = (k2 - 1) * hlyr2
3509 1091260338 : z2b = k2 * hlyr2
3510 :
3511 : ! calculate upper and lower boundary of added snow ice at top
3512 1091260338 : z1a = c0
3513 1091260338 : z1b = h2 - h1
3514 :
3515 : ! calculate overlap between added ice and new cell
3516 1091260338 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3517 :
3518 : ! aggregate added ice contribution to new cell
3519 1091260338 : trc2(k2) = trc2(k2) + overlap * trc0
3520 :
3521 : ! loop over old grid cells
3522 8730082704 : do k1 = 1, nilyr
3523 :
3524 : ! calculate upper and lower boundary of old cell
3525 7638822366 : z1a = (k1 - 1) * hlyr1 + h2 - h1
3526 7638822366 : z1b = k1 * hlyr1 + h2 - h1
3527 :
3528 : ! calculate overlap between old and new cell
3529 7638822366 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3530 :
3531 : ! aggregate old grid cell contribution to new cell
3532 8730082704 : trc2(k2) = trc2(k2) + overlap * trc(k1)
3533 :
3534 : enddo ! k1
3535 :
3536 : ! renormalize new grid cell
3537 1247154672 : trc2(k2) = trc2(k2) / hlyr2
3538 :
3539 : enddo ! k2
3540 :
3541 : ! update vertical tracer array with the adjusted tracer
3542 1247154672 : trc = trc2
3543 :
3544 155894334 : end subroutine update_vertical_tracers_ice
3545 :
3546 : !=======================================================================
3547 :
3548 : end module icepack_therm_mushy
3549 :
3550 : !=======================================================================
|
