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