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, &
29 : update_hbrine, &
30 : calculate_drho, &
31 : icepack_init_hbrine, &
32 : icepack_init_zsalinity ! deprecated
33 :
34 : real (kind=dbl_kind), parameter :: &
35 : maxhbr = 1.25_dbl_kind , & ! brine overflows if hbr > maxhbr*hin
36 : viscos = 2.1e-6_dbl_kind, & ! kinematic viscosity (m^2/s)
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)
40 : a3 = -0.012_dbl_kind, & ! (psu/C^3)
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 660884 : subroutine preflushing_changes (aicen, vicen, vsnon, &
62 : meltb, meltt, congel, &
63 : snoice, hice_old, dhice, &
64 : fbri, dhbr_top, dhbr_bot, &
65 : hbr_old, hin, hsn)
66 :
67 : real (kind=dbl_kind), intent(in) :: &
68 : aicen , & ! concentration of ice
69 : vicen , & ! volume per unit area of ice (m)
70 : vsnon , & ! volume per unit area of snow (m)
71 : meltb , & ! bottom ice melt (m)
72 : meltt , & ! top ice melt (m)
73 : congel , & ! bottom ice growth (m)
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)
81 : hsn , & ! snow thickness (m)
82 : dhice ! change due to sublimation (<0)/condensation (>0) (m)
83 :
84 : real (kind=dbl_kind), intent(inout) :: &
85 : fbri , & ! trcrn(nt_fbri)
86 : dhbr_top , & ! brine change in top for diagnostics (m)
87 : dhbr_bot , & ! brine change in bottom for diagnostics (m)
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 660884 : 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 660884 : hin = vicen / aicen
111 660884 : hsn = vsnon / aicen
112 660884 : hin_old = max(c0, hin + meltb + meltt - congel - snoice)
113 660884 : dhice = hin_old - hice_old ! change due to subl/cond
114 660884 : dhbr_top = meltt - snoice - dhice
115 660884 : dhbr_bot = congel - meltb
116 :
117 660884 : hbr_old = fbri * hice_old
118 :
119 660884 : 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 660884 : subroutine compute_microS_mushy (trcrn, hice_old, hbr_old, &
128 660884 : sss, sst, bTin, &
129 660884 : iTin, bphin, &
130 : kperm, bphi_min, &
131 660884 : bSin, brine_sal, brine_rho, &
132 660884 : iphin, ibrine_rho, ibrine_sal, &
133 660884 : iDin, iSin)
134 :
135 : real (kind=dbl_kind), intent(in) :: &
136 : hice_old , & ! previous timestep ice height (m)
137 : sss , & ! ocean salinity (ppt)
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)
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
155 : ibrine_rho , & ! brine rho on interface
156 : ibrine_sal , & ! brine sal on interface
157 : iTin , & ! Temperature on the igrid (oC)
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
162 : brine_sal , & ! equilibrium brine salinity (ppt)
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
167 : bphin ! porosity on bgrid
168 :
169 : ! local variables
170 :
171 : real (kind=dbl_kind), dimension (nilyr) :: &
172 1321768 : cSin , & ! bulk salinity (ppt)
173 1321768 : cqin ! enthalpy ()
174 :
175 : real (kind=dbl_kind), dimension (nblyr+2) :: &
176 : zTin , & ! Temperature of ice layers on bgrid (C)
177 : zSin , & ! Salinity of ice layers on bgrid (C)
178 660884 : bqin ! enthalpy on the bgrid ()
179 :
180 : real (kind=dbl_kind), dimension (nblyr+1) :: &
181 1321768 : 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
188 : hinc_old , & ! mean ice thickness before current melt/growth (m)
189 : hbrc_old ! mean brine thickness before current melt/growth (m)
190 :
191 : real (kind=dbl_kind), dimension (ntrcr+2) :: & ! nblyr+2)
192 1321768 : trtmp_s , & ! temporary, remapped tracers
193 1321768 : trtmp_q ! temporary, remapped tracers
194 :
195 : real (kind=dbl_kind), dimension(nblyr+1) :: &
196 660884 : 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 115765150 : trtmp_s(:) = c0
208 115765150 : trtmp_q(:) = c0
209 5947956 : iDin(:) = c0
210 :
211 5287072 : do k = 1, nilyr
212 4626188 : cSin(k) = trcrn(nt_sice+k-1)
213 5287072 : cqin(k) = trcrn(nt_qice+k-1)
214 : enddo
215 :
216 : ! map Sin and qin (cice) profiles to bgc grid
217 660884 : surface_S = min_salin
218 660884 : hinc_old = hice_old
219 :
220 : call remap_zbgc(nilyr, &
221 : nt_sice, &
222 : trcrn, trtmp_s, &
223 : 0, nblyr, &
224 : hinc_old, hinc_old, &
225 : cgrid(2:nilyr+1), &
226 660884 : bgrid(2:nblyr+1), surface_S )
227 660884 : if (icepack_warnings_aborted(subname)) return
228 :
229 : call remap_zbgc(nilyr, &
230 : nt_qice, &
231 : trcrn, trtmp_q, &
232 : 0, nblyr, &
233 : hinc_old, hinc_old, &
234 : cgrid(2:nilyr+1), &
235 660884 : bgrid(2:nblyr+1), surface_S )
236 660884 : if (icepack_warnings_aborted(subname)) return
237 :
238 5287072 : do k = 1, nblyr
239 4626188 : bqin (k+1) = min(c0, trtmp_q(nt_qice+k-1))
240 4626188 : bSin (k+1) = max(Smin, trtmp_s(nt_sice+k-1))
241 4626188 : bTin (k+1) = icepack_mushy_temperature_mush(bqin(k+1), bSin(k+1))
242 5287072 : bphin(k+1) = icepack_mushy_liquid_fraction (bTin(k+1), bSin(k+1))
243 : enddo ! k
244 :
245 660884 : bSin (1) = bSin(2)
246 660884 : bTin (1) = bTin(2)
247 660884 : bphin(1) = bphin(2)
248 660884 : bphin(nblyr+2) = c1
249 660884 : bSin (nblyr+2) = sss
250 660884 : bTin (nblyr+2) = sst
251 660884 : 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, &
259 : ibrine_sal, ibrine_rho, &
260 : bphin, iphin, &
261 : kperm, bphi_min, &
262 660884 : sss, iSin)
263 660884 : if (icepack_warnings_aborted(subname)) return
264 :
265 660884 : call calculate_drho(brine_rho, ibrine_rho, drho)
266 660884 : if (icepack_warnings_aborted(subname)) return
267 :
268 5287072 : do k= 2, nblyr+1
269 4626188 : ikin(k) = k_o*iphin(k)**exp_h
270 4626188 : if (hbr_old .GT. puny) iDin(k) = iphin(k)*Dm/hbr_old**2
271 4626188 : if (hbr_old .GE. Ra_c) &
272 : iDin(k) = iDin(k) &
273 5278833 : + 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 1321768 : brine_sal, brine_rho, &
282 660884 : ibrine_sal, ibrine_rho, &
283 1321768 : bphin, iphin, &
284 : kperm, bphi_min, &
285 660884 : sss, iSin)
286 :
287 : real (kind=dbl_kind), dimension (:), intent(in) :: &
288 : bSin , & ! salinity of ice layers on bio grid (ppt)
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)
293 : brine_rho , & ! internal brine density (kg/m^3)
294 : ibrine_rho , & ! brine density on interface (kg/m^3)
295 : ibrine_sal , & ! brine salinity on interface (ppt)
296 : iphin , & ! porosity on interface
297 : iTin , & ! Temperature on interface
298 : bphin , & ! porosity of layers
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)
306 : bphi_min ! minimum porosity
307 :
308 : ! local variables
309 :
310 : real (kind=dbl_kind), dimension(nblyr+1) :: &
311 1321768 : kin ! permeability
312 :
313 : real (kind=dbl_kind) :: &
314 : k_min, ktemp, &
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 5947956 : do k = 1, nblyr+1
327 :
328 5287072 : if (k == 1) then
329 660884 : igrm = 0
330 : else
331 4626188 : igrm = igrid(k) - igrid(k-1)
332 : endif
333 :
334 : brine_sal(k) = a1*bTin(k) &
335 : + a2*bTin(k)**2 &
336 5287072 : + a3*bTin(k)**3
337 5287072 : brine_rho(k) = b1 + b2*brine_sal(k)
338 : bphin (k) = max(puny, bSin(k)*rhosi &
339 5287072 : / (brine_sal(k)*brine_rho(k)))
340 5287072 : bphin (k) = min(c1, bphin(k))
341 5947956 : kin (k) = k_o*bphin(k)**exp_h
342 : enddo ! k
343 :
344 660884 : brine_sal (nblyr+2) = sss
345 660884 : brine_rho (nblyr+2) = rhow
346 660884 : bphin (nblyr+2) = c1
347 660884 : ibrine_sal(1) = brine_sal (2)
348 660884 : ibrine_sal(nblyr+1) = brine_sal (nblyr+2)
349 660884 : ibrine_rho(1) = brine_rho (2)
350 660884 : ibrine_rho(nblyr+1) = brine_rho (nblyr+2)
351 660884 : iTin (1) = bTin(2)
352 660884 : iTin (nblyr+1) = bTin(nblyr+1)
353 660884 : iSin (1) = bSin(2)
354 660884 : iSin (nblyr+1) = bSin(nblyr+1)
355 660884 : iphin (1) = bphin (2)
356 660884 : iphin (nblyr+1) = bphin (nblyr+1)
357 5947956 : k_min = MINVAL(kin(2:nblyr+1))
358 660884 : kperm = c0 ! initialize
359 660884 : ktemp = c0
360 660884 : 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 4626188 : do k = 2, nblyr
365 3965304 : if (k_min > c0) then
366 3965304 : ktemp = ktemp + c1/kin(k)
367 3965304 : kperm = k_min
368 : endif
369 :
370 3965304 : igrp = igrid(k+1) - igrid(k )
371 3965304 : igrm = igrid(k ) - igrid(k-1)
372 3965304 : rigr = c1 / (igrid(k+1)-igrid(k-1))
373 :
374 3965304 : ibrine_sal(k) = (brine_sal(k+1)*igrp + brine_sal(k)*igrm) * rigr
375 3965304 : ibrine_rho(k) = (brine_rho(k+1)*igrp + brine_rho(k)*igrm) * rigr
376 3965304 : iTin (k) = (bTin (k+1)*igrp + bTin (k)*igrm) * rigr
377 3965304 : iSin (k) = (bSin (k+1)*igrp + bSin (k)*igrm) * rigr
378 : iphin (k) = max(puny, &
379 3965304 : (bphin (k+1)*igrp + bphin (k)*igrm) * rigr)
380 4626188 : iphin (k) = min(c1, iphin (k))
381 : enddo ! k
382 :
383 660884 : if (k_min > c0) then
384 660884 : ktemp = ktemp + c1/kin(nblyr+1)
385 660884 : kperm = real(nblyr,kind=dbl_kind)/ktemp
386 : endif
387 :
388 660884 : 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 660884 : subroutine update_hbrine (meltt, &
404 : melts, dt, &
405 : hin, hsn, &
406 : hin_old, hbr, &
407 : hbr_old, &
408 : fbri, &
409 : dhS_top, dhS_bottom, &
410 : dh_top_chl, dh_bot_chl, &
411 : kperm, bphi_min, &
412 : darcy_V, darcy_V_chl, &
413 : bphin, aice0, &
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)
421 : melts, & ! true snow melt over dt (m)
422 : hin, & ! ice thickness (m)
423 : hsn, & ! snow thickness (m)
424 : hin_old, & ! past timestep ice thickness (m)
425 : hbr_old, & ! previous timestep hbr
426 : kperm, & ! avg ice permeability
427 : bphin, & ! upper brine porosity
428 : dhS_bottom, & ! change in bottom hbr initially before darcy flow
429 : aice0 ! open water area fraction
430 :
431 : real (kind=dbl_kind), intent(inout):: &
432 : darcy_V , & ! Darcy velocity: m/s
433 : darcy_V_chl, & ! Darcy velocity: m/s for bgc
434 : dhS_top , & ! change in top hbr before darcy flow
435 : dh_bot_chl , & ! change in bottom for algae
436 : dh_top_chl , & ! change in bottom for algae
437 : hbr , & ! thickness of brine (m)
438 : fbri , & ! brine height ratio tracer (hbr/hin)
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
448 : dhbr_hin , & ! hbr-hin
449 : hbrocn , & ! brine height above sea level (m) hbr-h_ocn
450 : dhbr , & ! change in brine surface
451 : h_ocn , & ! new ocean surface from ice bottom (m)
452 : darcy_coeff, & ! magnitude of the Darcy velocity/hbrocn (1/s)
453 : hbrocn_new , & ! hbrocn after flushing
454 : dhflood , & ! surface flooding by ocean
455 : exp_arg , & ! temporary exp value
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 660884 : hbrocn = c0
465 660884 : darcy_V = c0
466 660884 : darcy_V_chl = c0
467 660884 : hbrocn_new = c0
468 660884 : h_ocn = rhosi/rhow*hin + rhos/rhow*hsn
469 660884 : dh_direct = c0
470 :
471 660884 : if (hbr_old > thinS .AND. hin_old > thinS .AND. hin > thinS ) then
472 659684 : hbrmin = thinS
473 659684 : dhS_top = -max(c0, min(hin_old-hbr_old, meltt)) * rhoi/rhow
474 659684 : dhS_top = dhS_top - max(c0, melts) * rhos/rhow
475 659684 : dh_top_chl = dhS_top
476 659684 : dhbr = dhS_bottom - dhS_top
477 659684 : hbr = max(puny, hbr_old+dhbr)
478 659684 : hbrocn = hbr - h_ocn
479 659684 : darcy_coeff = max(c0, kperm*gravit/(viscos*hbr_old))
480 :
481 659684 : if (hbrocn > c0 .AND. hbr > thinS ) then
482 432411 : bphi_min = bphin
483 432411 : dhrunoff = -dhS_top*aice0
484 432411 : hbrocn = max(c0,hbrocn - dhrunoff)
485 432411 : exp_arg = darcy_coeff/bphi_min*dt
486 : ! tcraig avoids underflows
487 432411 : if (exp_arg > exp_argmax) then
488 6918 : hbrocn_new = c0
489 : else
490 425493 : hbrocn_new = hbrocn*exp(-exp_arg)
491 : endif
492 432411 : hbr = max(hbrmin, h_ocn + hbrocn_new)
493 432411 : hbrocn_new = hbr-h_ocn
494 432411 : darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn)
495 432411 : darcy_V_chl= darcy_V
496 432411 : dhS_top = dhS_top - darcy_V*dt/bphi_min + dhrunoff
497 432411 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min + dhrunoff
498 432411 : dh_direct = dhrunoff
499 227273 : elseif (hbrocn < c0 .AND. hbr > thinS) then
500 227233 : exp_arg = darcy_coeff/bphi_min*dt
501 : ! tcraig avoids underflows
502 227233 : if (exp_arg > exp_argmax) then
503 2235 : hbrocn_new = c0
504 : else
505 224998 : hbrocn_new = hbrocn*exp(-exp_arg)
506 : endif
507 227233 : dhflood = max(c0,hbrocn_new - hbrocn)*aice0
508 227233 : hbr = max(hbrmin, h_ocn + hbrocn_new)
509 227233 : darcy_V = -SIGN((hbrocn-hbrocn_new + dhflood)/dt*bphi_min, hbrocn)
510 227233 : darcy_V_chl= darcy_V
511 227233 : dhS_top = dhS_top - darcy_V*dt/bphi_min - dhflood
512 227233 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min - dhflood
513 227233 : dh_direct = -dhflood
514 : endif
515 :
516 659684 : dh_bot_chl = dhS_bottom
517 :
518 : else ! very thin brine height
519 1200 : hbrmin = min(thinS, hin)
520 1200 : hbr = max(hbrmin, hbr_old+dhS_bottom-dhS_top)
521 1200 : dhbr_hin = hbr - h_ocn
522 1200 : if (abs(dhbr_hin) > dh_min) &
523 531 : hbr = max(hbrmin, h_ocn + SIGN(dh_min,dhbr_hin))
524 1200 : dhS_top = hbr_old - hbr + dhS_bottom
525 1200 : dh_top_chl = dhS_top
526 1200 : dh_bot_chl = dhS_bottom
527 : endif
528 :
529 660884 : fbri = hbr/hin
530 :
531 660884 : end subroutine update_hbrine
532 :
533 : !==========================================================================================
534 : !
535 : ! Find density difference about interface grid points
536 : ! for gravity drainage parameterization
537 :
538 660884 : 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 1321768 : rho_a , & ! average brine density above grid point (kg/m^3)
559 1321768 : 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 1321768 : rho_b , & ! brine density above grid point (kg/m^3)
563 1321768 : rho_2b ! brine density above and below grid points (kg/m^3)
564 :
565 : character(len=*),parameter :: subname='(calculate_drho)'
566 :
567 5947956 : rho_a (:) = c0
568 5947956 : rho_2a(:) = c0
569 5947956 : rho_b (:) = c0
570 5947956 : rho_2b(:) = c0
571 5947956 : drho (:) = c0 ! surface is snow or atmosphere
572 :
573 5947956 : 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 5947956 : if (k == 2) then
582 : rho_a(2) = (brine_rho(2)*bgrid(2) &
583 : + (ibrine_rho(2) + brine_rho(2)) &
584 660884 : * p5*(igrid(2)-bgrid(2)) )/igrid(2)
585 660884 : rho_b(2) = brine_rho(2)
586 :
587 4626188 : 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)) &
590 3304420 : * p5*(igrid(k)-bgrid(k)))/igrid(k)
591 3304420 : 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)) + &
595 1321768 : brine_rho(nblyr+1)*(igrid(nblyr+1)-bgrid(nblyr+1)))/igrid(nblyr+1)
596 1321768 : rho_a(1) = brine_rho(2) !for k == 1 use grid point value
597 1321768 : rho_b(nblyr+1) = brine_rho(nblyr+1)
598 1321768 : 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 5947956 : do k = 1, nblyr+1 !igrid values
608 5287072 : if (k == 1) then
609 : rho_2a(1) = (rho_a(1)*bgrid(2) + p5*(brine_rho(2) + ibrine_rho(2)) &
610 660884 : * (igrid(2)-bgrid(2)))/igrid(2)
611 660884 : rho_2b(1) = brine_rho(2)
612 : else
613 4626188 : mstop = 2*(k-1) + 1
614 4626188 : if (mstop < nblyr+1) then
615 1982652 : rho_2a(k) = rho_a(mstop)
616 1982652 : mstart = 2
617 1982652 : mstop = 1
618 : else
619 2643536 : mstart = nblyr+2
620 2643536 : mstop = nblyr+3
621 : endif
622 :
623 9913260 : do mm = mstart,mstop
624 9913260 : rho_2a(k) =(rho_a(nblyr+1) + rhow*(c2*igrid(k)-c1))*p5/igrid(k)
625 : enddo
626 4626188 : 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 5947956 : c2*(brine_rho(k)-brine_rho(k+1))/real(nblyr,kind=dbl_kind)))
630 : enddo
631 :
632 660884 : end subroutine calculate_drho
633 :
634 : !=======================================================================
635 : !autodocument_start icepack_init_hbrine
636 : ! Initialize brine height tracer
637 :
638 9 : subroutine icepack_init_hbrine(bgrid_out, igrid_out, cgrid_out, &
639 9 : 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
652 : icgrid_out , & ! interface grid for CICE (shortwave variable)
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 9 : if (present(phi_snow)) then
670 9 : if (phi_snow .le. c0) phi_snow = c1-rhos/rhoi
671 : endif
672 :
673 9 : allocate(bgrid (nblyr+2))
674 9 : allocate(igrid (nblyr+1))
675 9 : allocate(cgrid (nilyr+1))
676 9 : allocate(icgrid(nilyr+1))
677 9 : allocate(swgrid(nilyr+1))
678 :
679 : !-----------------------------------------------------------------
680 : ! Calculate bio gridn: 0 to 1 corresponds to ice top to bottom
681 : !-----------------------------------------------------------------
682 :
683 90 : bgrid(:) = c0 ! zsalinity grid points
684 9 : bgrid(nblyr+2) = c1 ! bottom value
685 81 : igrid(:) = c0 ! bgc interface grid points
686 9 : igrid(1) = c0 ! ice top
687 9 : igrid(nblyr+1) = c1 ! ice bottom
688 :
689 9 : zspace = c1/max(c1,(real(nblyr,kind=dbl_kind)))
690 72 : do k = 2, nblyr+1
691 72 : bgrid(k) = zspace*(real(k,kind=dbl_kind) - c1p5)
692 : enddo
693 :
694 63 : do k = 2, nblyr
695 63 : 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 9 : cgrid(1) = c0 ! CICE vertical grid top point
703 9 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
704 :
705 72 : do k = 2, nilyr+1
706 72 : 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 9 : icgrid(1) = c0
714 9 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
715 :
716 72 : do k = 2, nilyr+1
717 72 : 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 9 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
727 9 : swgrid(1) = min(c1/60.0_dbl_kind, zspace*p25) !p5 to p25. NJ: allows thinner surface layers
728 9 : swgrid(2) = zspace/c2 !+ swgrid(1)
729 63 : do k = 3, nilyr+1
730 63 : swgrid(k) = zspace * (real(k,kind=dbl_kind)-c1p5)
731 : enddo
732 :
733 9 : if (present( bgrid_out)) bgrid_out=bgrid
734 9 : if (present( cgrid_out)) cgrid_out=cgrid
735 9 : if (present( igrid_out)) igrid_out=igrid
736 9 : if (present(icgrid_out)) icgrid_out=icgrid
737 9 : if (present(swgrid_out)) swgrid_out=swgrid
738 :
739 9 : 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 : !=======================================================================
|
