Line data Source code
1 : !=======================================================================
2 : !
3 : ! Computes ice microstructural information for use in biogeochemistry
4 : !
5 : ! authors: Nicole Jeffery, LANL
6 : !
7 : module icepack_brine
8 :
9 : use icepack_kinds
10 : use icepack_parameters, only: p01, p001, p5, c0, c1, c2, c1p5, puny
11 : use icepack_parameters, only: gravit, rhoi, rhow, rhos, depressT
12 : use icepack_parameters, only: salt_loss, min_salin, rhosi
13 : use icepack_parameters, only: dts_b, l_sk
14 : use icepack_tracers, only: ntrcr, nt_qice, nt_sice, nt_bgc_S
15 : use icepack_tracers, only: nt_Tsfc
16 : use icepack_zbgc_shared, only: k_o, exp_h, Dm, Ra_c, viscos_dynamic, thinS
17 : use icepack_zbgc_shared, only: remap_zbgc
18 : use icepack_warnings, only: warnstr, icepack_warnings_add
19 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
20 :
21 : use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction
22 : use icepack_therm_shared, only: calculate_Tin_from_qin
23 :
24 : implicit none
25 :
26 : private
27 : public :: preflushing_changes, &
28 : compute_microS_mushy, &
29 : update_hbrine, &
30 : compute_microS, &
31 : calculate_drho, &
32 : icepack_init_hbrine, &
33 : icepack_init_zsalinity
34 :
35 : real (kind=dbl_kind), parameter :: &
36 : maxhbr = 1.25_dbl_kind , & ! brine overflows if hbr > maxhbr*hin
37 : viscos = 2.1e-6_dbl_kind, & ! kinematic viscosity (m^2/s)
38 : ! Brine salinity as a cubic function of temperature
39 : a1 = -21.4_dbl_kind , & ! (psu/C)
40 : a2 = -0.886_dbl_kind, & ! (psu/C^2)
41 : a3 = -0.012_dbl_kind, & ! (psu/C^3)
42 : ! Brine density as a quadratic of brine salinity
43 : b1 = 1000.0_dbl_kind, & ! (kg/m^3)
44 : b2 = 0.8_dbl_kind ! (kg/m^3/ppt)
45 :
46 : real (kind=dbl_kind), parameter :: &
47 : exp_argmax = 30.0_dbl_kind ! maximum argument of exponential for underflow
48 :
49 : !=======================================================================
50 :
51 : contains
52 :
53 : !=======================================================================
54 : ! Computes the top and bottom brine boundary changes for flushing
55 : ! works for zsalinity and tr_salinity
56 : !
57 : ! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice with
58 : ! dynamic salinity or the height ratio = hbr/vicen*aicen, where hbr is the
59 : ! height of the brine surface relative to the bottom of the ice. This volume fraction
60 : ! may be > 1 in which case there is brine above the ice surface (meltponds).
61 :
62 207881 : subroutine preflushing_changes (aicen, vicen, vsnon, &
63 : meltb, meltt, congel, &
64 : snoice, hice_old, dhice, &
65 : fbri, dhbr_top, dhbr_bot, &
66 : hbr_old, hin,hsn, firstice )
67 :
68 : real (kind=dbl_kind), intent(in) :: &
69 : aicen , & ! concentration of ice
70 : vicen , & ! volume per unit area of ice (m)
71 : vsnon , & ! volume per unit area of snow (m)
72 : meltb , & ! bottom ice melt (m)
73 : meltt , & ! top ice melt (m)
74 : congel , & ! bottom ice growth (m)
75 : snoice ! top ice growth from flooding (m)
76 :
77 : real (kind=dbl_kind), intent(out) :: &
78 : hbr_old ! old brine height (m)
79 :
80 : real (kind=dbl_kind), intent(inout) :: &
81 : hin , & ! ice thickness (m)
82 : hsn , & ! snow thickness (m)
83 : dhice ! change due to sublimation (<0)/condensation (>0) (m)
84 :
85 : real (kind=dbl_kind), intent(inout) :: &
86 : fbri , & ! trcrn(nt_fbri)
87 : dhbr_top , & ! brine change in top for diagnostics (m)
88 : dhbr_bot , & ! brine change in bottom for diagnostics (m)
89 : hice_old ! old ice thickness (m)
90 :
91 : logical (kind=log_kind), intent(in) :: &
92 : firstice ! if true, initialized values should be used
93 :
94 : ! local variables
95 :
96 : real (kind=dbl_kind) :: &
97 56339 : hin_old ! ice thickness before current melt/growth (m)
98 :
99 : character(len=*),parameter :: subname='(preflushing_changes)'
100 :
101 : !-----------------------------------------------------------------
102 : ! initialize
103 : !-----------------------------------------------------------------
104 :
105 207881 : if (fbri <= c0) then
106 0 : write(warnstr, *) subname,'fbri, hice_old', fbri, hice_old
107 0 : call icepack_warnings_add(warnstr)
108 0 : write(warnstr, *) subname,'vicen, aicen', vicen, aicen
109 0 : call icepack_warnings_add(warnstr)
110 0 : call icepack_warnings_add(subname//' icepack_brine preflushing: fbri <= c0')
111 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
112 : endif
113 :
114 207881 : hin = vicen / aicen
115 207881 : hsn = vsnon / aicen
116 207881 : hin_old = max(c0, hin + meltb + meltt - congel - snoice)
117 207881 : dhice = hin_old - hice_old ! change due to subl/cond
118 207881 : dhbr_top = meltt - snoice - dhice
119 207881 : dhbr_bot = congel - meltb
120 :
121 207881 : if ((hice_old < puny) .OR. (hin_old < puny) .OR. firstice) then
122 3769 : hice_old = hin
123 3769 : dhbr_top = c0
124 3769 : dhbr_bot = c0
125 3769 : dhice = c0
126 3769 : fbri = c1
127 : endif
128 :
129 207881 : hbr_old = fbri * hice_old
130 :
131 207881 : end subroutine preflushing_changes
132 :
133 : !=======================================================================
134 :
135 : ! Computes ice microstructural properties for updating hbrine
136 : !
137 : ! NOTE: This subroutine uses thermosaline_vertical output to compute
138 : ! average ice permeability and the surface ice porosity
139 :
140 207881 : subroutine compute_microS_mushy (nilyr, nblyr, &
141 415762 : bgrid, cgrid, igrid, &
142 207881 : trcrn, hice_old, hbr_old, &
143 207881 : sss, sst, bTin, &
144 207881 : iTin, bphin, &
145 : kperm, bphi_min, &
146 415762 : bSin, brine_sal, brine_rho, &
147 207881 : iphin, ibrine_rho, ibrine_sal, &
148 207881 : sice_rho, iDin )
149 :
150 : integer (kind=int_kind), intent(in) :: &
151 : nilyr , & ! number of ice layers
152 : nblyr ! number of bio layers
153 :
154 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
155 : bgrid ! biology nondimensional vertical grid points
156 :
157 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
158 : igrid ! biology vertical interface points
159 :
160 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
161 : cgrid ! CICE vertical coordinate
162 :
163 : real (kind=dbl_kind), &
164 : intent(in) :: &
165 : hice_old , & ! previous timestep ice height (m)
166 : sss , & ! ocean salinity (ppt)
167 : sst ! ocean temperature (C)
168 :
169 : real (kind=dbl_kind), dimension(ntrcr), &
170 : intent(in) :: &
171 : trcrn
172 :
173 : real (kind=dbl_kind), intent(out) :: &
174 : kperm , & ! average ice permeability (m^2)
175 : bphi_min ! surface porosity
176 :
177 : real (kind=dbl_kind), intent(inout) :: &
178 : hbr_old ! previous timestep brine height (m)
179 :
180 : real (kind=dbl_kind), dimension (nblyr+1), &
181 : intent(inout) :: &
182 : iDin ! tracer diffusivity/h^2 (1/s) includes gravity drainage/molecular
183 :
184 : real (kind=dbl_kind), dimension (nblyr+1), &
185 : intent(inout) :: &
186 : iphin , & ! porosity on the igrid
187 : ibrine_rho , & ! brine rho on interface
188 : ibrine_sal , & ! brine sal on interface
189 : iTin ! Temperature on the igrid (oC)
190 :
191 : real (kind=dbl_kind), dimension (nblyr+2), &
192 : intent(inout) :: &
193 : bSin , & ! bulk salinity (ppt) on bgrid
194 : brine_sal , & ! equilibrium brine salinity (ppt)
195 : brine_rho ! internal brine density (kg/m^3)
196 :
197 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
198 : bTin , & ! Temperature on bgrid
199 : bphin ! porosity on bgrid
200 :
201 : real (kind=dbl_kind), intent(inout) :: &
202 : sice_rho ! average ice density
203 :
204 : real (kind=dbl_kind), dimension (nblyr+2) :: &
205 929586 : bqin ! enthalpy on the bgrid ()
206 :
207 : real (kind=dbl_kind), dimension (nblyr+1) :: &
208 760569 : ikin ! permeability (m^2)
209 :
210 : integer (kind=int_kind) :: &
211 : k ! vertical biology layer index
212 :
213 : real (kind=dbl_kind) :: &
214 56339 : surface_S , & ! salinity of ice above hin > hbr
215 56339 : hinc_old ! mean ice thickness before current melt/growth (m)
216 :
217 : real (kind=dbl_kind), dimension (ntrcr+2) :: & ! nblyr+2)
218 12057411 : trtmp_s , & ! temporary, remapped tracers
219 12057411 : trtmp_q ! temporary, remapped tracers
220 :
221 : real (kind=dbl_kind), dimension(nblyr+1) :: &
222 721705 : drho ! brine density difference (kg/m^3)
223 :
224 : real(kind=dbl_kind), parameter :: &
225 : Smin = p01
226 :
227 : character(len=*),parameter :: subname='(compute_microS_mushy)'
228 :
229 : !-----------------------------------------------------------------
230 : ! Define ice salinity and temperature on bgrid
231 : !-----------------------------------------------------------------
232 :
233 48882117 : trtmp_s(:) = c0
234 48882117 : trtmp_q(:) = c0
235 1821363 : iDin(:) = c0
236 :
237 : ! map Sin and qin (cice) profiles to bgc grid
238 207881 : surface_S = min_salin
239 207881 : hbr_old = min(hbr_old, maxhbr*hice_old)
240 207881 : hinc_old = hice_old
241 :
242 : call remap_zbgc(nilyr, &
243 : nt_sice, &
244 0 : trcrn, trtmp_s, &
245 : 0, nblyr, &
246 : hinc_old, hinc_old, &
247 0 : cgrid(2:nilyr+1), &
248 207881 : bgrid(2:nblyr+1), surface_S )
249 207881 : if (icepack_warnings_aborted(subname)) return
250 :
251 : call remap_zbgc(nilyr, &
252 : nt_qice, &
253 0 : trcrn, trtmp_q, &
254 : 0, nblyr, &
255 : hinc_old, hinc_old, &
256 0 : cgrid(2:nilyr+1), &
257 207881 : bgrid(2:nblyr+1), surface_S )
258 207881 : if (icepack_warnings_aborted(subname)) return
259 :
260 1613482 : do k = 1, nblyr
261 1405601 : bqin (k+1) = min(c0, trtmp_q(nt_qice+k-1))
262 1405601 : bSin (k+1) = max(Smin, trtmp_s(nt_sice+k-1))
263 1405601 : bTin (k+1) = icepack_mushy_temperature_mush(bqin(k+1), bSin(k+1))
264 1613482 : bphin(k+1) = icepack_mushy_liquid_fraction (bTin(k+1), bSin(k+1))
265 : enddo ! k
266 :
267 207881 : bSin (1) = bSin(2)
268 207881 : bTin (1) = bTin(2)
269 207881 : bphin(1) = bphin(2)
270 207881 : bphin(nblyr+2) = c1
271 207881 : bSin (nblyr+2) = sss
272 207881 : bTin (nblyr+2) = sst
273 207881 : bphin(nblyr+2) = c1
274 :
275 : !-----------------------------------------------------------------
276 : ! Define ice multiphase structure
277 : !-----------------------------------------------------------------
278 :
279 : call prepare_hbrine (nblyr, &
280 0 : bSin, bTin, iTin, &
281 0 : brine_sal, brine_rho, &
282 0 : ibrine_sal, ibrine_rho, &
283 : sice_rho, &
284 0 : bphin, iphin, &
285 : kperm, bphi_min, &
286 207881 : igrid, sss)
287 207881 : if (icepack_warnings_aborted(subname)) return
288 :
289 : call calculate_drho(nblyr, igrid, bgrid, &
290 207881 : brine_rho, ibrine_rho, drho)
291 207881 : if (icepack_warnings_aborted(subname)) return
292 :
293 1613482 : do k= 2, nblyr+1
294 1405601 : ikin(k) = k_o*iphin(k)**exp_h
295 1405601 : iDin(k) = iphin(k)*Dm/hbr_old**2
296 1405601 : if (hbr_old .GE. Ra_c) &
297 344807 : iDin(k) = iDin(k) &
298 1613482 : + l_sk*ikin(k)*gravit/viscos_dynamic*drho(k)/hbr_old**2
299 : enddo ! k
300 :
301 : end subroutine compute_microS_mushy
302 :
303 : !=======================================================================
304 :
305 207881 : subroutine prepare_hbrine (nblyr, &
306 623643 : bSin, bTin, iTin, &
307 415762 : brine_sal, brine_rho, &
308 207881 : ibrine_sal, ibrine_rho, &
309 415762 : sice_rho, bphin, iphin,&
310 : kperm, bphi_min, &
311 207881 : i_grid, sss)
312 :
313 : integer (kind=int_kind), intent(in) :: &
314 : nblyr ! number of bio layers
315 :
316 : real (kind=dbl_kind), dimension (:), &
317 : intent(in) :: &
318 : bSin , & ! salinity of ice layers on bio grid (ppt)
319 : bTin , & ! temperature of ice layers on bio grid for history (C)
320 : i_grid ! biology grid interface points
321 :
322 : real (kind=dbl_kind), dimension (:), &
323 : intent(inout) :: &
324 : brine_sal , & ! equilibrium brine salinity (ppt)
325 : brine_rho , & ! internal brine density (kg/m^3)
326 : ibrine_rho , & ! brine density on interface (kg/m^3)
327 : ibrine_sal , & ! brine salinity on interface (ppt)
328 : iphin , & ! porosity on interface
329 : iTin , & ! Temperature on interface
330 : bphin ! porosity of layers
331 :
332 : real (kind=dbl_kind), intent(in) :: &
333 : sss ! sea surface salinity (ppt)
334 :
335 : real (kind=dbl_kind), intent(out) :: &
336 : kperm , & ! harmonic average permeability (m^2)
337 : bphi_min ! minimum porosity
338 :
339 : real (kind=dbl_kind), intent(inout) :: &
340 : sice_rho ! avg sea ice density
341 :
342 : ! local variables
343 :
344 : real (kind=dbl_kind), dimension(nblyr+1) :: &
345 760569 : kin ! permeability
346 :
347 : real (kind=dbl_kind) :: &
348 112678 : k_min, ktemp, &
349 56339 : igrp, igrm, rigr ! grid finite differences
350 :
351 : integer (kind=int_kind) :: &
352 : k ! layer index
353 :
354 : character(len=*),parameter :: subname='(prepare_hbrine)'
355 :
356 : !-----------------------------------------------------------------
357 : ! calculate equilibrium brine density and gradients
358 : !-----------------------------------------------------------------
359 :
360 207881 : sice_rho = c0
361 :
362 1821363 : do k = 1, nblyr+1
363 :
364 1613482 : if (k == 1) then
365 207881 : igrm = 0
366 : else
367 1405601 : igrm = i_grid(k) - i_grid(k-1)
368 : endif
369 :
370 401146 : brine_sal(k) = a1*bTin(k) &
371 401146 : + a2*bTin(k)**2 &
372 1613482 : + a3*bTin(k)**3
373 1613482 : brine_rho(k) = b1 + b2*brine_sal(k)
374 401146 : bphin (k) = max(puny, bSin(k)*rhosi &
375 1613482 : / (brine_sal(k)*brine_rho(k)))
376 1613482 : bphin (k) = min(c1, bphin(k))
377 1613482 : kin (k) = k_o*bphin(k)**exp_h
378 401146 : sice_rho = sice_rho + (rhoi*(c1-bphin(k)) &
379 1821363 : + brine_rho(k)*bphin(k))*igrm
380 : enddo ! k
381 :
382 207881 : brine_sal (nblyr+2) = sss
383 207881 : brine_rho (nblyr+2) = rhow
384 207881 : bphin (nblyr+2) = c1
385 207881 : ibrine_sal(1) = brine_sal (2)
386 207881 : ibrine_sal(nblyr+1) = brine_sal (nblyr+2)
387 207881 : ibrine_rho(1) = brine_rho (2)
388 207881 : ibrine_rho(nblyr+1) = brine_rho (nblyr+2)
389 207881 : iTin (1) = bTin(2)
390 207881 : iTin (nblyr+1) = bTin(nblyr+1)
391 207881 : iphin (1) = bphin (2)
392 207881 : iphin (nblyr+1) = bphin (nblyr+1)
393 1821363 : k_min = MINVAL(kin(2:nblyr+1))
394 207881 : kperm = c0 ! initialize
395 207881 : ktemp = c0
396 207881 : bphi_min = bphin (1)
397 : ! bphi_min = max(bphin(1),bSin(2)*rhosi/bphin(2) &
398 : ! / (brine_sal(1)*brine_rho(1))*phi_snow)
399 :
400 1405601 : do k = 2, nblyr
401 1197720 : if (k_min > c0) then
402 1197720 : ktemp = ktemp + c1/kin(k)
403 1197720 : kperm = k_min
404 : endif
405 :
406 1197720 : igrp = i_grid(k+1) - i_grid(k )
407 1197720 : igrm = i_grid(k ) - i_grid(k-1)
408 1197720 : rigr = c1 / (i_grid(k+1)-i_grid(k-1))
409 :
410 1197720 : ibrine_sal(k) = (brine_sal(k+1)*igrp + brine_sal(k)*igrm) * rigr
411 1197720 : ibrine_rho(k) = (brine_rho(k+1)*igrp + brine_rho(k)*igrm) * rigr
412 1197720 : iTin (k) = (bTin (k+1)*igrp + bTin (k)*igrm) * rigr
413 288468 : iphin (k) = max(puny, &
414 1197720 : (bphin (k+1)*igrp + bphin (k)*igrm) * rigr)
415 1405601 : iphin (k) = min(c1, iphin (k))
416 : enddo ! k
417 :
418 207881 : if (k_min > c0) then
419 207881 : ktemp = ktemp + c1/kin(nblyr+1)
420 207881 : kperm = real(nblyr,kind=dbl_kind)/ktemp
421 : endif
422 :
423 207881 : end subroutine prepare_hbrine
424 :
425 : !=======================================================================
426 :
427 : ! Changes include brine height increases from ice and snow surface melt,
428 : ! congelation growth, and upward pressure driven flow from snow loading.
429 : !
430 : ! Decreases arise from downward flushing and bottom melt.
431 : !
432 : ! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice
433 : ! with dynamic salinity or the height ratio == hbr/vicen*aicen, where
434 : ! hbr is the height of the brine surface relative to the bottom of the
435 : ! ice. This volume fraction may be > 1 in which case there is brine
436 : ! above the ice surface (ponds).
437 :
438 207881 : subroutine update_hbrine (meltt, &
439 : melts, dt, &
440 : hin, hsn, &
441 : hin_old, hbr, &
442 : hbr_old, &
443 : fbri, &
444 : dhS_top, dhS_bottom, &
445 : dh_top_chl, dh_bot_chl, &
446 : kperm, bphi_min, &
447 : darcy_V, darcy_V_chl, &
448 : bphin, aice0, &
449 : dh_direct)
450 :
451 : real (kind=dbl_kind), intent(in) :: &
452 : dt ! timestep
453 :
454 : real (kind=dbl_kind), intent(in):: &
455 : meltt, & ! true top melt over dt (m)
456 : melts, & ! true snow melt over dt (m)
457 : hin, & ! ice thickness (m)
458 : hsn, & ! snow thickness (m)
459 : hin_old, & ! past timestep ice thickness (m)
460 : hbr_old, & ! previous timestep hbr
461 : kperm, & ! avg ice permeability
462 : bphin, & ! upper brine porosity
463 : dhS_bottom, & ! change in bottom hbr initially before darcy flow
464 : aice0 ! open water area fraction
465 :
466 : real (kind=dbl_kind), intent(inout):: &
467 : darcy_V , & ! Darcy velocity: m/s
468 : darcy_V_chl, & ! Darcy velocity: m/s for bgc
469 : dhS_top , & ! change in top hbr before darcy flow
470 : dh_bot_chl , & ! change in bottom for algae
471 : dh_top_chl , & ! change in bottom for algae
472 : hbr , & ! thickness of brine (m)
473 : fbri , & ! brine height ratio tracer (hbr/hin)
474 : bphi_min ! surface porosity
475 :
476 : real (kind=dbl_kind), intent(out):: &
477 : dh_direct ! surface flooding or runoff (m)
478 :
479 : ! local variables
480 :
481 : real (kind=dbl_kind) :: &
482 56339 : hbrmin , & ! thinS or hin
483 56339 : dhbr_hin , & ! hbr-hin
484 56339 : hbrocn , & ! brine height above sea level (m) hbr-h_ocn
485 56339 : dhbr , & ! change in brine surface
486 56339 : h_ocn , & ! new ocean surface from ice bottom (m)
487 56339 : darcy_coeff, & ! magnitude of the Darcy velocity/hbrocn (1/s)
488 56339 : hbrocn_new , & ! hbrocn after flushing
489 56339 : dhflood , & ! surface flooding by ocean
490 56339 : exp_arg , & ! temporary exp value
491 56339 : dhrunoff ! direct runoff to ocean
492 :
493 : real (kind=dbl_kind), parameter :: &
494 : dh_min = p001 ! brine remains within dh_min of sea level
495 : ! when ice thickness is less than thinS
496 :
497 : character(len=*),parameter :: subname='(update_hbrine)'
498 :
499 207881 : hbrocn = c0
500 207881 : darcy_V = c0
501 207881 : darcy_V_chl = c0
502 207881 : hbrocn_new = c0
503 207881 : h_ocn = rhosi/rhow*hin + rhos/rhow*hsn
504 207881 : dh_direct = c0
505 :
506 207881 : if (hbr_old > thinS .AND. hin_old > thinS .AND. hin > thinS ) then
507 207881 : hbrmin = thinS
508 207881 : dhS_top = -max(c0, min(hin_old-hbr_old, meltt)) * rhoi/rhow
509 207881 : dhS_top = dhS_top - max(c0, melts) * rhos/rhow
510 207881 : dh_top_chl = dhS_top
511 207881 : dhbr = dhS_bottom - dhS_top
512 207881 : hbr = max(puny, hbr_old+dhbr)
513 207881 : hbrocn = hbr - h_ocn
514 207881 : darcy_coeff = max(c0, kperm*gravit/(viscos*hbr_old))
515 :
516 207881 : if (hbrocn > c0 .AND. hbr > thinS ) then
517 61751 : bphi_min = bphin
518 61751 : dhrunoff = -dhS_top*aice0
519 61751 : hbrocn = max(c0,hbrocn - dhrunoff)
520 61751 : exp_arg = darcy_coeff/bphi_min*dt
521 : ! tcraig avoids underflows
522 61751 : if (exp_arg > exp_argmax) then
523 1831 : hbrocn_new = c0
524 : else
525 59920 : hbrocn_new = hbrocn*exp(-exp_arg)
526 : endif
527 61751 : hbr = max(hbrmin, h_ocn + hbrocn_new)
528 61751 : hbrocn_new = hbr-h_ocn
529 61751 : darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn)
530 61751 : darcy_V_chl= darcy_V
531 61751 : dhS_top = dhS_top - darcy_V*dt/bphi_min + dhrunoff
532 61751 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min + dhrunoff
533 61751 : dh_direct = dhrunoff
534 146130 : elseif (hbrocn < c0 .AND. hbr > thinS) then
535 146130 : exp_arg = darcy_coeff/bphi_min*dt
536 : ! tcraig avoids underflows
537 146130 : if (exp_arg > exp_argmax) then
538 1612 : hbrocn_new = c0
539 : else
540 144518 : hbrocn_new = hbrocn*exp(-exp_arg)
541 : endif
542 146130 : dhflood = max(c0,hbrocn_new - hbrocn)*aice0
543 146130 : hbr = max(hbrmin, h_ocn + hbrocn_new)
544 146130 : darcy_V = -SIGN((hbrocn-hbrocn_new + dhflood)/dt*bphi_min, hbrocn)
545 146130 : darcy_V_chl= darcy_V
546 146130 : dhS_top = dhS_top - darcy_V*dt/bphi_min - dhflood
547 146130 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min - dhflood
548 146130 : dh_direct = -dhflood
549 : endif
550 :
551 207881 : dh_bot_chl = dhS_bottom
552 :
553 : else ! very thin brine height
554 0 : hbrmin = min(thinS, hin)
555 0 : hbr = max(hbrmin, hbr_old+dhS_bottom-dhS_top)
556 0 : dhbr_hin = hbr - h_ocn
557 0 : if (abs(dhbr_hin) > dh_min) &
558 0 : hbr = max(hbrmin, h_ocn + SIGN(dh_min,dhbr_hin))
559 0 : dhS_top = hbr_old - hbr + dhS_bottom
560 0 : dh_top_chl = dhS_top
561 0 : dh_bot_chl = dhS_bottom
562 : endif
563 :
564 207881 : fbri = hbr/hin
565 :
566 207881 : end subroutine update_hbrine
567 :
568 : !=======================================================================
569 : !
570 : ! Computes ice microstructural properties for zbgc and zsalinity
571 : !
572 : ! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice with
573 : ! dynamic salinity or the height ratio == hbr/vicen*aicen, where hbr is the
574 : ! height of the brine surface relative to the bottom of the ice.
575 : ! This volume fraction
576 : ! may be > 1 in which case there is brine above the ice surface (meltponds).
577 : !
578 0 : subroutine compute_microS (n_cat, nilyr, nblyr, &
579 0 : bgrid, cgrid, igrid, &
580 0 : trcrn, hice_old, &
581 : hbr_old, sss, sst, &
582 0 : bTin, iTin, bphin, &
583 : kperm, bphi_min, &
584 : Rayleigh_criteria, firstice, &
585 0 : bSin, brine_sal, &
586 0 : brine_rho, iphin, ibrine_rho, &
587 0 : ibrine_sal, sice_rho, sloss)
588 :
589 : integer (kind=int_kind), intent(in) :: &
590 : n_cat , & ! ice category
591 : nilyr , & ! number of ice layers
592 : nblyr ! number of bio layers
593 :
594 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
595 : bgrid ! biology nondimensional vertical grid points
596 :
597 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
598 : igrid ! biology vertical interface points
599 :
600 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
601 : cgrid ! CICE vertical coordinate
602 :
603 : real (kind=dbl_kind), intent(in) :: &
604 : hice_old , & ! previous timestep ice height (m)
605 : sss , & ! ocean salinity (ppt)
606 : sst ! ocean temperature (oC)
607 :
608 : real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: &
609 : trcrn
610 :
611 : real (kind=dbl_kind), intent(inout) :: &
612 : hbr_old , & ! old brine height
613 : sice_rho ! average ice density
614 :
615 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
616 : bTin , & ! Temperature of ice layers on bio grid for history file (^oC)
617 : bphin , & ! Porosity of layers
618 : brine_sal , & ! equilibrium brine salinity (ppt)
619 : brine_rho ! Internal brine density (kg/m^3)
620 :
621 : real (kind=dbl_kind), dimension (nblyr+2), intent(out) :: &
622 : bSin
623 :
624 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
625 : iTin ! Temperature on the interface grid
626 :
627 : real (kind=dbl_kind), intent(out) :: &
628 : bphi_min , & ! surface porosity
629 : kperm , & ! average ice permeability (m^2)
630 : sloss ! (g/m^2) salt from brine runoff for hbr > maxhbr*hin
631 :
632 : logical (kind=log_kind), intent(inout) :: &
633 : Rayleigh_criteria ! .true. if ice exceeded a minimum thickness hin >= Ra_c
634 :
635 : logical (kind=log_kind), intent(in) :: &
636 : firstice ! .true. if ice is newly formed
637 :
638 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
639 : iphin , & ! porosity on the igrid
640 : ibrine_rho , & ! brine rho on interface
641 : ibrine_sal ! brine sal on interface
642 :
643 : ! local variables
644 :
645 : integer (kind=int_kind) :: &
646 : k ! vertical biology layer index
647 :
648 : real (kind=dbl_kind) :: &
649 0 : surface_S , & ! salinity of ice above hin > hbr
650 0 : hinc_old ! ice thickness (cell quantity) before current melt/growth (m)
651 :
652 : ! logical (kind=log_kind) :: &
653 : ! Rayleigh ! .true. if ice exceeded a minimum thickness hin >= Ra_c
654 :
655 : real (kind=dbl_kind), dimension (ntrcr+2) :: &
656 0 : trtmp0 , & ! temporary, remapped tracers
657 0 : trtmp ! temporary, remapped tracers
658 :
659 : real (kind=dbl_kind) :: &
660 0 : Tmlts ! melting temperature
661 :
662 : character(len=*),parameter :: subname='(compute_microS)'
663 :
664 : !-----------------------------------------------------------------
665 : ! Initialize
666 : !-----------------------------------------------------------------
667 :
668 0 : sloss = c0
669 0 : bTin(:) = c0
670 0 : bSin(:) = c0
671 :
672 0 : trtmp(:) = c0
673 0 : surface_S = min_salin
674 :
675 0 : hinc_old = hice_old
676 :
677 : !-----------------------------------------------------------------
678 : ! Rayleigh condition for salinity and bgc:
679 : ! Implemented as a minimum thickness criteria for category 1 ice only.
680 : ! When hin >= Ra_c (m), pressure flow is allowed.
681 : ! Turn off by putting Ra_c = 0 in ice_in namelist.
682 : !-----------------------------------------------------------------
683 :
684 : ! Rayleigh = .true.
685 : ! if (n_cat == 1 .AND. hbr_old < Ra_c) then
686 : ! Rayleigh = Rayleigh_criteria ! only category 1 ice can be false
687 : ! endif
688 :
689 : !-----------------------------------------------------------------
690 : ! Define ice salinity on Sin
691 : !-----------------------------------------------------------------
692 :
693 0 : if (firstice) then
694 :
695 0 : do k = 1, nilyr
696 0 : trcrn(nt_sice+k-1) = sss*salt_loss
697 : enddo
698 :
699 : call remap_zbgc(nilyr, &
700 : nt_sice, &
701 0 : trcrn, trtmp, &
702 : 0, nblyr, &
703 : hinc_old, hinc_old, &
704 0 : cgrid(2:nilyr+1), &
705 0 : bgrid(2:nblyr+1), surface_S )
706 0 : if (icepack_warnings_aborted(subname)) return
707 :
708 0 : do k = 1, nblyr
709 0 : trcrn(nt_bgc_S+k-1) = max(min_salin,trtmp(nt_sice+k-1))
710 0 : bSin(k+1) = max(min_salin,trcrn(nt_bgc_S+k-1))
711 0 : if (trcrn(nt_bgc_S+k-1) < min_salin-puny) &
712 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
713 : enddo ! k
714 :
715 0 : bSin(1) = bSin(2)
716 0 : bSin(nblyr+2) = sss
717 :
718 0 : elseif (hbr_old > maxhbr*hice_old) then
719 :
720 : call remap_zbgc(nblyr, &
721 : nt_bgc_S, &
722 0 : trcrn, trtmp, &
723 : 0, nblyr, &
724 : hbr_old, &
725 : maxhbr*hinc_old, &
726 0 : bgrid(2:nblyr+1), &
727 0 : bgrid(2:nblyr+1), surface_S )
728 0 : if (icepack_warnings_aborted(subname)) return
729 :
730 0 : do k = 1, nblyr
731 0 : bSin(k+1) = max(min_salin,trtmp(nt_bgc_S+k-1))
732 0 : sloss = sloss + rhosi*(hbr_old*trcrn(nt_bgc_S+k-1) &
733 0 : - maxhbr*hice_old*bSin(k+1))*(igrid(k+1)-igrid(k))
734 0 : trcrn(nt_bgc_S+k-1) = bSin(k+1)
735 0 : if (trcrn(nt_bgc_S+k-1) < min_salin-puny) &
736 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
737 : enddo ! k
738 :
739 0 : bSin(1) = bSin(2)
740 0 : bSin(nblyr+2) = sss
741 0 : hbr_old = maxhbr*hinc_old
742 :
743 : else ! old, thin ice
744 :
745 0 : do k = 1, nblyr
746 0 : trcrn(nt_bgc_S+k-1) = max(min_salin,trcrn(nt_bgc_S+k-1))
747 0 : bSin (k+1) = trcrn(nt_bgc_S+k-1)
748 : enddo ! k
749 :
750 0 : bSin (1) = bSin(2)
751 0 : bSin (nblyr+2) = sss
752 :
753 : endif ! ice type
754 :
755 : !-----------------------------------------------------------------
756 : ! sea ice temperature for bio grid
757 : !-----------------------------------------------------------------
758 :
759 0 : do k = 1, nilyr
760 0 : Tmlts = -trcrn(nt_sice+k-1)*depressT
761 0 : trtmp0(nt_qice+k-1) = calculate_Tin_from_qin(trcrn(nt_qice+k-1),Tmlts)
762 : enddo ! k
763 :
764 0 : trtmp(:) = c0
765 :
766 : ! CICE to Bio: remap temperatures
767 : call remap_zbgc (nilyr, nt_qice, &
768 0 : trtmp0(1:ntrcr), trtmp, &
769 : 0, nblyr, &
770 : hinc_old, hbr_old, &
771 0 : cgrid(2:nilyr+1), &
772 0 : bgrid(2:nblyr+1), surface_S )
773 0 : if (icepack_warnings_aborted(subname)) return
774 :
775 0 : do k = 1, nblyr
776 0 : Tmlts = -bSin(k+1) * depressT
777 0 : bTin (k+1) = min(Tmlts,trtmp(nt_qice+k-1))
778 : enddo !k
779 :
780 0 : Tmlts = -min_salin* depressT
781 0 : bTin (1) = min(Tmlts,(bTin(2) + trcrn(nt_Tsfc))*p5)
782 0 : Tmlts = -bSin(nblyr+2)* depressT
783 0 : bTin (nblyr+2) = sst
784 :
785 : !-----------------------------------------------------------------
786 : ! Define ice multiphase structure
787 : !-----------------------------------------------------------------
788 :
789 : call prepare_hbrine (nblyr, &
790 0 : bSin, bTin, iTin, &
791 0 : brine_sal, brine_rho, &
792 0 : ibrine_sal, ibrine_rho, &
793 : sice_rho, &
794 0 : bphin, iphin, &
795 : kperm, bphi_min, &
796 0 : igrid, sss)
797 0 : if (icepack_warnings_aborted(subname)) return
798 :
799 : end subroutine compute_microS
800 :
801 : !==========================================================================================
802 : !
803 : ! Find density difference about interface grid points
804 : ! for gravity drainage parameterization
805 :
806 623643 : subroutine calculate_drho (nblyr, i_grid, b_grid, &
807 415762 : brine_rho, ibrine_rho, drho)
808 :
809 : integer (kind=int_kind), intent(in) :: &
810 : nblyr ! number of bio layers
811 :
812 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
813 : b_grid ! biology nondimensional grid layer points
814 :
815 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
816 : i_grid ! biology grid interface points
817 :
818 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
819 : brine_rho ! Internal brine density (kg/m^3)
820 :
821 : real (kind=dbl_kind), dimension (nblyr + 1), intent(in) :: &
822 : ibrine_rho ! Internal brine density (kg/m^3)
823 :
824 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
825 : drho ! brine difference about grid point (kg/m^3)
826 :
827 : ! local variables
828 :
829 : integer (kind=int_kind) :: &
830 : k, mm ! indices
831 :
832 : integer (kind=int_kind) :: &
833 : mstop, mstart
834 :
835 : real (kind=dbl_kind), dimension (nblyr + 1) :: & !on the zbgc vertical grid
836 760569 : rho_a , & ! average brine density above grid point (kg/m^3)
837 760569 : rho_2a ! average brine density above and below grid points (kg/m^3)
838 :
839 : real (kind=dbl_kind), dimension (nblyr + 1) :: & !on the zbgc vertical grid
840 760569 : rho_b , & ! brine density above grid point (kg/m^3)
841 816908 : rho_2b ! brine density above and below grid points (kg/m^3)
842 :
843 : character(len=*),parameter :: subname='(calculate_drho)'
844 :
845 1821363 : rho_a (:) = c0
846 1821363 : rho_2a(:) = c0
847 1821363 : rho_b (:) = c0
848 1821363 : rho_2b(:) = c0
849 1821363 : drho (:) = c0 ! surface is snow or atmosphere
850 :
851 1821363 : do k = 1, nblyr+1 ! i_grid values
852 :
853 : !----------------------------------------------
854 : ! h_avg(k) = i_grid(k)
855 : ! Calculate rho_a(k), ie average rho above i_grid(k)
856 : ! first part is good
857 : !----------------------------------------------
858 :
859 1821363 : if (k == 2) then
860 56339 : rho_a(2) = (brine_rho(2)*b_grid(2) &
861 56339 : + (ibrine_rho(2) + brine_rho(2)) &
862 207881 : * p5*(i_grid(2)-b_grid(2)) )/i_grid(2)
863 207881 : rho_b(2) = brine_rho(2)
864 :
865 1405601 : elseif (k > 2 .AND. k < nblyr+1) then
866 961560 : rho_a(k) = (rho_a(k-1)*i_grid(k-1) + (ibrine_rho(k-1) + brine_rho(k)) &
867 480780 : * p5*(b_grid(k)-i_grid(k-1)) + (ibrine_rho(k ) + brine_rho(k)) &
868 1959660 : * p5*(i_grid(k)-b_grid(k)))/i_grid(k)
869 998100 : rho_b(k) = brine_rho(k)
870 : else
871 417668 : rho_a(nblyr+1) = (rho_a(nblyr)*i_grid(nblyr) + (ibrine_rho(nblyr) + &
872 313251 : brine_rho(nblyr+1))*p5*(b_grid(nblyr+1)-i_grid(nblyr)) + &
873 1138420 : brine_rho(nblyr+1)*(i_grid(nblyr+1)-b_grid(nblyr+1)))/i_grid(nblyr+1)
874 407501 : rho_a(1) = brine_rho(2) !for k == 1 use grid point value
875 407501 : rho_b(nblyr+1) = brine_rho(nblyr+1)
876 407501 : rho_b(1) = brine_rho(2)
877 : endif
878 :
879 : enddo !k
880 :
881 : !----------------------------------------------
882 : ! Calculate average above and below k rho_2a
883 : !----------------------------------------------
884 :
885 1821363 : do k = 1, nblyr+1 !i_grid values
886 1613482 : if (k == 1) then
887 56339 : rho_2a(1) = (rho_a(1)*b_grid(2) + p5*(brine_rho(2) + ibrine_rho(2)) &
888 207881 : * (i_grid(2)-b_grid(2)))/i_grid(2)
889 207881 : rho_2b(1) = brine_rho(2)
890 : else
891 1405601 : mstop = 2*(k-1) + 1
892 1405601 : if (mstop < nblyr+1) then
893 598860 : rho_2a(k) = rho_a(mstop)
894 598860 : mstart = 2
895 598860 : mstop = 1
896 : else
897 806741 : mstart = nblyr+2
898 806741 : mstop = nblyr+3
899 : endif
900 :
901 3019083 : do mm = mstart,mstop
902 3019083 : rho_2a(k) =(rho_a(nblyr+1) + rhow*(c2*i_grid(k)-c1))*p5/i_grid(k)
903 : enddo
904 1405601 : rho_2b(k) = brine_rho(k+1)
905 : endif
906 401146 : drho(k) = max(rho_b(k) - rho_2b(k),max(c0,c2*(rho_a(k)-rho_2a(k)), &
907 1821363 : c2*(brine_rho(k)-brine_rho(k+1))/real(nblyr,kind=dbl_kind)))
908 : enddo
909 :
910 207881 : end subroutine calculate_drho
911 :
912 : !=======================================================================
913 : !autodocument_start icepack_init_hbrine
914 : ! Initialize brine height tracer
915 :
916 5 : subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
917 5 : icgrid, swgrid, nblyr, nilyr, phi_snow)
918 :
919 : integer (kind=int_kind), intent(in) :: &
920 : nilyr, & ! number of ice layers
921 : nblyr ! number of bio layers
922 :
923 : real (kind=dbl_kind), intent(inout) :: &
924 : phi_snow !porosity at the ice-snow interface
925 :
926 : real (kind=dbl_kind), dimension (nblyr+2), intent(out) :: &
927 : bgrid ! biology nondimensional vertical grid points
928 :
929 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
930 : igrid ! biology vertical interface points
931 :
932 : real (kind=dbl_kind), dimension (nilyr+1), intent(out) :: &
933 : cgrid , & ! CICE vertical coordinate
934 : icgrid , & ! interface grid for CICE (shortwave variable)
935 : swgrid ! grid for ice tracers used in dEdd scheme
936 :
937 : !autodocument_end
938 :
939 : ! local variables
940 :
941 : integer (kind=int_kind) :: &
942 : k ! vertical index
943 :
944 : real (kind=dbl_kind) :: &
945 3 : zspace ! grid spacing for CICE vertical grid
946 :
947 : character(len=*),parameter :: subname='(icepack_init_hbrine)'
948 :
949 :
950 5 : if (phi_snow .le. c0) phi_snow = c1-rhos/rhoi
951 :
952 : !-----------------------------------------------------------------
953 : ! Calculate bio gridn: 0 to 1 corresponds to ice top to bottom
954 : !-----------------------------------------------------------------
955 :
956 44 : bgrid(:) = c0 ! zsalinity grid points
957 5 : bgrid(nblyr+2) = c1 ! bottom value
958 39 : igrid(:) = c0 ! bgc interface grid points
959 5 : igrid(1) = c0 ! ice top
960 5 : igrid(nblyr+1) = c1 ! ice bottom
961 :
962 5 : zspace = c1/max(c1,(real(nblyr,kind=dbl_kind)))
963 34 : do k = 2, nblyr+1
964 34 : bgrid(k) = zspace*(real(k,kind=dbl_kind) - c1p5)
965 : enddo
966 :
967 29 : do k = 2, nblyr
968 29 : igrid(k) = p5*(bgrid(k+1)+bgrid(k))
969 : enddo
970 :
971 : !-----------------------------------------------------------------
972 : ! Calculate CICE cgrid for interpolation ice top (0) to bottom (1)
973 : !-----------------------------------------------------------------
974 :
975 5 : cgrid(1) = c0 ! CICE vertical grid top point
976 5 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
977 :
978 40 : do k = 2, nilyr+1
979 40 : cgrid(k) = zspace * (real(k,kind=dbl_kind) - c1p5)
980 : enddo
981 :
982 : !-----------------------------------------------------------------
983 : ! Calculate CICE icgrid for ishortwave interpolation top(0) , bottom (1)
984 : !-----------------------------------------------------------------
985 :
986 5 : icgrid(1) = c0
987 5 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
988 :
989 40 : do k = 2, nilyr+1
990 40 : icgrid(k) = zspace * (real(k,kind=dbl_kind)-c1)
991 : enddo
992 :
993 : !------------------------------------------------------------------------
994 : ! Calculate CICE swgrid for dEdd ice: top of ice (0) , bottom of ice (1)
995 : ! Does not include snow
996 : ! see icepack_shortwave.F90
997 : ! swgrid represents the layer index of the delta-eddington ice layer index
998 : !------------------------------------------------------------------------
999 5 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
1000 5 : swgrid(1) = min(c1/60.0_dbl_kind, zspace/c2)
1001 5 : swgrid(2) = zspace/c2 !+ swgrid(1)
1002 35 : do k = 3, nilyr+1
1003 35 : swgrid(k) = zspace * (real(k,kind=dbl_kind)-c1p5)
1004 : enddo
1005 :
1006 5 : end subroutine icepack_init_hbrine
1007 :
1008 : !=======================================================================
1009 : !autodocument_start icepack_init_zsalinity
1010 : ! Initialize zSalinity
1011 :
1012 0 : subroutine icepack_init_zsalinity(nblyr,ntrcr_o, Rayleigh_criteria, &
1013 0 : Rayleigh_real, trcrn_bgc, nt_bgc_S, ncat, sss)
1014 :
1015 : integer (kind=int_kind), intent(in) :: &
1016 : nblyr, & ! number of biolayers
1017 : ntrcr_o, & ! number of non bio tracers
1018 : ncat , & ! number of categories
1019 : nt_bgc_S ! zsalinity index
1020 :
1021 : logical (kind=log_kind), intent(inout) :: &
1022 : Rayleigh_criteria
1023 :
1024 : real (kind=dbl_kind), intent(inout):: &
1025 : Rayleigh_real
1026 :
1027 : real (kind=dbl_kind), intent(in):: &
1028 : sss
1029 :
1030 : real (kind=dbl_kind), dimension(:,:), intent(inout):: &
1031 : trcrn_bgc ! bgc subset of trcrn
1032 :
1033 : !autodocument_end
1034 :
1035 : ! local variables
1036 :
1037 : integer (kind=int_kind) :: &
1038 : k, n
1039 :
1040 : character(len=*),parameter :: subname='(icepack_init_zsalinity)'
1041 :
1042 0 : if (nblyr .LE. 7) then
1043 0 : dts_b = 300.0_dbl_kind
1044 : else
1045 0 : dts_b = 50.0_dbl_kind
1046 : endif
1047 :
1048 0 : Rayleigh_criteria = .false. ! no ice initial condition
1049 0 : Rayleigh_real = c0
1050 0 : do n = 1,ncat
1051 0 : do k = 1,nblyr
1052 0 : trcrn_bgc(nt_bgc_S+k-1-ntrcr_o,n) = sss*salt_loss
1053 : enddo ! k
1054 : enddo ! n
1055 :
1056 0 : end subroutine icepack_init_zsalinity
1057 :
1058 : !=======================================================================
1059 :
1060 : end module icepack_brine
1061 :
1062 : !=======================================================================
|
