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 609317819 : subroutine temperature_changes_salinity(dt, &
42 : nilyr, nslyr, &
43 : rhoa, flw, &
44 : potT, Qa, &
45 : shcoef, lhcoef, &
46 : fswsfc, fswint, &
47 609317819 : Sswabs, Iswabs, &
48 : hilyr, hslyr, &
49 : apond, hpond, &
50 1218635638 : zqin, zTin, &
51 1218635638 : zqsn, zTsn, &
52 609317819 : 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 1413344656 : zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep
118 1413344656 : zSin0 , & ! internal ice layer salinities (ppt) at start of timestep
119 1413344656 : phi , & ! liquid fraction
120 1413344656 : km , & ! ice conductivity (W m-1 K-1)
121 1478247662 : dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
122 :
123 : real(kind=dbl_kind), dimension(1:nilyr+1) :: &
124 1445796159 : Sbr , & ! brine salinity (ppt)
125 1445796159 : qbr ! brine enthalpy (J m-3)
126 :
127 : real(kind=dbl_kind), dimension(0:nilyr) :: &
128 1445796159 : q ! upward interface vertical Darcy flow (m s-1)
129 :
130 : real(kind=dbl_kind), dimension(1:nslyr) :: &
131 1218635638 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
132 706672328 : ks ! snow conductivity (W m-1 K-1)
133 :
134 : real(kind=dbl_kind) :: &
135 32451503 : Tsf0 , & ! ice/snow surface temperature (C) at start of timestep
136 32451503 : hin , & ! ice thickness (m)
137 32451503 : hsn , & ! snow thickness (m)
138 32451503 : hslyr_min , & ! minimum snow layer thickness (m)
139 32451503 : w , & ! vertical flushing Darcy velocity (m/s)
140 32451503 : qocn , & ! ocean brine enthalpy (J m-3)
141 32451503 : qpond , & ! melt pond brine enthalpy (J m-3)
142 32451503 : 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 609317819 : fadvheat = c0
153 609317819 : snoice = c0
154 :
155 609317819 : Tsf0 = Tsf
156 1218635638 : zqsn0 = zqsn
157 4874542552 : zqin0 = zqin
158 4874542552 : zSin0 = zSin
159 :
160 609317819 : Spond = c0
161 609317819 : qpond = enthalpy_brine(c0)
162 :
163 609317819 : hslyr_min = hs_min / real(nslyr, dbl_kind)
164 :
165 609317819 : lsnow = (hslyr > hslyr_min)
166 :
167 609317819 : hin = hilyr * real(nilyr,dbl_kind)
168 :
169 609317819 : qocn = enthalpy_brine(Tbot)
170 :
171 609317819 : if (lsnow) then
172 531599226 : hsn = hslyr * real(nslyr,dbl_kind)
173 : else
174 77718593 : hsn = c0
175 : endif
176 :
177 4874542552 : do k = 1, nilyr
178 4874542552 : 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 609317819 : dt, w)
188 609317819 : 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 609317819 : hilyr, hin)
198 609317819 : if (icepack_warnings_aborted(subname)) return
199 :
200 : ! calculate the conductivities
201 609317819 : call conductivity_mush_array(nilyr, zqin0, zSin0, km)
202 609317819 : if (icepack_warnings_aborted(subname)) return
203 :
204 609317819 : if (lsnow) then
205 : ! case with snow
206 :
207 : ! calculate the snow conductivities
208 531599226 : call conductivity_snow_array(ks)
209 531599226 : 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 531599226 : flatn, fsurfn )
235 :
236 531599226 : 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 1063198452 : do k = 1, nslyr
244 1063198452 : zTsn(k) = temperature_snow(zqsn(k))
245 : enddo ! k
246 :
247 4252793808 : do k = 1, nilyr
248 3721194582 : zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
249 3721194582 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
250 4252793808 : 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 77718593 : flatn, fsurfn )
280 :
281 77718593 : 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 621748744 : do k = 1, nilyr
289 544030151 : zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
290 544030151 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
291 621748744 : qbr(k) = enthalpy_brine(zTin(k))
292 : enddo ! k
293 :
294 : endif
295 :
296 : ! drain ponds from flushing
297 609317819 : call flush_pond(w, hpond, apond, dt)
298 609317819 : 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 609317819 : snoice, fadvheat)
309 609317819 : if (icepack_warnings_aborted(subname)) return
310 :
311 : end subroutine temperature_changes_salinity
312 :
313 : !=======================================================================
314 :
315 531599226 : subroutine two_stage_solver_snow(nilyr, nslyr, &
316 : Tsf, Tsf0, &
317 1063198452 : zqsn, zqsn0, &
318 1063198452 : zqin, zqin0, &
319 1063198452 : zSin, zSin0, &
320 531599226 : zTsn, &
321 531599226 : zTin, &
322 531599226 : phi, Tbot, &
323 531599226 : km, ks, &
324 531599226 : q, dSdt, &
325 : w, dt, &
326 : fswint, fswsfc, &
327 : rhoa, flw, &
328 : potT, Qa, &
329 : shcoef, lhcoef, &
330 531599226 : 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 30172660 : fcondtop1 , & ! first stage downward cond flux at top surface (W m-2)
411 30172660 : fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop
412 30172660 : 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 531599226 : 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 511866832 : w )
445 1043466058 : if (icepack_warnings_aborted(subname)) return
446 :
447 : ! halt if solver failed
448 511866832 : if (icepack_warnings_aborted(subname)) return
449 :
450 : ! check if solution is consistent - surface should still be cold
451 511866832 : if (Tsf < dTemp_errmax) then
452 :
453 : ! solution is consistent - have solution so finish
454 477647204 : return
455 :
456 : else
457 :
458 : ! solution is inconsistent - surface is warmer than melting
459 : ! resolve assuming surface is melting
460 34219628 : Tsf1 = Tsf
461 :
462 : ! reset the solution to initial values
463 34219628 : Tsf = c0
464 68439256 : zqsn = zqsn0
465 273757024 : zqin = zqin0
466 273757024 : 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 34219628 : w )
491 34219628 : if (icepack_warnings_aborted(subname)) return
492 :
493 : ! halt if solver failed
494 34219628 : 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 34219628 : if (fcondtop - fsurfn < ferrmax) then
500 :
501 : ! solution is consistent - have solution so finish
502 34219628 : 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 19732394 : 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 19732394 : w )
545 :
546 19732394 : if (icepack_warnings_aborted(subname)) return
547 :
548 : ! halt if solver failed
549 19732394 : 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 19732394 : if (fcondtop - fsurfn < ferrmax) then
555 :
556 : ! solution is consistent - have solution so finish
557 19082042 : return
558 :
559 : else
560 :
561 : ! solution is inconsistent - resolve assuming other surface condition
562 : ! assume surface is cold
563 650352 : fcondtop1 = fcondtop
564 650352 : fsurfn1 = fsurfn
565 :
566 : ! reset the solution to initial values
567 650352 : Tsf = Tsf0
568 1300704 : zqsn = zqsn0
569 5202816 : zqin = zqin0
570 5202816 : 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 650352 : w )
595 650352 : if (icepack_warnings_aborted(subname)) return
596 :
597 : ! halt if solver failed
598 650352 : if (icepack_warnings_aborted(subname)) return
599 :
600 : ! check if solution is consistent - surface should be cold
601 650352 : if (Tsf < dTemp_errmax) then
602 :
603 : ! solution is consistent - have solution so finish
604 650352 : 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 77718593 : subroutine two_stage_solver_nosnow(nilyr, nslyr, &
627 : Tsf, Tsf0, &
628 77718593 : zqsn, &
629 155437186 : zqin, zqin0, &
630 155437186 : zSin, zSin0, &
631 77718593 : zTsn, &
632 77718593 : zTin, &
633 77718593 : phi, Tbot, &
634 77718593 : km, ks, &
635 77718593 : q, dSdt, &
636 : w, dt, &
637 : fswint, fswsfc, &
638 : rhoa, flw, &
639 : potT, Qa, &
640 : shcoef, lhcoef, &
641 77718593 : 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 2278843 : Tmlt , & ! upper ice layer melting temperature (C)
721 2278843 : fcondtop1 , & ! first stage downward cond flux at top surface (W m-2)
722 2278843 : fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop
723 2278843 : Tsf1 ! first stage ice surface temperature (C)
724 :
725 : character(len=*),parameter :: subname='(two_stage_solver_nosnow)'
726 :
727 : ! initial surface melting temperature
728 77718593 : Tmlt = liquidus_temperature_mush(zSin0(1))
729 :
730 : ! determine if surface is initially cold or melting
731 77718593 : 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 62331792 : w )
758 140050385 : if (icepack_warnings_aborted(subname)) return
759 :
760 : ! halt if solver failed
761 62331792 : if (icepack_warnings_aborted(subname)) return
762 :
763 : ! check if solution is consistent - surface should still be cold
764 62331792 : if (Tsf < Tmlt + dTemp_errmax) then
765 :
766 : ! solution is consistent - have solution so finish
767 41135217 : return
768 :
769 : else
770 : ! solution is inconsistent - surface is warmer than melting
771 : ! resolve assuming surface is melting
772 21196575 : Tsf1 = Tsf
773 :
774 : ! reset the solution to initial values
775 21196575 : Tsf = liquidus_temperature_mush(zSin0(1))
776 169572600 : zqin = zqin0
777 169572600 : 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 21196575 : w )
802 21196575 : if (icepack_warnings_aborted(subname)) return
803 :
804 : ! halt if solver failed
805 21196575 : 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 21196575 : if (fcondtop - fsurfn < ferrmax) then
811 :
812 : ! solution is consistent - have solution so finish
813 21196575 : 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 15386801 : 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 15386801 : w )
856 15386801 : if (icepack_warnings_aborted(subname)) return
857 :
858 : ! halt if solver failed
859 15386801 : 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 15386801 : if (fcondtop - fsurfn < ferrmax) then
865 :
866 : ! solution is consistent - have solution so finish
867 15223802 : return
868 :
869 : else
870 :
871 : ! solution is inconsistent - resolve assuming other surface condition
872 : ! assume surface is cold
873 162999 : fcondtop1 = fcondtop
874 162999 : fsurfn1 = fsurfn
875 :
876 : ! reset the solution to initial values
877 162999 : Tsf = Tsf0
878 1303992 : zqin = zqin0
879 1303992 : 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 162999 : w )
904 162999 : if (icepack_warnings_aborted(subname)) return
905 :
906 : ! halt if solver failed
907 162999 : if (icepack_warnings_aborted(subname)) return
908 :
909 : ! check if solution is consistent - surface should be cold
910 162999 : if (Tsf < Tmlt + dTemp_errmax) then
911 :
912 : ! solution is consistent - have solution so finish
913 162999 : 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 665547373 : subroutine prep_picard(nilyr, nslyr, &
1005 665547373 : lsnow, zqsn, &
1006 665547373 : zqin, zSin, &
1007 : hilyr, hslyr, &
1008 665547373 : km, ks, &
1009 1331094746 : zTin, zTsn, &
1010 665547373 : Sbr, phi, &
1011 1331094746 : 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 5324378984 : do k = 1, nilyr
1049 4658831611 : zTin(k) = icepack_mushy_temperature_mush(zqin(k), zSin(k))
1050 4658831611 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
1051 5324378984 : phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
1052 : enddo ! k
1053 :
1054 665547373 : if (lsnow) then
1055 :
1056 1132938412 : do k = 1, nslyr
1057 1132938412 : zTsn(k) = temperature_snow(zqsn(k))
1058 : enddo ! k
1059 :
1060 : endif ! lsnow
1061 :
1062 : ! interface distances
1063 665547373 : call calc_intercell_thickness(nilyr, nslyr, lsnow, hilyr, hslyr, dxp)
1064 665547373 : if (icepack_warnings_aborted(subname)) return
1065 :
1066 : ! interface conductivities
1067 : call calc_intercell_conductivity(lsnow, nilyr, nslyr, &
1068 665547373 : km, ks, hilyr, hslyr, kcstar)
1069 665547373 : 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 665547373 : einit)
1077 665547373 : if (icepack_warnings_aborted(subname)) return
1078 :
1079 : end subroutine prep_picard
1080 :
1081 : !=======================================================================
1082 :
1083 665547373 : subroutine picard_solver(nilyr, nslyr, &
1084 : lsnow, lcold, &
1085 665547373 : Tsf, zqsn, &
1086 665547373 : zqin, zSin, &
1087 665547373 : zTin, zTsn, &
1088 665547373 : phi, dt, &
1089 : hilyr, hslyr, &
1090 665547373 : km, ks, &
1091 665547373 : 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 1331094746 : 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 1551935642 : Sbr , & ! ice layer brine salinity (ppt)
1169 1551935642 : qbr , & ! ice layer brine enthalpy (J m-3)
1170 1551935642 : zTin0 , & ! ice layer temperature (C) at start of timestep
1171 1551935642 : zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep
1172 1551935642 : zSin0 , & ! ice layer bulk salinity (ppt) at start of timestep
1173 1551935642 : zTin_prev ! ice layer temperature at previous iteration
1174 :
1175 : real(kind=dbl_kind), dimension(nslyr) :: &
1176 1331094746 : zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
1177 1331094746 : zTsn0 , & ! snow layer temperature (C) at start of timestep
1178 1331094746 : zTsn_prev ! snow layer temperature at previous iteration
1179 :
1180 : real(kind=dbl_kind), dimension(nslyr+nilyr+1) :: &
1181 1699162906 : dxp , & ! distances between grid points (m)
1182 1070422349 : kcstar ! interface conductivities (W m-1 K-1)
1183 :
1184 : real(kind=dbl_kind) :: &
1185 36806816 : Tsf0 , & ! snow surface temperature (C) at start of timestep
1186 36806816 : dfsurfn_dTsf , & ! derivative of net flux to top surface, excluding fcondtopn
1187 36806816 : dflwoutn_dTsf , & ! derivative of longwave flux wrt surface temperature
1188 36806816 : dfsensn_dTsf , & ! derivative of sensible heat flux wrt surface temperature
1189 36806816 : dflatn_dTsf , & ! derivative of latent heat flux wrt surface temperature
1190 36806816 : Tsf_prev , & ! snow surface temperature at previous iteration
1191 36806816 : einit , & ! initial total energy (J)
1192 36806816 : 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 665547373 : 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 665547373 : einit)
1217 665547373 : if (icepack_warnings_aborted(subname)) return
1218 :
1219 665547373 : Tsf0 = Tsf
1220 5324378984 : zqin0 = zqin
1221 1331094746 : zqsn0 = zqsn
1222 5324378984 : zTin0 = zTin
1223 1331094746 : zTsn0 = zTsn
1224 5324378984 : zSin0 = zSin
1225 :
1226 : ! set prev variables
1227 665547373 : Tsf_prev = Tsf
1228 1331094746 : zTsn_prev = zTsn
1229 5324378984 : zTin_prev = zTin
1230 :
1231 : ! picard iteration
1232 1490253356 : 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 1490253356 : flatn, fsurfn)
1241 1490253356 : 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 1490253356 : dfsensn_dTsf, dflatn_dTsf)
1248 1490253356 : 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 1490253356 : zTin, zTsn)
1263 1490253356 : if (icepack_warnings_aborted(subname)) return
1264 :
1265 : ! update brine enthalpy
1266 1490253356 : call picard_updates_enthalpy(nilyr, zTin, qbr)
1267 1490253356 : if (icepack_warnings_aborted(subname)) return
1268 :
1269 : ! drainage fluxes
1270 : call picard_drainage_fluxes(fadvheat_nit, q, &
1271 0 : qbr, qocn, &
1272 1490253356 : nilyr)
1273 1490253356 : if (icepack_warnings_aborted(subname)) return
1274 :
1275 : ! flushing fluxes
1276 : call picard_flushing_fluxes(nilyr, &
1277 : fadvheat_nit, w, &
1278 0 : qbr, &
1279 1490253356 : qpond)
1280 1490253356 : 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 1490253356 : fadvheat_nit)
1297 1490253356 : if (icepack_warnings_aborted(subname)) return
1298 :
1299 1490253356 : if (lconverged) exit
1300 :
1301 824705983 : Tsf_prev = Tsf
1302 1649411966 : zTsn_prev = zTsn
1303 7263195237 : zTin_prev = zTin
1304 :
1305 : enddo picard
1306 :
1307 665547373 : fadvheat = fadvheat_nit
1308 :
1309 : ! update the picard iterants
1310 0 : call picard_updates(nilyr, zTin, &
1311 665547373 : Sbr, qbr)
1312 665547373 : 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 665547373 : dt, nilyr)
1320 665547373 : 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 665547373 : flatn, fsurfn)
1329 665547373 : if (icepack_warnings_aborted(subname)) return
1330 :
1331 : ! if not converged
1332 665547373 : 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 1490253356 : subroutine check_picard_convergence(nilyr, nslyr, &
1413 : lsnow, &
1414 : lconverged, &
1415 : Tsf, Tsf_prev, &
1416 2980506712 : zTin, zTin_prev,&
1417 2980506712 : zTsn, zTsn_prev,&
1418 1490253356 : phi, Tbot, &
1419 1490253356 : zqin, zqsn, &
1420 1490253356 : 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 82004691 : ferr , & ! energy flux error
1471 82004691 : efinal , & ! initial total energy (J) at iteration
1472 82004691 : dzTsn , & ! change in snow temperature (C) between iterations
1473 82004691 : dzTin , & ! change in ice temperature (C) between iterations
1474 82004691 : 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 1490253356 : phi)
1483 1490253356 : if (icepack_warnings_aborted(subname)) return
1484 :
1485 : call total_energy_content(lsnow, &
1486 : nilyr, nslyr, &
1487 0 : zqin, zqsn, &
1488 : hilyr, hslyr, &
1489 1490253356 : efinal)
1490 1490253356 : 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 1490253356 : zTin, zTin_prev, dzTin)
1496 1490253356 : if (icepack_warnings_aborted(subname)) return
1497 :
1498 1490253356 : fcondbot = c2 * km(nilyr) * ((zTin(nilyr) - Tbot) / hilyr)
1499 :
1500 1490253356 : if (lsnow) then
1501 1267640144 : fcondtop = c2 * ks(1) * ((Tsf - zTsn(1)) / hslyr)
1502 : else
1503 222613212 : fcondtop = c2 * km(1) * ((Tsf - zTin(1)) / hilyr)
1504 : endif
1505 :
1506 1490253356 : 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 1490253356 : abs(ferr) < 0.9_dbl_kind*ferrmax)
1512 :
1513 : end subroutine check_picard_convergence
1514 :
1515 : !=======================================================================
1516 :
1517 1490253356 : subroutine picard_drainage_fluxes(fadvheat, q, &
1518 1490253356 : 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 1490253356 : fadvheat = c0
1542 :
1543 : ! calculate fluxes from base upwards
1544 10431773492 : do k = 1, nilyr-1
1545 :
1546 10431773492 : fadvheat = fadvheat - q(k) * (qbr(k+1) - qbr(k))
1547 :
1548 : enddo ! k
1549 :
1550 1490253356 : k = nilyr
1551 :
1552 1490253356 : fadvheat = fadvheat - q(k) * (qocn - qbr(k))
1553 :
1554 1490253356 : end subroutine picard_drainage_fluxes
1555 :
1556 : !=======================================================================
1557 :
1558 1490253356 : subroutine picard_flushing_fluxes(nilyr, &
1559 : fadvheat, w, &
1560 1490253356 : 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 1490253356 : fadvheat = fadvheat + w * (qbr(nilyr) - qpond)
1579 :
1580 1490253356 : end subroutine picard_flushing_fluxes
1581 :
1582 : !=======================================================================
1583 :
1584 1490253356 : subroutine maximum_variables_changes(lsnow, &
1585 : Tsf, Tsf_prev, dTsf, &
1586 1490253356 : zTsn, zTsn_prev, dzTsn, &
1587 1490253356 : 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 1490253356 : dTsf = abs(Tsf - Tsf_prev)
1610 :
1611 1490253356 : if (lsnow) then
1612 2535280288 : dzTsn = maxval(abs(zTsn - zTsn_prev))
1613 : else ! lsnow
1614 222613212 : dzTsn = c0
1615 : endif ! lsnow
1616 :
1617 11922026848 : dzTin = maxval(abs(zTin - zTin_prev))
1618 :
1619 1490253356 : end subroutine maximum_variables_changes
1620 :
1621 : !=======================================================================
1622 :
1623 2155800729 : subroutine total_energy_content(lsnow, &
1624 : nilyr, nslyr, &
1625 2155800729 : 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 2155800729 : energy = c0
1653 :
1654 2155800729 : if (lsnow) then
1655 :
1656 3668218700 : do k = 1, nslyr
1657 :
1658 3668218700 : energy = energy + hslyr * zqsn(k)
1659 :
1660 : enddo ! k
1661 :
1662 : endif ! lsnow
1663 :
1664 17246405832 : do k = 1, nilyr
1665 :
1666 17246405832 : energy = energy + hilyr * zqin(k)
1667 :
1668 : enddo ! k
1669 :
1670 2155800729 : end subroutine total_energy_content
1671 :
1672 : !=======================================================================
1673 :
1674 665547373 : subroutine picard_updates(nilyr, zTin, &
1675 1331094746 : 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 5324378984 : do k = 1, nilyr
1695 :
1696 4658831611 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
1697 5324378984 : qbr(k) = enthalpy_brine(zTin(k))
1698 :
1699 : enddo ! k
1700 :
1701 665547373 : end subroutine picard_updates
1702 :
1703 : !=======================================================================
1704 :
1705 1490253356 : 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 11922026848 : do k = 1, nilyr
1724 :
1725 11922026848 : qbr(k) = enthalpy_brine(zTin(k))
1726 :
1727 : enddo ! k
1728 :
1729 1490253356 : end subroutine picard_updates_enthalpy
1730 :
1731 : !=======================================================================
1732 :
1733 1490253356 : subroutine picard_final(lsnow, &
1734 : nilyr, nslyr, &
1735 2980506712 : zqin, zqsn, &
1736 2980506712 : zTin, zTsn, &
1737 1490253356 : 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 11922026848 : do k = 1, nilyr
1761 11922026848 : zqin(k) = enthalpy_mush_liquid_fraction(zTin(k), phi(k))
1762 : enddo ! k
1763 :
1764 1490253356 : if (lsnow) then
1765 :
1766 2535280288 : do k = 1, nslyr
1767 2535280288 : zqsn(k) = enthalpy_snow(zTsn(k))
1768 : enddo ! k
1769 :
1770 : endif ! lsnow
1771 :
1772 1490253356 : end subroutine picard_final
1773 :
1774 : !=======================================================================
1775 :
1776 665547373 : 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 665547373 : if (lsnow) then
1798 :
1799 566469206 : dxp(1) = hslyr / c2
1800 :
1801 566469206 : do l = 2, nslyr
1802 :
1803 566469206 : dxp(l) = hslyr
1804 :
1805 : enddo ! l
1806 :
1807 566469206 : dxp(nslyr+1) = (hilyr + hslyr) / c2
1808 :
1809 3965284442 : do l = nslyr+2, nilyr+nslyr
1810 :
1811 3965284442 : dxp(l) = hilyr
1812 :
1813 : enddo ! l
1814 :
1815 566469206 : dxp(nilyr+nslyr+1) = hilyr / c2
1816 :
1817 : else ! lsnow
1818 :
1819 99078167 : dxp(1) = hilyr / c2
1820 :
1821 693547169 : do l = 2, nilyr
1822 :
1823 693547169 : dxp(l) = hilyr
1824 :
1825 : enddo ! l
1826 :
1827 99078167 : dxp(nilyr+1) = hilyr / c2
1828 :
1829 198156334 : do l = nilyr+2, nilyr+nslyr+1
1830 :
1831 198156334 : dxp(l) = c0
1832 :
1833 : enddo ! l
1834 :
1835 : endif ! lsnow
1836 :
1837 665547373 : end subroutine calc_intercell_thickness
1838 :
1839 : !=======================================================================
1840 :
1841 665547373 : subroutine calc_intercell_conductivity(lsnow, &
1842 : nilyr, nslyr, &
1843 665547373 : km, ks, &
1844 : hilyr, hslyr, &
1845 665547373 : 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 36806816 : 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 665547373 : if (lsnow) then
1875 :
1876 566469206 : kcstar(1) = ks(1)
1877 :
1878 566469206 : do l = 2, nslyr
1879 :
1880 0 : k = l
1881 566469206 : kcstar(l) = (c2 * ks(k) * ks(k-1)) / (ks(k) + ks(k-1))
1882 :
1883 : enddo ! l
1884 :
1885 566469206 : fe = hilyr / (hilyr + hslyr)
1886 566469206 : kcstar(nslyr+1) = c1 / ((c1 - fe) / ks(nslyr) + fe / km(1))
1887 :
1888 3965284442 : do k = 2, nilyr
1889 :
1890 3398815236 : l = k + nslyr
1891 3965284442 : kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
1892 :
1893 : enddo ! k
1894 :
1895 566469206 : kcstar(nilyr+nslyr+1) = km(nilyr)
1896 :
1897 : else ! lsnow
1898 :
1899 99078167 : kcstar(1) = km(1)
1900 :
1901 693547169 : do k = 2, nilyr
1902 :
1903 594469002 : l = k
1904 693547169 : kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
1905 :
1906 : enddo ! k
1907 :
1908 99078167 : kcstar(nilyr+1) = km(nilyr)
1909 :
1910 198156334 : do l = nilyr+2, nilyr+nslyr+1
1911 :
1912 198156334 : kcstar(l) = c0
1913 :
1914 : enddo ! l
1915 :
1916 : endif ! lsnow
1917 :
1918 665547373 : end subroutine calc_intercell_conductivity
1919 :
1920 : !=======================================================================
1921 :
1922 1490253356 : subroutine solve_heat_conduction(lsnow, lcold, &
1923 : nilyr, nslyr, &
1924 : Tsf, Tbot, &
1925 1490253356 : zqin0, zqsn0, &
1926 1490253356 : phi, dt, &
1927 : qpond, qocn, &
1928 1490253356 : q, w, &
1929 : hilyr, hslyr, &
1930 1490253356 : dxp, kcstar, &
1931 1490253356 : Iswabs, Sswabs, &
1932 : fsurfn, dfsurfn_dTsf, &
1933 1490253356 : 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 3636544240 : Ap , & ! diagonal of tridiagonal matrix
1979 3636544240 : As , & ! lower off-diagonal of tridiagonal matrix
1980 3800553622 : An , & ! upper off-diagonal of tridiagonal matrix
1981 3636544240 : b , & ! right hand side of matrix solve
1982 2310300266 : 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 1490253356 : if (lsnow) then
1991 :
1992 1267640144 : 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 1159736259 : dt)
2006 1159736259 : 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 107903885 : dt)
2021 107903885 : if (icepack_warnings_aborted(subname)) return
2022 :
2023 : endif ! lcold
2024 :
2025 : else ! lsnow
2026 :
2027 222613212 : 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 149447219 : dt)
2041 149447219 : 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 73165993 : dt)
2056 73165993 : 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 1490253356 : An(1:nyn), Ap(1:nyn), As(1:nyn), b(1:nyn), T(1:nyn), nyn)
2065 1490253356 : if (icepack_warnings_aborted(subname)) return
2066 :
2067 : call update_temperatures(lsnow, lcold, &
2068 : nilyr, nslyr, &
2069 0 : T, Tsf, &
2070 1490253356 : zTin, zTsn)
2071 1490253356 : if (icepack_warnings_aborted(subname)) return
2072 :
2073 : end subroutine solve_heat_conduction
2074 :
2075 : !=======================================================================
2076 :
2077 1490253356 : subroutine update_temperatures(lsnow, lcold, &
2078 : nilyr, nslyr, &
2079 1490253356 : T, Tsf, &
2080 1490253356 : 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 1490253356 : if (lsnow) then
2107 :
2108 1267640144 : if (lcold) then
2109 :
2110 1159736259 : Tsf = T(1)
2111 :
2112 2319472518 : do k = 1, nslyr
2113 1159736259 : l = k + 1
2114 2319472518 : zTsn(k) = T(l)
2115 : enddo ! k
2116 :
2117 9277890072 : do k = 1, nilyr
2118 8118153813 : l = k + nslyr + 1
2119 9277890072 : zTin(k) = T(l)
2120 : enddo ! k
2121 :
2122 : else ! lcold
2123 :
2124 215807770 : do k = 1, nslyr
2125 107903885 : l = k
2126 215807770 : zTsn(k) = T(l)
2127 : enddo ! k
2128 :
2129 863231080 : do k = 1, nilyr
2130 755327195 : l = k + nslyr
2131 863231080 : zTin(k) = T(l)
2132 : enddo ! k
2133 :
2134 : endif ! lcold
2135 :
2136 : else ! lsnow
2137 :
2138 222613212 : if (lcold) then
2139 :
2140 149447219 : Tsf = T(1)
2141 :
2142 1195577752 : do k = 1, nilyr
2143 1046130533 : l = k + 1
2144 1195577752 : zTin(k) = T(l)
2145 : enddo ! k
2146 :
2147 : else ! lcold
2148 :
2149 585327944 : do k = 1, nilyr
2150 512161951 : l = k
2151 585327944 : zTin(k) = T(l)
2152 : enddo ! k
2153 :
2154 : endif ! lcold
2155 :
2156 : endif ! lsnow
2157 :
2158 1490253356 : end subroutine update_temperatures
2159 :
2160 : !=======================================================================
2161 :
2162 146331986 : subroutine matrix_elements_nosnow_melt(Ap, As, An, b, nyn, &
2163 : nilyr, &
2164 : Tsf, Tbot, &
2165 73165993 : zqin0, &
2166 : qpond, qocn, &
2167 73165993 : phi, q, &
2168 : w, &
2169 : hilyr, &
2170 73165993 : dxp, kcstar, &
2171 73165993 : 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 73165993 : k = 1
2215 73165993 : l = k
2216 :
2217 4403370 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2218 8806740 : kcstar(k+1) / dxp(k+1) + &
2219 4403370 : kcstar(k) / dxp(k) + &
2220 4403370 : q(k) * cp_ocn * rhow + &
2221 77569363 : w * cp_ocn * rhow
2222 13210110 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2223 81972733 : q(k) * cp_ocn * rhow
2224 73165993 : An(l) = c0
2225 4403370 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2226 4403370 : (kcstar(k) / dxp(k)) * Tsf + &
2227 73165993 : w * qpond
2228 :
2229 : ! interior ice layers
2230 438995958 : do k = 2, nilyr-1
2231 :
2232 365829965 : l = k
2233 :
2234 22016850 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2235 44033700 : kcstar(k+1) / dxp(k+1) + &
2236 22016850 : kcstar(k) / dxp(k) + &
2237 22016850 : q(k) * cp_ocn * rhow + &
2238 387846815 : w * cp_ocn * rhow
2239 66050550 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2240 409863665 : q(k) * cp_ocn * rhow
2241 22016850 : An(l) = -kcstar(k) / dxp(k) - &
2242 365829965 : w * cp_ocn * rhow
2243 438995958 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2244 :
2245 : enddo ! k
2246 :
2247 : ! bottom layer
2248 73165993 : k = nilyr
2249 73165993 : l = k
2250 :
2251 4403370 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2252 8806740 : kcstar(k+1) / dxp(k+1) + &
2253 4403370 : kcstar(k) / dxp(k) + &
2254 4403370 : q(k) * cp_ocn * rhow + &
2255 77569363 : w * cp_ocn * rhow
2256 73165993 : As(l) = c0
2257 4403370 : An(l) = -kcstar(k) / dxp(k) - &
2258 73165993 : w * cp_ocn * rhow
2259 4403370 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2260 8806740 : (kcstar(k+1) * Tbot) / dxp(k+1) + &
2261 81972733 : q(k) * qocn
2262 :
2263 73165993 : nyn = nilyr
2264 :
2265 73165993 : end subroutine matrix_elements_nosnow_melt
2266 :
2267 : !=======================================================================
2268 :
2269 298894438 : subroutine matrix_elements_nosnow_cold(Ap, As, An, b, nyn, &
2270 : nilyr, &
2271 : Tsf, Tbot, &
2272 149447219 : zqin0, &
2273 : qpond, qocn, &
2274 149447219 : phi, q, &
2275 : w, &
2276 : hilyr, &
2277 149447219 : dxp, kcstar, &
2278 149447219 : 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 149447219 : l = 1
2325 149447219 : Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
2326 149447219 : As(l) = kcstar(1) / dxp(1)
2327 149447219 : An(l) = c0
2328 149447219 : b (l) = dfsurfn_dTsf * Tsf - fsurfn
2329 :
2330 : ! surface layer
2331 149447219 : k = 1
2332 149447219 : l = k + 1
2333 :
2334 3722439 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2335 7444878 : kcstar(k+1) / dxp(k+1) + &
2336 3722439 : kcstar(k) / dxp(k) + &
2337 3722439 : q(k) * cp_ocn * rhow + &
2338 153169658 : w * cp_ocn * rhow
2339 11167317 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2340 156892097 : q(k) * cp_ocn * rhow
2341 149447219 : An(l) = -kcstar(k) / dxp(k)
2342 3722439 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2343 149447219 : w * qpond
2344 :
2345 : ! interior ice layers
2346 896683314 : do k = 2, nilyr-1
2347 :
2348 747236095 : l = k + 1
2349 :
2350 18612195 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2351 37224390 : kcstar(k+1) / dxp(k+1) + &
2352 18612195 : kcstar(k) / dxp(k) + &
2353 18612195 : q(k) * cp_ocn * rhow + &
2354 765848290 : w * cp_ocn * rhow
2355 55836585 : As(l) = -kcstar(k+1) / dxp(k+1) - &
2356 784460485 : q(k) * cp_ocn * rhow
2357 18612195 : An(l) = -kcstar(k) / dxp(k) - &
2358 747236095 : w * cp_ocn * rhow
2359 896683314 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2360 :
2361 : enddo ! k
2362 :
2363 : ! bottom layer
2364 149447219 : k = nilyr
2365 149447219 : l = k + 1
2366 :
2367 3722439 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2368 7444878 : kcstar(k+1) / dxp(k+1) + &
2369 3722439 : kcstar(k) / dxp(k) + &
2370 3722439 : q(k) * cp_ocn * rhow + &
2371 153169658 : w * cp_ocn * rhow
2372 149447219 : As(l) = c0
2373 3722439 : An(l) = -kcstar(k) / dxp(k) - &
2374 149447219 : w * cp_ocn * rhow
2375 3722439 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2376 7444878 : (kcstar(k+1) * Tbot) / dxp(k+1) + &
2377 156892097 : q(k) * qocn
2378 :
2379 149447219 : nyn = nilyr + 1
2380 :
2381 149447219 : end subroutine matrix_elements_nosnow_cold
2382 :
2383 : !=======================================================================
2384 :
2385 215807770 : subroutine matrix_elements_snow_melt(Ap, As, An, b, nyn, &
2386 : nilyr, nslyr, &
2387 : Tsf, Tbot, &
2388 215807770 : zqin0, zqsn0, &
2389 : qpond, qocn, &
2390 107903885 : phi, q, &
2391 : w, &
2392 : hilyr, hslyr, &
2393 107903885 : dxp, kcstar, &
2394 215807770 : 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 107903885 : k = 1
2442 107903885 : l = k
2443 :
2444 9898576 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2445 19797152 : kcstar(l+1) / dxp(l+1) + &
2446 127701037 : kcstar(l) / dxp(l)
2447 107903885 : As(l) = -kcstar(l+1) / dxp(l+1)
2448 107903885 : An(l) = c0
2449 9898576 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k) + &
2450 107903885 : (kcstar(l) * Tsf) / dxp(l)
2451 :
2452 : ! interior snow layers
2453 107903885 : 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 107903885 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2463 :
2464 : enddo ! k
2465 :
2466 : ! top ice layer
2467 107903885 : k = 1
2468 107903885 : l = nslyr + k
2469 :
2470 9898576 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2471 19797152 : kcstar(l+1) / dxp(l+1) + &
2472 9898576 : kcstar(l) / dxp(l) + &
2473 9898576 : q(k) * cp_ocn * rhow + &
2474 117802461 : w * cp_ocn * rhow
2475 29695728 : As(l) = -kcstar(l+1) / dxp(l+1) - &
2476 127701037 : q(k) * cp_ocn * rhow
2477 107903885 : An(l) = -kcstar(l) / dxp(l)
2478 9898576 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2479 107903885 : w * qpond
2480 :
2481 : ! interior ice layers
2482 647423310 : do k = 2, nilyr-1
2483 :
2484 539519425 : l = nslyr + k
2485 :
2486 49492880 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2487 98985760 : kcstar(l+1) / dxp(l+1) + &
2488 49492880 : kcstar(l) / dxp(l) + &
2489 49492880 : q(k) * cp_ocn * rhow + &
2490 589012305 : w * cp_ocn * rhow
2491 148478640 : As(l) = -kcstar(l+1) / dxp(l+1) - &
2492 638505185 : q(k) * cp_ocn * rhow
2493 49492880 : An(l) = -kcstar(l) / dxp(l) - &
2494 539519425 : w * cp_ocn * rhow
2495 647423310 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2496 :
2497 : enddo ! k
2498 :
2499 : ! bottom layer
2500 107903885 : k = nilyr
2501 107903885 : l = nilyr + nslyr
2502 :
2503 9898576 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2504 19797152 : kcstar(l+1) / dxp(l+1) + &
2505 9898576 : kcstar(l) / dxp(l) + &
2506 9898576 : q(k) * cp_ocn * rhow + &
2507 117802461 : w * cp_ocn * rhow
2508 107903885 : As(l) = c0
2509 9898576 : An(l) = -kcstar(l) / dxp(l) - &
2510 107903885 : w * cp_ocn * rhow
2511 9898576 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2512 19797152 : (kcstar(l+1) * Tbot) / dxp(l+1) + &
2513 127701037 : q(k) * qocn
2514 :
2515 107903885 : nyn = nilyr + nslyr
2516 :
2517 107903885 : end subroutine matrix_elements_snow_melt
2518 :
2519 : !=======================================================================
2520 :
2521 2319472518 : subroutine matrix_elements_snow_cold(Ap, As, An, b, nyn, &
2522 : nilyr, nslyr, &
2523 : Tsf, Tbot, &
2524 2319472518 : zqin0, zqsn0, &
2525 : qpond, qocn, &
2526 1159736259 : phi, q, &
2527 : w, &
2528 : hilyr, hslyr, &
2529 1159736259 : dxp, kcstar, &
2530 2319472518 : 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 1159736259 : l = 1
2582 1159736259 : Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
2583 1159736259 : As(l) = kcstar(1) / dxp(1)
2584 1159736259 : An(l) = c0
2585 1159736259 : b (l) = dfsurfn_dTsf * Tsf - fsurfn
2586 :
2587 : ! surface layer
2588 1159736259 : k = 1
2589 1159736259 : l = k + 1
2590 1159736259 : m = 1
2591 :
2592 63980306 : Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
2593 127960612 : kcstar(m+1) / dxp(m+1) + &
2594 1287696871 : kcstar(m) / dxp(m)
2595 1159736259 : As(l) = -kcstar(m+1) / dxp(m+1)
2596 1159736259 : An(l) = -kcstar(m) / dxp(m)
2597 1159736259 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2598 :
2599 : ! interior snow layers
2600 1159736259 : 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 1159736259 : b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
2611 :
2612 : enddo ! k
2613 :
2614 : ! top ice layer
2615 1159736259 : k = 1
2616 1159736259 : l = nslyr + k + 1
2617 1159736259 : m = k + nslyr
2618 :
2619 63980306 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2620 127960612 : kcstar(m+1) / dxp(m+1) + &
2621 63980306 : kcstar(m) / dxp(m) + &
2622 63980306 : q(k) * cp_ocn * rhow + &
2623 1223716565 : w * cp_ocn * rhow
2624 191940918 : As(l) = -kcstar(m+1) / dxp(m+1) - &
2625 1287696871 : q(k) * cp_ocn * rhow
2626 1159736259 : An(l) = -kcstar(m) / dxp(m)
2627 63980306 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2628 1159736259 : w * qpond
2629 :
2630 : ! interior ice layers
2631 6958417554 : do k = 2, nilyr-1
2632 :
2633 5798681295 : l = nslyr + k + 1
2634 5798681295 : m = k + nslyr
2635 :
2636 319901530 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2637 639803060 : kcstar(m+1) / dxp(m+1) + &
2638 319901530 : kcstar(m) / dxp(m) + &
2639 319901530 : q(k) * cp_ocn * rhow + &
2640 6118582825 : w * cp_ocn * rhow
2641 959704590 : As(l) = -kcstar(m+1) / dxp(m+1) - &
2642 6438484355 : q(k) * cp_ocn * rhow
2643 319901530 : An(l) = -kcstar(m) / dxp(m) - &
2644 5798681295 : w * cp_ocn * rhow
2645 6958417554 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
2646 :
2647 : enddo ! k
2648 :
2649 : ! bottom layer
2650 1159736259 : k = nilyr
2651 1159736259 : l = nilyr + nslyr + 1
2652 1159736259 : m = k + nslyr
2653 :
2654 63980306 : Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
2655 127960612 : kcstar(m+1) / dxp(m+1) + &
2656 63980306 : kcstar(m) / dxp(m) + &
2657 63980306 : q(k) * cp_ocn * rhow + &
2658 1223716565 : w * cp_ocn * rhow
2659 1159736259 : As(l) = c0
2660 63980306 : An(l) = -kcstar(m) / dxp(m) - &
2661 1159736259 : w * cp_ocn * rhow
2662 63980306 : b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
2663 127960612 : (kcstar(m+1) * Tbot) / dxp(m+1) + &
2664 1287696871 : q(k) * qocn
2665 :
2666 1159736259 : nyn = nilyr + nslyr + 1
2667 :
2668 1159736259 : end subroutine matrix_elements_snow_cold
2669 :
2670 : !=======================================================================
2671 :
2672 1331094746 : subroutine solve_salinity(zSin, Sbr, &
2673 : Spond, sss, &
2674 1331094746 : 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 960001901 : zSin0
2706 :
2707 : character(len=*),parameter :: subname='(solve_salinity)'
2708 :
2709 5324378984 : zSin0 = zSin
2710 :
2711 665547373 : k = 1
2712 36806816 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2713 73613632 : ((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + &
2714 36806816 : dSdt(k) + &
2715 702354189 : (w * (Spond - Sbr(k))) / hilyr) * dt)
2716 :
2717 3993284238 : do k = 2, nilyr-1
2718 :
2719 184034080 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2720 368068160 : ((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + &
2721 184034080 : dSdt(k) + &
2722 4177318318 : (w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
2723 :
2724 : enddo ! k
2725 :
2726 665547373 : k = nilyr
2727 36806816 : zSin(k) = zSin(k) + max(S_min - zSin(k), &
2728 36806816 : ((q(k) * (sss - Sbr(k))) / hilyr + &
2729 36806816 : dSdt(k) + &
2730 665547373 : (w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
2731 :
2732 :
2733 5324378984 : 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 1490253356 : 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 3636544240 : cp , & ! modified upper off-diagonal vector
2779 3636544240 : 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 1490253356 : cp(1) = c(1) / b(1)
2788 11518343758 : do i = 2, n-1
2789 11518343758 : cp(i) = c(i) / (b(i) - cp(i-1)*a(i))
2790 : enddo
2791 :
2792 1490253356 : dp(1) = d(1) / b(1)
2793 13008597114 : do i = 2, n
2794 13008597114 : dp(i) = (d(i) - dp(i-1)*a(i)) / (b(i) - cp(i-1)*a(i))
2795 : enddo
2796 :
2797 : ! back substitution
2798 1490253356 : x(n) = dp(n)
2799 13008597114 : do i = n-1,1,-1
2800 13008597114 : x(i) = dp(i) - cp(i)*x(i+1)
2801 : enddo
2802 :
2803 1490253356 : end subroutine tdma_solve_sparse
2804 :
2805 : !=======================================================================
2806 : ! Effect of salinity
2807 : !=======================================================================
2808 :
2809 8493988986 : 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 8493988986 : perm = 3.0e-8_dbl_kind * max(phi - phic, c0)**3
2826 :
2827 8493988986 : end function permeability
2828 :
2829 : !=======================================================================
2830 :
2831 609317819 : subroutine explicit_flow_velocities(nilyr, zSin, &
2832 609317819 : zTin, Tsf, &
2833 609317819 : Tbot, q, &
2834 1218635638 : dSdt, Sbr, &
2835 609317819 : 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 1413344656 : phi ! ice layer liquid fraction
2876 :
2877 : real(kind=dbl_kind), dimension(0:nilyr) :: &
2878 1445796159 : rho ! ice layer brine density (kg m-3)
2879 :
2880 : real(kind=dbl_kind) :: &
2881 32451503 : rho_ocn , & ! ocean density (kg m-3)
2882 32451503 : perm_min , & ! minimum permeability from layer to ocean (m2)
2883 32451503 : perm_harm , & ! harmonic mean of permeability from layer to ocean (m2)
2884 32451503 : rho_sum , & ! sum of the brine densities from layer to ocean (kg m-3)
2885 32451503 : rho_pipe , & ! density of the brine in the channel (kg m-3)
2886 32451503 : z , & ! distance to layer from top surface (m)
2887 32451503 : perm , & ! ice layer permeability (m2)
2888 32451503 : drho , & ! brine density difference between layer and ocean (kg m-3)
2889 32451503 : Ra , & ! local mush Rayleigh number
2890 32451503 : rn , & ! real value of number of layers considered
2891 32451503 : L , & ! thickness of the layers considered (m)
2892 32451503 : dx , & ! horizontal size of convective flow (m)
2893 32451503 : dx2 , & ! square of the horizontal size of convective flow (m2)
2894 32451503 : Am , & ! A parameter for mush
2895 32451503 : Bm , & ! B parameter for mush
2896 32451503 : Ap , & ! A parameter for channel
2897 32451503 : Bp , & ! B parameter for channel
2898 32451503 : qlimit , & ! limit to vertical Darcy flow for numerical stability
2899 32451503 : dS_guess , & ! expected bulk salinity without limits
2900 32451503 : alpha , & ! desalination limiting factor
2901 32451503 : 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 4874542552 : do k = 1, nilyr
2910 :
2911 4265224733 : Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
2912 4265224733 : phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
2913 4265224733 : qbr(k) = enthalpy_brine(zTin(k))
2914 4874542552 : rho(k) = icepack_mushy_density_brine(Sbr(k))
2915 :
2916 : enddo ! k
2917 :
2918 609317819 : rho(0) = rho(1)
2919 :
2920 : ! ocean conditions
2921 609317819 : Sbr(nilyr+1) = sss
2922 609317819 : qbr(nilyr+1) = qocn
2923 609317819 : rho_ocn = icepack_mushy_density_brine(sss)
2924 :
2925 : ! initialize accumulated quantities
2926 609317819 : perm_min = bignum
2927 609317819 : perm_harm = c0
2928 609317819 : rho_sum = c0
2929 :
2930 : ! limit to q for numerical stability
2931 609317819 : qlimit = (fracmax * hilyr) / dt
2932 :
2933 : ! no flow through ice top surface
2934 609317819 : q(0) = c0
2935 :
2936 609317819 : ra_constants = gravit / (viscosity_dyn * kappal)
2937 : ! first iterate over layers going up
2938 4874542552 : do k = nilyr, 1, -1
2939 :
2940 : ! vertical position from ice top surface
2941 4265224733 : z = ((real(k, dbl_kind) - p5) / real(nilyr, dbl_kind)) * hin
2942 :
2943 : ! permeabilities
2944 4265224733 : perm = permeability(phi(k))
2945 4265224733 : perm_min = min(perm_min,perm)
2946 4265224733 : perm_harm = perm_harm + (c1 / max(perm,1.0e-30_dbl_kind))
2947 :
2948 : ! densities
2949 4265224733 : rho_sum = rho_sum + rho(k)
2950 : !rho_pipe = rho(k)
2951 4265224733 : rho_pipe = p5 * (rho(k) + rho(k-1))
2952 4265224733 : drho = max(rho(k) - rho_ocn, c0)
2953 :
2954 : ! mush Rayleigh number
2955 4265224733 : Ra = drho * (hin-z) * perm_min * ra_constants
2956 :
2957 : ! height of mush layer to layer k
2958 4265224733 : rn = real(nilyr-k+1,dbl_kind)
2959 4265224733 : L = rn * hilyr
2960 :
2961 : ! horizontal size of convection
2962 4265224733 : dx = L * c2 * aspect_rapid_mode
2963 4265224733 : dx2 = dx**2
2964 :
2965 : ! determine vertical Darcy flow
2966 4265224733 : Am = (dx2 * rn) / (viscosity_dyn * perm_harm)
2967 4265224733 : Bm = (-gravit * rho_sum) / rn
2968 :
2969 4265224733 : Ap = (pi * a_rapid_mode**4) / (c8 * viscosity_dyn)
2970 4265224733 : Bp = -rho_pipe * gravit
2971 :
2972 4265224733 : 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 4265224733 : q(k) = min(q(k) * (max(Ra - Rac_rapid_mode, c0) / (Ra+puny)), qlimit)
2976 :
2977 : ! late stage drainage
2978 227160521 : dSdt(k) = dSdt_slow_mode * (max((zSin(k) - phi_c_slow_mode*Sbr(k)), c0) &
2979 4265224733 : * max((Tbot - Tsf), c0)) / (hin + 0.001_dbl_kind)
2980 :
2981 4265224733 : dSdt(k) = max(dSdt(k), (-zSin(k) * 0.5_dbl_kind) / dt)
2982 :
2983 : ! restrict flows to prevent too much salt loss
2984 4265224733 : dS_guess = (((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + dSdt(k)) * dt) * safety_factor
2985 :
2986 4265224733 : if (abs(dS_guess) < puny) then
2987 1863561863 : alpha = c1
2988 : else
2989 2401662870 : alpha = (zSin_min - zSin(k)) / dS_guess
2990 : endif
2991 :
2992 4265224733 : if (alpha < c0 .or. alpha > c1) alpha = c1
2993 :
2994 4265224733 : q(k) = q(k) * alpha
2995 4874542552 : dSdt(k) = dSdt(k) * alpha
2996 :
2997 : enddo ! k
2998 :
2999 609317819 : end subroutine explicit_flow_velocities
3000 :
3001 : !=======================================================================
3002 : ! Flushing
3003 : !=======================================================================
3004 :
3005 609317819 : subroutine flushing_velocity(zTin, &
3006 609317819 : 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 32451503 : perm , & ! ice layer permeability (m2)
3038 32451503 : ice_mass , & ! mass of ice (kg m-2)
3039 32451503 : perm_harm , & ! harmonic mean of ice permeability (m2)
3040 32451503 : hocn , & ! ocean surface height above ice base (m)
3041 32451503 : hbrine , & ! brine surface height above ice base (m)
3042 32451503 : w_down_max , & ! maximum downward flushing Darcy flow rate (m s-1)
3043 32451503 : phi_min , & ! minimum porosity in the mush
3044 32451503 : wlimit , & ! limit to w to avoid advecting all brine in layer
3045 32451503 : 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 609317819 : w = c0
3054 :
3055 : ! only flush if ponds are active
3056 609317819 : if (tr_pond) then
3057 :
3058 604109179 : ice_mass = c0
3059 604109179 : perm_harm = c0
3060 604109179 : phi_min = c1
3061 :
3062 4832873432 : do k = 1, nilyr
3063 :
3064 : ! liquid fraction
3065 : !phi = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
3066 4228764253 : phi_min = min(phi_min,phi(k))
3067 :
3068 : ! permeability
3069 4228764253 : perm = permeability(phi(k))
3070 :
3071 : ! ice mass
3072 222602961 : ice_mass = ice_mass + phi(k) * icepack_mushy_density_brine(liquidus_brine_salinity_mush(zTin(k))) + &
3073 4228764253 : (c1 - phi(k)) * rhoi
3074 :
3075 : ! permeability harmonic mean
3076 4832873432 : perm_harm = perm_harm + c1 / (perm + 1e-30_dbl_kind)
3077 :
3078 : enddo ! k
3079 :
3080 604109179 : ice_mass = ice_mass * hilyr
3081 :
3082 604109179 : perm_harm = real(nilyr,dbl_kind) / perm_harm
3083 :
3084 : ! calculate ocean surface height above bottom of ice
3085 604109179 : hocn = (ice_mass + hpond * apond * rhow + hsn * rhos) / rhow
3086 :
3087 : ! calculate brine height above bottom of ice
3088 604109179 : hbrine = hin + hpond
3089 :
3090 : ! pressure head
3091 604109179 : dhhead = max(hbrine - hocn,c0)
3092 :
3093 : ! darcy flow through ice
3094 604109179 : w = (perm_harm * rhow * gravit * (dhhead / hin)) / viscosity_dyn
3095 :
3096 : ! maximum down flow to drain pond
3097 604109179 : w_down_max = (hpond * apond) / dt
3098 :
3099 : ! limit flow
3100 604109179 : w = min(w,w_down_max)
3101 :
3102 : ! limit amount of brine that can be advected out of any particular layer
3103 604109179 : wlimit = (advection_limit * phi_min * hilyr) / dt
3104 :
3105 604109179 : if (abs(w) > puny) then
3106 127192674 : w = w * max(min(abs(wlimit/w),c1),c0)
3107 : else
3108 476916505 : w = c0
3109 : endif
3110 :
3111 604109179 : w = max(w, c0)
3112 :
3113 : endif
3114 :
3115 609317819 : end subroutine flushing_velocity
3116 :
3117 : !=======================================================================
3118 :
3119 609317819 : 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 609317819 : if (tr_pond) then
3138 604109179 : if (apond > c0 .and. hpond > c0) then
3139 :
3140 : ! flush pond through mush
3141 386503755 : hpond = hpond - w * dt / apond
3142 :
3143 386503755 : hpond = max(hpond, c0)
3144 :
3145 : ! exponential decay of pond
3146 386503755 : hpond = hpond - lambda_pond * dt * (hpond + hpond0)
3147 :
3148 386503755 : hpond = max(hpond, c0)
3149 :
3150 : endif
3151 : endif
3152 :
3153 609317819 : end subroutine flush_pond
3154 :
3155 : !=======================================================================
3156 :
3157 609317819 : subroutine flood_ice(hsn, hin, &
3158 : nslyr, nilyr, &
3159 : hslyr, hilyr, &
3160 609317819 : zqsn, zqin, &
3161 609317819 : phi, dt, &
3162 1218635638 : 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 32451503 : hin2 , & ! new ice thickness (m)
3201 32451503 : hsn2 , & ! new snow thickness (m)
3202 32451503 : hilyr2 , & ! new ice layer thickness (m)
3203 32451503 : hslyr2 , & ! new snow layer thickness (m)
3204 32451503 : dh , & ! thickness of snowice formation (m)
3205 32451503 : phi_snowice , & ! liquid fraction of new snow ice
3206 32451503 : rho_snowice , & ! density of snowice (kg m-3)
3207 32451503 : zSin_snowice , & ! bulk salinity of new snowice (ppt)
3208 32451503 : zqin_snowice , & ! ice enthalpy of new snowice (J m-2)
3209 32451503 : zqsn_snowice , & ! snow enthalpy of snow thats becoming snowice (J m-2)
3210 32451503 : freeboard_density , & ! negative of ice surface freeboard times the ocean density (kg m-2)
3211 32451503 : ice_mass , & ! mass of the ice (kg m-2)
3212 32451503 : rho_ocn , & ! density of the ocean (kg m-3)
3213 32451503 : ice_density , & ! density of ice layer (kg m-3)
3214 32451503 : hadded , & ! thickness rate of water used from ocean (m/s)
3215 32451503 : wadded , & ! mass rate of water used from ocean (kg/m^2/s)
3216 32451503 : 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 609317819 : snoice = c0
3227 :
3228 : ! check we have snow
3229 609317819 : if (hsn > puny) then
3230 :
3231 531599226 : rho_ocn = icepack_mushy_density_brine(sss)
3232 :
3233 : ! ice mass
3234 531599226 : ice_mass = c0
3235 4252793808 : do k = 1, nilyr
3236 3721194582 : ice_density = min(phi(k) * icepack_mushy_density_brine(Sbr(k)) + (c1 - phi(k)) * rhoi,rho_ocn)
3237 4252793808 : ice_mass = ice_mass + ice_density
3238 : enddo ! k
3239 531599226 : ice_mass = ice_mass * hilyr
3240 :
3241 : ! negative freeboard times ocean density
3242 531599226 : freeboard_density = max(ice_mass + hsn * rhos - hin * rho_ocn, c0)
3243 :
3244 : ! check if have flooded ice
3245 531599226 : if (freeboard_density > c0) then
3246 :
3247 : ! sea ice fraction of newly formed snow ice
3248 51963859 : phi_snowice = (c1 - rhos / rhoi)
3249 :
3250 : ! density of newly formed snowice
3251 51963859 : rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
3252 :
3253 : ! calculate thickness of new ice added
3254 51963859 : dh = freeboard_density / (rho_ocn - rho_snowice + rhos)
3255 51963859 : dh = max(min(dh,hsn),c0)
3256 :
3257 : ! enthalpy of snow that becomes snowice
3258 51963859 : call enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
3259 51963859 : if (icepack_warnings_aborted(subname)) return
3260 :
3261 : ! change thicknesses
3262 51963859 : hin2 = hin + dh
3263 51963859 : hsn2 = hsn - dh
3264 :
3265 51963859 : hilyr2 = hin2 / real(nilyr,dbl_kind)
3266 51963859 : hslyr2 = hsn2 / real(nslyr,dbl_kind)
3267 :
3268 : ! properties of new snow ice
3269 51963859 : zSin_snowice = phi_snowice * sss
3270 51963859 : zqin_snowice = phi_snowice * qocn + zqsn_snowice
3271 :
3272 : ! change snow properties
3273 51963859 : call update_vertical_tracers_snow(nslyr, zqsn, hslyr, hslyr2)
3274 51963859 : if (icepack_warnings_aborted(subname)) return
3275 :
3276 : ! change ice properties
3277 0 : call update_vertical_tracers_ice(nilyr, zqin, hilyr, hilyr2, &
3278 51963859 : hin, hin2, zqin_snowice)
3279 51963859 : if (icepack_warnings_aborted(subname)) return
3280 0 : call update_vertical_tracers_ice(nilyr, zSin, hilyr, hilyr2, &
3281 51963859 : hin, hin2, zSin_snowice)
3282 51963859 : if (icepack_warnings_aborted(subname)) return
3283 0 : call update_vertical_tracers_ice(nilyr, phi, hilyr, hilyr2, &
3284 51963859 : hin, hin2, phi_snowice)
3285 51963859 : if (icepack_warnings_aborted(subname)) return
3286 :
3287 : ! change thicknesses
3288 51963859 : hilyr = hilyr2
3289 51963859 : hslyr = hslyr2
3290 51963859 : snoice = dh
3291 :
3292 51963859 : hadded = (dh * phi_snowice) / dt
3293 51963859 : wadded = hadded * rhoi
3294 51963859 : eadded = hadded * qocn
3295 : ! sadded = wadded * ice_ref_salinity * p001
3296 :
3297 : ! conservation
3298 51963859 : fadvheat = fadvheat - eadded
3299 :
3300 : endif
3301 :
3302 : endif
3303 :
3304 : end subroutine flood_ice
3305 :
3306 : !=======================================================================
3307 :
3308 51963859 : 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 4452112 : 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 51963859 : zqsn_snowice = c0
3335 :
3336 : ! snow depth and snow layers affected by snowice formation
3337 51963859 : if (hsn > puny) then
3338 51963859 : rnlyr = (dh / hsn) * nslyr
3339 51963859 : 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 51963859 : do k = nslyr, nslyr-nlyr+1, -1
3344 51963859 : zqsn_snowice = zqsn_snowice + zqsn(k) / rnlyr
3345 : enddo ! k
3346 :
3347 : ! partially converted snow layer
3348 : zqsn_snowice = zqsn_snowice + &
3349 51963859 : ((rnlyr - real(nlyr,dbl_kind)) / rnlyr) * zqsn(nslyr-nlyr)
3350 : endif
3351 :
3352 51963859 : end subroutine enthalpy_snow_snowice
3353 :
3354 : !=======================================================================
3355 :
3356 51963859 : 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 103927718 : 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 4452112 : z1a , & ! lower boundary of old cell
3380 4452112 : z1b , & ! upper boundary of old cell
3381 4452112 : z2a , & ! lower boundary of new cell
3382 4452112 : z2b , & ! upper boundary of new cell
3383 4452112 : 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 103927718 : do k2 = 1, nslyr
3389 :
3390 : ! initialize new tracer
3391 51963859 : trc2(k2) = c0
3392 :
3393 : ! calculate upper and lower boundary of new cell
3394 51963859 : z2a = (k2 - 1) * hlyr2
3395 51963859 : z2b = k2 * hlyr2
3396 :
3397 : ! loop over old grid cells
3398 103927718 : do k1 = 1, nslyr
3399 :
3400 : ! calculate upper and lower boundary of old cell
3401 51963859 : z1a = (k1 - 1) * hlyr1
3402 51963859 : z1b = k1 * hlyr1
3403 :
3404 : ! calculate overlap between old and new cell
3405 51963859 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3406 :
3407 : ! aggregate old grid cell contribution to new cell
3408 103927718 : trc2(k2) = trc2(k2) + overlap * trc(k1)
3409 :
3410 : enddo ! k1
3411 :
3412 : ! renormalize new grid cell
3413 103927718 : trc2(k2) = trc2(k2) / hlyr2
3414 :
3415 : enddo ! k2
3416 :
3417 : ! update vertical tracer array with the adjusted tracer
3418 103927718 : trc = trc2
3419 :
3420 51963859 : end subroutine update_vertical_tracers_snow
3421 :
3422 : !=======================================================================
3423 :
3424 155891577 : 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 391921170 : 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 13356336 : z1a , & ! lower boundary of old cell
3451 13356336 : z1b , & ! upper boundary of old cell
3452 13356336 : z2a , & ! lower boundary of new cell
3453 13356336 : z2b , & ! upper boundary of new cell
3454 13356336 : 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 1247132616 : do k2 = 1, nilyr
3460 :
3461 : ! initialize new tracer
3462 1091241039 : trc2(k2) = c0
3463 :
3464 : ! calculate upper and lower boundary of new cell
3465 1091241039 : z2a = (k2 - 1) * hlyr2
3466 1091241039 : z2b = k2 * hlyr2
3467 :
3468 : ! calculate upper and lower boundary of added snow ice at top
3469 1091241039 : z1a = c0
3470 1091241039 : z1b = h2 - h1
3471 :
3472 : ! calculate overlap between added ice and new cell
3473 1091241039 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3474 :
3475 : ! aggregate added ice contribution to new cell
3476 1091241039 : trc2(k2) = trc2(k2) + overlap * trc0
3477 :
3478 : ! loop over old grid cells
3479 8729928312 : do k1 = 1, nilyr
3480 :
3481 : ! calculate upper and lower boundary of old cell
3482 7638687273 : z1a = (k1 - 1) * hlyr1 + h2 - h1
3483 7638687273 : z1b = k1 * hlyr1 + h2 - h1
3484 :
3485 : ! calculate overlap between old and new cell
3486 7638687273 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
3487 :
3488 : ! aggregate old grid cell contribution to new cell
3489 8729928312 : trc2(k2) = trc2(k2) + overlap * trc(k1)
3490 :
3491 : enddo ! k1
3492 :
3493 : ! renormalize new grid cell
3494 1247132616 : trc2(k2) = trc2(k2) / hlyr2
3495 :
3496 : enddo ! k2
3497 :
3498 : ! update vertical tracer array with the adjusted tracer
3499 1247132616 : trc = trc2
3500 :
3501 155891577 : end subroutine update_vertical_tracers_ice
3502 :
3503 : !=======================================================================
3504 :
3505 : end module icepack_therm_mushy
3506 :
3507 : !=======================================================================
|
