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
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, & ! LCOV_EXCL_LINE
29 : update_hbrine, & ! LCOV_EXCL_LINE
30 : calculate_drho, & ! LCOV_EXCL_LINE
31 : icepack_init_hbrine, & ! LCOV_EXCL_LINE
32 : icepack_init_zsalinity ! deprecated
33 :
34 : real (kind=dbl_kind), parameter :: &
35 : maxhbr = 1.25_dbl_kind , & ! brine overflows if hbr > maxhbr*hin ! LCOV_EXCL_LINE
36 : viscos = 2.1e-6_dbl_kind, & ! kinematic viscosity (m^2/s) ! LCOV_EXCL_LINE
37 : ! Brine salinity as a cubic function of temperature
38 : a1 = -21.4_dbl_kind , & ! (psu/C)
39 : a2 = -0.886_dbl_kind, & ! (psu/C^2) ! LCOV_EXCL_LINE
40 : a3 = -0.012_dbl_kind, & ! (psu/C^3) ! LCOV_EXCL_LINE
41 : ! Brine density as a quadratic of brine salinity
42 : b1 = 1000.0_dbl_kind, & ! (kg/m^3)
43 : b2 = 0.8_dbl_kind ! (kg/m^3/ppt)
44 :
45 : real (kind=dbl_kind), parameter :: &
46 : exp_argmax = 30.0_dbl_kind ! maximum argument of exponential for underflow
47 :
48 : !=======================================================================
49 :
50 : contains
51 :
52 : !=======================================================================
53 : ! Computes the top and bottom brine boundary changes for flushing
54 : ! works for zsalinity and tr_salinity
55 : !
56 : ! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice with
57 : ! dynamic salinity or the height ratio = hbr/vicen*aicen, where hbr is the
58 : ! height of the brine surface relative to the bottom of the ice. This volume fraction
59 : ! may be > 1 in which case there is brine above the ice surface (meltponds).
60 :
61 0 : subroutine preflushing_changes (aicen, vicen, vsnon, &
62 : meltb, meltt, congel, & ! LCOV_EXCL_LINE
63 : snoice, hice_old, dhice, & ! LCOV_EXCL_LINE
64 : fbri, dhbr_top, dhbr_bot, & ! LCOV_EXCL_LINE
65 : hbr_old, hin,hsn, firstice )
66 :
67 : real (kind=dbl_kind), intent(in) :: &
68 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
69 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
70 : vsnon , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
71 : meltb , & ! bottom ice melt (m) ! LCOV_EXCL_LINE
72 : meltt , & ! top ice melt (m) ! LCOV_EXCL_LINE
73 : congel , & ! bottom ice growth (m) ! LCOV_EXCL_LINE
74 : snoice ! top ice growth from flooding (m)
75 :
76 : real (kind=dbl_kind), intent(out) :: &
77 : hbr_old ! old brine height (m)
78 :
79 : real (kind=dbl_kind), intent(inout) :: &
80 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
81 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
82 : dhice ! change due to sublimation (<0)/condensation (>0) (m)
83 :
84 : real (kind=dbl_kind), intent(inout) :: &
85 : fbri , & ! trcrn(nt_fbri) ! LCOV_EXCL_LINE
86 : dhbr_top , & ! brine change in top for diagnostics (m) ! LCOV_EXCL_LINE
87 : dhbr_bot , & ! brine change in bottom for diagnostics (m) ! LCOV_EXCL_LINE
88 : hice_old ! old ice thickness (m)
89 :
90 : logical (kind=log_kind), intent(in) :: &
91 : firstice ! if true, initialized values should be used
92 :
93 : ! local variables
94 :
95 : real (kind=dbl_kind) :: &
96 0 : hin_old ! ice thickness before current melt/growth (m)
97 :
98 : character(len=*),parameter :: subname='(preflushing_changes)'
99 :
100 : !-----------------------------------------------------------------
101 : ! initialize
102 : !-----------------------------------------------------------------
103 :
104 0 : if (fbri <= c0) then
105 0 : write(warnstr, *) subname,'fbri, hice_old', fbri, hice_old
106 0 : call icepack_warnings_add(warnstr)
107 0 : write(warnstr, *) subname,'vicen, aicen', vicen, aicen
108 0 : call icepack_warnings_add(warnstr)
109 0 : call icepack_warnings_add(subname//' icepack_brine preflushing: fbri <= c0')
110 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
111 : endif
112 :
113 0 : hin = vicen / aicen
114 0 : hsn = vsnon / aicen
115 0 : hin_old = max(c0, hin + meltb + meltt - congel - snoice)
116 0 : dhice = hin_old - hice_old ! change due to subl/cond
117 0 : dhbr_top = meltt - snoice - dhice
118 0 : dhbr_bot = congel - meltb
119 :
120 0 : if ((hice_old < puny) .OR. (hin_old < puny) .OR. firstice) then
121 0 : hice_old = hin
122 0 : dhbr_top = c0
123 0 : dhbr_bot = c0
124 0 : dhice = c0
125 0 : fbri = c1
126 : endif
127 :
128 0 : hbr_old = fbri * hice_old
129 :
130 0 : end subroutine preflushing_changes
131 :
132 : !=======================================================================
133 : ! Computes ice microstructural properties for updating hbrine
134 : !
135 : ! NOTE: This subroutine uses thermosaline_vertical output to compute
136 : ! average ice permeability and the surface ice porosity
137 :
138 0 : subroutine compute_microS_mushy (nilyr, nblyr, &
139 : bgrid, cgrid, igrid, & ! LCOV_EXCL_LINE
140 : trcrn, hice_old, hbr_old, & ! LCOV_EXCL_LINE
141 : sss, sst, bTin, & ! LCOV_EXCL_LINE
142 : iTin, bphin, & ! LCOV_EXCL_LINE
143 : kperm, bphi_min, & ! LCOV_EXCL_LINE
144 : bSin, brine_sal, brine_rho, & ! LCOV_EXCL_LINE
145 : iphin, ibrine_rho, ibrine_sal, & ! LCOV_EXCL_LINE
146 0 : sice_rho, iDin )
147 :
148 : integer (kind=int_kind), intent(in) :: &
149 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
150 : nblyr ! number of bio layers
151 :
152 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
153 : bgrid ! biology nondimensional vertical grid points
154 :
155 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
156 : igrid ! biology vertical interface points
157 :
158 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
159 : cgrid ! CICE vertical coordinate
160 :
161 : real (kind=dbl_kind), intent(in) :: &
162 : hice_old , & ! previous timestep ice height (m) ! LCOV_EXCL_LINE
163 : sss , & ! ocean salinity (ppt) ! LCOV_EXCL_LINE
164 : sst ! ocean temperature (C)
165 :
166 : real (kind=dbl_kind), dimension(ntrcr), intent(in) :: &
167 : trcrn
168 :
169 : real (kind=dbl_kind), intent(out) :: &
170 : kperm , & ! average ice permeability (m^2) ! LCOV_EXCL_LINE
171 : bphi_min ! surface porosity
172 :
173 : real (kind=dbl_kind), intent(inout) :: &
174 : hbr_old ! previous timestep brine height (m)
175 :
176 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
177 : iDin ! tracer diffusivity/h^2 (1/s) includes gravity drainage/molecular
178 :
179 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
180 : iphin , & ! porosity on the igrid ! LCOV_EXCL_LINE
181 : ibrine_rho , & ! brine rho on interface ! LCOV_EXCL_LINE
182 : ibrine_sal , & ! brine sal on interface ! LCOV_EXCL_LINE
183 : iTin ! Temperature on the igrid (oC)
184 :
185 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
186 : bSin , & ! bulk salinity (ppt) on bgrid ! LCOV_EXCL_LINE
187 : brine_sal , & ! equilibrium brine salinity (ppt) ! LCOV_EXCL_LINE
188 : brine_rho ! internal brine density (kg/m^3)
189 :
190 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
191 : bTin , & ! Temperature on bgrid ! LCOV_EXCL_LINE
192 : bphin ! porosity on bgrid
193 :
194 : real (kind=dbl_kind), intent(inout) :: &
195 : sice_rho ! average ice density
196 :
197 : real (kind=dbl_kind), dimension (nblyr+2) :: &
198 0 : bqin ! enthalpy on the bgrid ()
199 :
200 : real (kind=dbl_kind), dimension (nblyr+1) :: &
201 0 : ikin ! permeability (m^2)
202 :
203 : integer (kind=int_kind) :: &
204 : k ! vertical biology layer index
205 :
206 : real (kind=dbl_kind) :: &
207 : surface_S , & ! salinity of ice above hin > hbr ! LCOV_EXCL_LINE
208 0 : hinc_old ! mean ice thickness before current melt/growth (m)
209 :
210 : real (kind=dbl_kind), dimension (ntrcr+2) :: & ! nblyr+2)
211 : trtmp_s , & ! temporary, remapped tracers ! LCOV_EXCL_LINE
212 0 : trtmp_q ! temporary, remapped tracers
213 :
214 : real (kind=dbl_kind), dimension(nblyr+1) :: &
215 0 : drho ! brine density difference (kg/m^3)
216 :
217 : real(kind=dbl_kind), parameter :: &
218 : Smin = p01
219 :
220 : character(len=*),parameter :: subname='(compute_microS_mushy)'
221 :
222 : !-----------------------------------------------------------------
223 : ! Define ice salinity and temperature on bgrid
224 : !-----------------------------------------------------------------
225 :
226 0 : trtmp_s(:) = c0
227 0 : trtmp_q(:) = c0
228 0 : iDin(:) = c0
229 :
230 : ! map Sin and qin (cice) profiles to bgc grid
231 0 : surface_S = min_salin
232 0 : hbr_old = min(hbr_old, maxhbr*hice_old)
233 0 : hinc_old = hice_old
234 :
235 : call remap_zbgc(nilyr, &
236 : nt_sice, & ! LCOV_EXCL_LINE
237 : trcrn, trtmp_s, & ! LCOV_EXCL_LINE
238 : 0, nblyr, & ! LCOV_EXCL_LINE
239 : hinc_old, hinc_old, & ! LCOV_EXCL_LINE
240 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
241 0 : bgrid(2:nblyr+1), surface_S )
242 0 : if (icepack_warnings_aborted(subname)) return
243 :
244 : call remap_zbgc(nilyr, &
245 : nt_qice, & ! LCOV_EXCL_LINE
246 : trcrn, trtmp_q, & ! LCOV_EXCL_LINE
247 : 0, nblyr, & ! LCOV_EXCL_LINE
248 : hinc_old, hinc_old, & ! LCOV_EXCL_LINE
249 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
250 0 : bgrid(2:nblyr+1), surface_S )
251 0 : if (icepack_warnings_aborted(subname)) return
252 :
253 0 : do k = 1, nblyr
254 0 : bqin (k+1) = min(c0, trtmp_q(nt_qice+k-1))
255 0 : bSin (k+1) = max(Smin, trtmp_s(nt_sice+k-1))
256 0 : bTin (k+1) = icepack_mushy_temperature_mush(bqin(k+1), bSin(k+1))
257 0 : bphin(k+1) = icepack_mushy_liquid_fraction (bTin(k+1), bSin(k+1))
258 : enddo ! k
259 :
260 0 : bSin (1) = bSin(2)
261 0 : bTin (1) = bTin(2)
262 0 : bphin(1) = bphin(2)
263 0 : bphin(nblyr+2) = c1
264 0 : bSin (nblyr+2) = sss
265 0 : bTin (nblyr+2) = sst
266 0 : bphin(nblyr+2) = c1
267 :
268 : !-----------------------------------------------------------------
269 : ! Define ice multiphase structure
270 : !-----------------------------------------------------------------
271 :
272 : call prepare_hbrine (nblyr, &
273 : bSin, bTin, iTin, & ! LCOV_EXCL_LINE
274 : brine_sal, brine_rho, & ! LCOV_EXCL_LINE
275 : ibrine_sal, ibrine_rho, & ! LCOV_EXCL_LINE
276 : sice_rho, & ! LCOV_EXCL_LINE
277 : bphin, iphin, & ! LCOV_EXCL_LINE
278 : kperm, bphi_min, & ! LCOV_EXCL_LINE
279 0 : igrid, sss)
280 0 : if (icepack_warnings_aborted(subname)) return
281 :
282 : call calculate_drho(nblyr, igrid, bgrid, &
283 0 : brine_rho, ibrine_rho, drho)
284 0 : if (icepack_warnings_aborted(subname)) return
285 :
286 0 : do k= 2, nblyr+1
287 0 : ikin(k) = k_o*iphin(k)**exp_h
288 0 : if (hbr_old .GT. puny) iDin(k) = iphin(k)*Dm/hbr_old**2
289 0 : if (hbr_old .GE. Ra_c) &
290 : iDin(k) = iDin(k) & ! LCOV_EXCL_LINE
291 0 : + l_sk*ikin(k)*gravit/viscos_dynamic*drho(k)/hbr_old**2
292 : enddo ! k
293 :
294 : end subroutine compute_microS_mushy
295 :
296 : !=======================================================================
297 :
298 0 : subroutine prepare_hbrine (nblyr, &
299 : bSin, bTin, iTin, & ! LCOV_EXCL_LINE
300 : brine_sal, brine_rho, & ! LCOV_EXCL_LINE
301 : ibrine_sal, ibrine_rho, & ! LCOV_EXCL_LINE
302 : sice_rho, bphin, iphin,& ! LCOV_EXCL_LINE
303 : kperm, bphi_min, & ! LCOV_EXCL_LINE
304 0 : i_grid, sss)
305 :
306 : integer (kind=int_kind), intent(in) :: &
307 : nblyr ! number of bio layers
308 :
309 : real (kind=dbl_kind), dimension (:), intent(in) :: &
310 : bSin , & ! salinity of ice layers on bio grid (ppt) ! LCOV_EXCL_LINE
311 : bTin , & ! temperature of ice layers on bio grid for history (C) ! LCOV_EXCL_LINE
312 : i_grid ! biology grid interface points
313 :
314 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
315 : brine_sal , & ! equilibrium brine salinity (ppt) ! LCOV_EXCL_LINE
316 : brine_rho , & ! internal brine density (kg/m^3) ! LCOV_EXCL_LINE
317 : ibrine_rho , & ! brine density on interface (kg/m^3) ! LCOV_EXCL_LINE
318 : ibrine_sal , & ! brine salinity on interface (ppt) ! LCOV_EXCL_LINE
319 : iphin , & ! porosity on interface ! LCOV_EXCL_LINE
320 : iTin , & ! Temperature on interface ! LCOV_EXCL_LINE
321 : bphin ! porosity of layers
322 :
323 : real (kind=dbl_kind), intent(in) :: &
324 : sss ! sea surface salinity (ppt)
325 :
326 : real (kind=dbl_kind), intent(out) :: &
327 : kperm , & ! harmonic average permeability (m^2) ! LCOV_EXCL_LINE
328 : bphi_min ! minimum porosity
329 :
330 : real (kind=dbl_kind), intent(inout) :: &
331 : sice_rho ! avg sea ice density
332 :
333 : ! local variables
334 :
335 : real (kind=dbl_kind), dimension(nblyr+1) :: &
336 0 : kin ! permeability
337 :
338 : real (kind=dbl_kind) :: &
339 : k_min, ktemp, & ! LCOV_EXCL_LINE
340 0 : igrp, igrm, rigr ! grid finite differences
341 :
342 : integer (kind=int_kind) :: &
343 : k ! layer index
344 :
345 : character(len=*),parameter :: subname='(prepare_hbrine)'
346 :
347 : !-----------------------------------------------------------------
348 : ! calculate equilibrium brine density and gradients
349 : !-----------------------------------------------------------------
350 :
351 0 : sice_rho = c0
352 :
353 0 : do k = 1, nblyr+1
354 :
355 0 : if (k == 1) then
356 0 : igrm = 0
357 : else
358 0 : igrm = i_grid(k) - i_grid(k-1)
359 : endif
360 :
361 0 : brine_sal(k) = a1*bTin(k) &
362 : + a2*bTin(k)**2 & ! LCOV_EXCL_LINE
363 0 : + a3*bTin(k)**3
364 0 : brine_rho(k) = b1 + b2*brine_sal(k)
365 0 : bphin (k) = max(puny, bSin(k)*rhosi &
366 0 : / (brine_sal(k)*brine_rho(k)))
367 0 : bphin (k) = min(c1, bphin(k))
368 0 : kin (k) = k_o*bphin(k)**exp_h
369 0 : sice_rho = sice_rho + (rhoi*(c1-bphin(k)) &
370 0 : + brine_rho(k)*bphin(k))*igrm
371 : enddo ! k
372 :
373 0 : brine_sal (nblyr+2) = sss
374 0 : brine_rho (nblyr+2) = rhow
375 0 : bphin (nblyr+2) = c1
376 0 : ibrine_sal(1) = brine_sal (2)
377 0 : ibrine_sal(nblyr+1) = brine_sal (nblyr+2)
378 0 : ibrine_rho(1) = brine_rho (2)
379 0 : ibrine_rho(nblyr+1) = brine_rho (nblyr+2)
380 0 : iTin (1) = bTin(2)
381 0 : iTin (nblyr+1) = bTin(nblyr+1)
382 0 : iphin (1) = bphin (2)
383 0 : iphin (nblyr+1) = bphin (nblyr+1)
384 0 : k_min = MINVAL(kin(2:nblyr+1))
385 0 : kperm = c0 ! initialize
386 0 : ktemp = c0
387 0 : bphi_min = bphin (1)
388 : ! bphi_min = max(bphin(1),bSin(2)*rhosi/bphin(2) &
389 : ! / (brine_sal(1)*brine_rho(1))*phi_snow)
390 :
391 0 : do k = 2, nblyr
392 0 : if (k_min > c0) then
393 0 : ktemp = ktemp + c1/kin(k)
394 0 : kperm = k_min
395 : endif
396 :
397 0 : igrp = i_grid(k+1) - i_grid(k )
398 0 : igrm = i_grid(k ) - i_grid(k-1)
399 0 : rigr = c1 / (i_grid(k+1)-i_grid(k-1))
400 :
401 0 : ibrine_sal(k) = (brine_sal(k+1)*igrp + brine_sal(k)*igrm) * rigr
402 0 : ibrine_rho(k) = (brine_rho(k+1)*igrp + brine_rho(k)*igrm) * rigr
403 0 : iTin (k) = (bTin (k+1)*igrp + bTin (k)*igrm) * rigr
404 0 : iphin (k) = max(puny, &
405 0 : (bphin (k+1)*igrp + bphin (k)*igrm) * rigr)
406 0 : iphin (k) = min(c1, iphin (k))
407 : enddo ! k
408 :
409 0 : if (k_min > c0) then
410 0 : ktemp = ktemp + c1/kin(nblyr+1)
411 0 : kperm = real(nblyr,kind=dbl_kind)/ktemp
412 : endif
413 :
414 0 : end subroutine prepare_hbrine
415 :
416 : !=======================================================================
417 :
418 : ! Changes include brine height increases from ice and snow surface melt,
419 : ! congelation growth, and upward pressure driven flow from snow loading.
420 : !
421 : ! Decreases arise from downward flushing and bottom melt.
422 : !
423 : ! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice
424 : ! with dynamic salinity or the height ratio == hbr/vicen*aicen, where
425 : ! hbr is the height of the brine surface relative to the bottom of the
426 : ! ice. This volume fraction may be > 1 in which case there is brine
427 : ! above the ice surface (ponds).
428 :
429 0 : subroutine update_hbrine (meltt, &
430 : melts, dt, & ! LCOV_EXCL_LINE
431 : hin, hsn, & ! LCOV_EXCL_LINE
432 : hin_old, hbr, & ! LCOV_EXCL_LINE
433 : hbr_old, & ! LCOV_EXCL_LINE
434 : fbri, & ! LCOV_EXCL_LINE
435 : dhS_top, dhS_bottom, & ! LCOV_EXCL_LINE
436 : dh_top_chl, dh_bot_chl, & ! LCOV_EXCL_LINE
437 : kperm, bphi_min, & ! LCOV_EXCL_LINE
438 : darcy_V, darcy_V_chl, & ! LCOV_EXCL_LINE
439 : bphin, aice0, & ! LCOV_EXCL_LINE
440 : dh_direct)
441 :
442 : real (kind=dbl_kind), intent(in) :: &
443 : dt ! timestep
444 :
445 : real (kind=dbl_kind), intent(in):: &
446 : meltt, & ! true top melt over dt (m) ! LCOV_EXCL_LINE
447 : melts, & ! true snow melt over dt (m) ! LCOV_EXCL_LINE
448 : hin, & ! ice thickness (m) ! LCOV_EXCL_LINE
449 : hsn, & ! snow thickness (m) ! LCOV_EXCL_LINE
450 : hin_old, & ! past timestep ice thickness (m) ! LCOV_EXCL_LINE
451 : hbr_old, & ! previous timestep hbr ! LCOV_EXCL_LINE
452 : kperm, & ! avg ice permeability ! LCOV_EXCL_LINE
453 : bphin, & ! upper brine porosity ! LCOV_EXCL_LINE
454 : dhS_bottom, & ! change in bottom hbr initially before darcy flow ! LCOV_EXCL_LINE
455 : aice0 ! open water area fraction
456 :
457 : real (kind=dbl_kind), intent(inout):: &
458 : darcy_V , & ! Darcy velocity: m/s ! LCOV_EXCL_LINE
459 : darcy_V_chl, & ! Darcy velocity: m/s for bgc ! LCOV_EXCL_LINE
460 : dhS_top , & ! change in top hbr before darcy flow ! LCOV_EXCL_LINE
461 : dh_bot_chl , & ! change in bottom for algae ! LCOV_EXCL_LINE
462 : dh_top_chl , & ! change in bottom for algae ! LCOV_EXCL_LINE
463 : hbr , & ! thickness of brine (m) ! LCOV_EXCL_LINE
464 : fbri , & ! brine height ratio tracer (hbr/hin) ! LCOV_EXCL_LINE
465 : bphi_min ! surface porosity
466 :
467 : real (kind=dbl_kind), intent(out):: &
468 : dh_direct ! surface flooding or runoff (m)
469 :
470 : ! local variables
471 :
472 : real (kind=dbl_kind) :: &
473 : hbrmin , & ! thinS or hin ! LCOV_EXCL_LINE
474 : dhbr_hin , & ! hbr-hin ! LCOV_EXCL_LINE
475 : hbrocn , & ! brine height above sea level (m) hbr-h_ocn ! LCOV_EXCL_LINE
476 : dhbr , & ! change in brine surface ! LCOV_EXCL_LINE
477 : h_ocn , & ! new ocean surface from ice bottom (m) ! LCOV_EXCL_LINE
478 : darcy_coeff, & ! magnitude of the Darcy velocity/hbrocn (1/s) ! LCOV_EXCL_LINE
479 : hbrocn_new , & ! hbrocn after flushing ! LCOV_EXCL_LINE
480 : dhflood , & ! surface flooding by ocean ! LCOV_EXCL_LINE
481 : exp_arg , & ! temporary exp value ! LCOV_EXCL_LINE
482 0 : dhrunoff ! direct runoff to ocean
483 :
484 : real (kind=dbl_kind), parameter :: &
485 : dh_min = p001 ! brine remains within dh_min of sea level
486 : ! when ice thickness is less than thinS
487 :
488 : character(len=*),parameter :: subname='(update_hbrine)'
489 :
490 0 : hbrocn = c0
491 0 : darcy_V = c0
492 0 : darcy_V_chl = c0
493 0 : hbrocn_new = c0
494 0 : h_ocn = rhosi/rhow*hin + rhos/rhow*hsn
495 0 : dh_direct = c0
496 :
497 0 : if (hbr_old > thinS .AND. hin_old > thinS .AND. hin > thinS ) then
498 0 : hbrmin = thinS
499 0 : dhS_top = -max(c0, min(hin_old-hbr_old, meltt)) * rhoi/rhow
500 0 : dhS_top = dhS_top - max(c0, melts) * rhos/rhow
501 0 : dh_top_chl = dhS_top
502 0 : dhbr = dhS_bottom - dhS_top
503 0 : hbr = max(puny, hbr_old+dhbr)
504 0 : hbrocn = hbr - h_ocn
505 0 : darcy_coeff = max(c0, kperm*gravit/(viscos*hbr_old))
506 :
507 0 : if (hbrocn > c0 .AND. hbr > thinS ) then
508 0 : bphi_min = bphin
509 0 : dhrunoff = -dhS_top*aice0
510 0 : hbrocn = max(c0,hbrocn - dhrunoff)
511 0 : exp_arg = darcy_coeff/bphi_min*dt
512 : ! tcraig avoids underflows
513 0 : if (exp_arg > exp_argmax) then
514 0 : hbrocn_new = c0
515 : else
516 0 : hbrocn_new = hbrocn*exp(-exp_arg)
517 : endif
518 0 : hbr = max(hbrmin, h_ocn + hbrocn_new)
519 0 : hbrocn_new = hbr-h_ocn
520 0 : darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn)
521 0 : darcy_V_chl= darcy_V
522 0 : dhS_top = dhS_top - darcy_V*dt/bphi_min + dhrunoff
523 0 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min + dhrunoff
524 0 : dh_direct = dhrunoff
525 0 : elseif (hbrocn < c0 .AND. hbr > thinS) then
526 0 : exp_arg = darcy_coeff/bphi_min*dt
527 : ! tcraig avoids underflows
528 0 : if (exp_arg > exp_argmax) then
529 0 : hbrocn_new = c0
530 : else
531 0 : hbrocn_new = hbrocn*exp(-exp_arg)
532 : endif
533 0 : dhflood = max(c0,hbrocn_new - hbrocn)*aice0
534 0 : hbr = max(hbrmin, h_ocn + hbrocn_new)
535 0 : darcy_V = -SIGN((hbrocn-hbrocn_new + dhflood)/dt*bphi_min, hbrocn)
536 0 : darcy_V_chl= darcy_V
537 0 : dhS_top = dhS_top - darcy_V*dt/bphi_min - dhflood
538 0 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min - dhflood
539 0 : dh_direct = -dhflood
540 : endif
541 :
542 0 : dh_bot_chl = dhS_bottom
543 :
544 : else ! very thin brine height
545 0 : hbrmin = min(thinS, hin)
546 0 : hbr = max(hbrmin, hbr_old+dhS_bottom-dhS_top)
547 0 : dhbr_hin = hbr - h_ocn
548 0 : if (abs(dhbr_hin) > dh_min) &
549 0 : hbr = max(hbrmin, h_ocn + SIGN(dh_min,dhbr_hin))
550 0 : dhS_top = hbr_old - hbr + dhS_bottom
551 0 : dh_top_chl = dhS_top
552 0 : dh_bot_chl = dhS_bottom
553 : endif
554 :
555 0 : fbri = hbr/hin
556 :
557 0 : end subroutine update_hbrine
558 :
559 : !==========================================================================================
560 : !
561 : ! Find density difference about interface grid points
562 : ! for gravity drainage parameterization
563 :
564 0 : subroutine calculate_drho (nblyr, i_grid, b_grid, &
565 0 : brine_rho, ibrine_rho, drho)
566 :
567 : integer (kind=int_kind), intent(in) :: &
568 : nblyr ! number of bio layers
569 :
570 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
571 : b_grid ! biology nondimensional grid layer points
572 :
573 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
574 : i_grid ! biology grid interface points
575 :
576 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
577 : brine_rho ! Internal brine density (kg/m^3)
578 :
579 : real (kind=dbl_kind), dimension (nblyr + 1), intent(in) :: &
580 : ibrine_rho ! Internal brine density (kg/m^3)
581 :
582 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
583 : drho ! brine difference about grid point (kg/m^3)
584 :
585 : ! local variables
586 :
587 : integer (kind=int_kind) :: &
588 : k, mm ! indices
589 :
590 : integer (kind=int_kind) :: &
591 : mstop, mstart
592 :
593 : real (kind=dbl_kind), dimension (nblyr + 1) :: & !on the zbgc vertical grid
594 : rho_a , & ! average brine density above grid point (kg/m^3) ! LCOV_EXCL_LINE
595 0 : rho_2a ! average brine density above and below grid points (kg/m^3)
596 :
597 : real (kind=dbl_kind), dimension (nblyr + 1) :: & !on the zbgc vertical grid
598 : rho_b , & ! brine density above grid point (kg/m^3) ! LCOV_EXCL_LINE
599 0 : rho_2b ! brine density above and below grid points (kg/m^3)
600 :
601 : character(len=*),parameter :: subname='(calculate_drho)'
602 :
603 0 : rho_a (:) = c0
604 0 : rho_2a(:) = c0
605 0 : rho_b (:) = c0
606 0 : rho_2b(:) = c0
607 0 : drho (:) = c0 ! surface is snow or atmosphere
608 :
609 0 : do k = 1, nblyr+1 ! i_grid values
610 :
611 : !----------------------------------------------
612 : ! h_avg(k) = i_grid(k)
613 : ! Calculate rho_a(k), ie average rho above i_grid(k)
614 : ! first part is good
615 : !----------------------------------------------
616 :
617 0 : if (k == 2) then
618 0 : rho_a(2) = (brine_rho(2)*b_grid(2) &
619 : + (ibrine_rho(2) + brine_rho(2)) & ! LCOV_EXCL_LINE
620 0 : * p5*(i_grid(2)-b_grid(2)) )/i_grid(2)
621 0 : rho_b(2) = brine_rho(2)
622 :
623 0 : elseif (k > 2 .AND. k < nblyr+1) then
624 0 : rho_a(k) = (rho_a(k-1)*i_grid(k-1) + (ibrine_rho(k-1) + brine_rho(k)) &
625 : * p5*(b_grid(k)-i_grid(k-1)) + (ibrine_rho(k ) + brine_rho(k)) & ! LCOV_EXCL_LINE
626 0 : * p5*(i_grid(k)-b_grid(k)))/i_grid(k)
627 0 : rho_b(k) = brine_rho(k)
628 : else
629 0 : rho_a(nblyr+1) = (rho_a(nblyr)*i_grid(nblyr) + (ibrine_rho(nblyr) + &
630 : brine_rho(nblyr+1))*p5*(b_grid(nblyr+1)-i_grid(nblyr)) + & ! LCOV_EXCL_LINE
631 0 : brine_rho(nblyr+1)*(i_grid(nblyr+1)-b_grid(nblyr+1)))/i_grid(nblyr+1)
632 0 : rho_a(1) = brine_rho(2) !for k == 1 use grid point value
633 0 : rho_b(nblyr+1) = brine_rho(nblyr+1)
634 0 : rho_b(1) = brine_rho(2)
635 : endif
636 :
637 : enddo !k
638 :
639 : !----------------------------------------------
640 : ! Calculate average above and below k rho_2a
641 : !----------------------------------------------
642 :
643 0 : do k = 1, nblyr+1 !i_grid values
644 0 : if (k == 1) then
645 0 : rho_2a(1) = (rho_a(1)*b_grid(2) + p5*(brine_rho(2) + ibrine_rho(2)) &
646 0 : * (i_grid(2)-b_grid(2)))/i_grid(2)
647 0 : rho_2b(1) = brine_rho(2)
648 : else
649 0 : mstop = 2*(k-1) + 1
650 0 : if (mstop < nblyr+1) then
651 0 : rho_2a(k) = rho_a(mstop)
652 0 : mstart = 2
653 0 : mstop = 1
654 : else
655 0 : mstart = nblyr+2
656 0 : mstop = nblyr+3
657 : endif
658 :
659 0 : do mm = mstart,mstop
660 0 : rho_2a(k) =(rho_a(nblyr+1) + rhow*(c2*i_grid(k)-c1))*p5/i_grid(k)
661 : enddo
662 0 : rho_2b(k) = brine_rho(k+1)
663 : endif
664 0 : drho(k) = max(rho_b(k) - rho_2b(k),max(c0,c2*(rho_a(k)-rho_2a(k)), &
665 0 : c2*(brine_rho(k)-brine_rho(k+1))/real(nblyr,kind=dbl_kind)))
666 : enddo
667 :
668 0 : end subroutine calculate_drho
669 :
670 : !=======================================================================
671 : !autodocument_start icepack_init_hbrine
672 : ! Initialize brine height tracer
673 :
674 0 : subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
675 0 : icgrid, swgrid, nblyr, nilyr, phi_snow)
676 :
677 : integer (kind=int_kind), intent(in) :: &
678 : nilyr, & ! number of ice layers ! LCOV_EXCL_LINE
679 : nblyr ! number of bio layers
680 :
681 : real (kind=dbl_kind), intent(inout) :: &
682 : phi_snow ! porosity at the ice-snow interface
683 :
684 : real (kind=dbl_kind), dimension (nblyr+2), intent(out) :: &
685 : bgrid ! biology nondimensional vertical grid points
686 :
687 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
688 : igrid ! biology vertical interface points
689 :
690 : real (kind=dbl_kind), dimension (nilyr+1), intent(out) :: &
691 : cgrid , & ! CICE vertical coordinate ! LCOV_EXCL_LINE
692 : icgrid , & ! interface grid for CICE (shortwave variable) ! LCOV_EXCL_LINE
693 : swgrid ! grid for ice tracers used in dEdd scheme
694 :
695 : !autodocument_end
696 :
697 : ! local variables
698 :
699 : integer (kind=int_kind) :: &
700 : k ! vertical index
701 :
702 : real (kind=dbl_kind) :: &
703 0 : zspace ! grid spacing for CICE vertical grid
704 :
705 : character(len=*),parameter :: subname='(icepack_init_hbrine)'
706 :
707 :
708 0 : if (phi_snow .le. c0) phi_snow = c1-rhos/rhoi
709 :
710 : !-----------------------------------------------------------------
711 : ! Calculate bio gridn: 0 to 1 corresponds to ice top to bottom
712 : !-----------------------------------------------------------------
713 :
714 0 : bgrid(:) = c0 ! zsalinity grid points
715 0 : bgrid(nblyr+2) = c1 ! bottom value
716 0 : igrid(:) = c0 ! bgc interface grid points
717 0 : igrid(1) = c0 ! ice top
718 0 : igrid(nblyr+1) = c1 ! ice bottom
719 :
720 0 : zspace = c1/max(c1,(real(nblyr,kind=dbl_kind)))
721 0 : do k = 2, nblyr+1
722 0 : bgrid(k) = zspace*(real(k,kind=dbl_kind) - c1p5)
723 : enddo
724 :
725 0 : do k = 2, nblyr
726 0 : igrid(k) = p5*(bgrid(k+1)+bgrid(k))
727 : enddo
728 :
729 : !-----------------------------------------------------------------
730 : ! Calculate CICE cgrid for interpolation ice top (0) to bottom (1)
731 : !-----------------------------------------------------------------
732 :
733 0 : cgrid(1) = c0 ! CICE vertical grid top point
734 0 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
735 :
736 0 : do k = 2, nilyr+1
737 0 : cgrid(k) = zspace * (real(k,kind=dbl_kind) - c1p5)
738 : enddo
739 :
740 : !-----------------------------------------------------------------
741 : ! Calculate CICE icgrid for ishortwave interpolation top(0) , bottom (1)
742 : !-----------------------------------------------------------------
743 :
744 0 : icgrid(1) = c0
745 0 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
746 :
747 0 : do k = 2, nilyr+1
748 0 : icgrid(k) = zspace * (real(k,kind=dbl_kind)-c1)
749 : enddo
750 :
751 : !------------------------------------------------------------------------
752 : ! Calculate CICE swgrid for dEdd ice: top of ice (0) , bottom of ice (1)
753 : ! Does not include snow
754 : ! see icepack_shortwave.F90
755 : ! swgrid represents the layer index of the delta-eddington ice layer index
756 : !------------------------------------------------------------------------
757 0 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
758 0 : swgrid(1) = min(c1/60.0_dbl_kind, zspace/c2)
759 0 : swgrid(2) = zspace/c2 !+ swgrid(1)
760 0 : do k = 3, nilyr+1
761 0 : swgrid(k) = zspace * (real(k,kind=dbl_kind)-c1p5)
762 : enddo
763 :
764 0 : end subroutine icepack_init_hbrine
765 :
766 : !=======================================================================
767 : !autodocument_start icepack_init_zsalinity
768 : ! **DEPRECATED**, all code removed
769 : ! Interface provided for backwards compatibility
770 :
771 0 : subroutine icepack_init_zsalinity(nblyr,ntrcr_o, Rayleigh_criteria, &
772 : Rayleigh_real, trcrn_bgc, nt_bgc_S, ncat, sss)
773 :
774 : integer (kind=int_kind), intent(in) :: &
775 : nblyr , & ! number of biolayers ! LCOV_EXCL_LINE
776 : ntrcr_o, & ! number of non bio tracers ! LCOV_EXCL_LINE
777 : ncat , & ! number of categories ! LCOV_EXCL_LINE
778 : nt_bgc_S ! zsalinity index
779 :
780 : logical (kind=log_kind), intent(inout) :: &
781 : Rayleigh_criteria
782 :
783 : real (kind=dbl_kind), intent(inout):: &
784 : Rayleigh_real
785 :
786 : real (kind=dbl_kind), intent(in):: &
787 : sss
788 :
789 : real (kind=dbl_kind), dimension(:,:), intent(inout):: &
790 : trcrn_bgc ! bgc subset of trcrn
791 :
792 : !autodocument_end
793 :
794 : ! local variables
795 :
796 : character(len=*),parameter :: subname='(icepack_init_zsalinity)'
797 :
798 0 : call icepack_warnings_add(subname//' DEPRECATED, do not use')
799 : ! call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
800 :
801 0 : end subroutine icepack_init_zsalinity
802 : !=======================================================================
803 :
804 : end module icepack_brine
805 :
806 : !=======================================================================
|
