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