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