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