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, & ! LCOV_EXCL_LINE
29 : update_hbrine, & ! LCOV_EXCL_LINE
30 : compute_microS, & ! LCOV_EXCL_LINE
31 : calculate_drho, & ! LCOV_EXCL_LINE
32 : icepack_init_hbrine, & ! LCOV_EXCL_LINE
33 : icepack_init_zsalinity
34 :
35 : real (kind=dbl_kind), parameter :: &
36 : maxhbr = 1.25_dbl_kind , & ! brine overflows if hbr > maxhbr*hin ! LCOV_EXCL_LINE
37 : viscos = 2.1e-6_dbl_kind, & ! kinematic viscosity (m^2/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
41 : a3 = -0.012_dbl_kind, & ! (psu/C^3) ! LCOV_EXCL_LINE
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 90745709 : subroutine preflushing_changes (aicen, vicen, vsnon, &
63 : meltb, meltt, congel, & ! LCOV_EXCL_LINE
64 : snoice, hice_old, dhice, & ! LCOV_EXCL_LINE
65 : fbri, dhbr_top, dhbr_bot, & ! LCOV_EXCL_LINE
66 : hbr_old, hin,hsn, firstice )
67 :
68 : real (kind=dbl_kind), intent(in) :: &
69 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
70 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
71 : vsnon , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
72 : meltb , & ! bottom ice melt (m) ! LCOV_EXCL_LINE
73 : meltt , & ! top ice melt (m) ! LCOV_EXCL_LINE
74 : congel , & ! bottom ice growth (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
82 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
83 : dhice ! change due to sublimation (<0)/condensation (>0) (m)
84 :
85 : real (kind=dbl_kind), intent(inout) :: &
86 : fbri , & ! trcrn(nt_fbri) ! LCOV_EXCL_LINE
87 : dhbr_top , & ! brine change in top for diagnostics (m) ! LCOV_EXCL_LINE
88 : dhbr_bot , & ! brine change in bottom for diagnostics (m) ! LCOV_EXCL_LINE
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 1220076 : hin_old ! ice thickness before current melt/growth (m)
98 :
99 : character(len=*),parameter :: subname='(preflushing_changes)'
100 :
101 : !-----------------------------------------------------------------
102 : ! initialize
103 : !-----------------------------------------------------------------
104 :
105 90745709 : 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 90745709 : hin = vicen / aicen
115 90745709 : hsn = vsnon / aicen
116 90745709 : hin_old = max(c0, hin + meltb + meltt - congel - snoice)
117 90745709 : dhice = hin_old - hice_old ! change due to subl/cond
118 90745709 : dhbr_top = meltt - snoice - dhice
119 90745709 : dhbr_bot = congel - meltb
120 :
121 90745709 : if ((hice_old < puny) .OR. (hin_old < puny) .OR. firstice) then
122 309674 : hice_old = hin
123 309674 : dhbr_top = c0
124 309674 : dhbr_bot = c0
125 309674 : dhice = c0
126 309674 : fbri = c1
127 : endif
128 :
129 90745709 : hbr_old = fbri * hice_old
130 :
131 90745709 : 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 86433963 : subroutine compute_microS_mushy (nilyr, nblyr, &
141 : bgrid, cgrid, igrid, & ! LCOV_EXCL_LINE
142 : trcrn, hice_old, hbr_old, & ! LCOV_EXCL_LINE
143 : sss, sst, bTin, & ! LCOV_EXCL_LINE
144 : iTin, bphin, & ! LCOV_EXCL_LINE
145 : kperm, bphi_min, & ! LCOV_EXCL_LINE
146 : bSin, brine_sal, brine_rho, & ! LCOV_EXCL_LINE
147 : iphin, ibrine_rho, ibrine_sal, & ! LCOV_EXCL_LINE
148 86433963 : sice_rho, iDin )
149 :
150 : integer (kind=int_kind), intent(in) :: &
151 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
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) :: & ! LCOV_EXCL_LINE
165 : hice_old , & ! previous timestep ice height (m) ! LCOV_EXCL_LINE
166 : sss , & ! ocean salinity (ppt) ! LCOV_EXCL_LINE
167 : sst ! ocean temperature (C)
168 :
169 : real (kind=dbl_kind), dimension(ntrcr), &
170 : intent(in) :: & ! LCOV_EXCL_LINE
171 : trcrn
172 :
173 : real (kind=dbl_kind), intent(out) :: &
174 : kperm , & ! average ice permeability (m^2) ! LCOV_EXCL_LINE
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) :: & ! LCOV_EXCL_LINE
182 : iDin ! tracer diffusivity/h^2 (1/s) includes gravity drainage/molecular
183 :
184 : real (kind=dbl_kind), dimension (nblyr+1), &
185 : intent(inout) :: & ! LCOV_EXCL_LINE
186 : iphin , & ! porosity on the igrid ! LCOV_EXCL_LINE
187 : ibrine_rho , & ! brine rho on interface ! LCOV_EXCL_LINE
188 : ibrine_sal , & ! brine sal on interface ! LCOV_EXCL_LINE
189 : iTin ! Temperature on the igrid (oC)
190 :
191 : real (kind=dbl_kind), dimension (nblyr+2), &
192 : intent(inout) :: & ! LCOV_EXCL_LINE
193 : bSin , & ! bulk salinity (ppt) on bgrid ! LCOV_EXCL_LINE
194 : brine_sal , & ! equilibrium brine salinity (ppt) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
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 174130068 : bqin ! enthalpy on the bgrid ()
206 :
207 : real (kind=dbl_kind), dimension (nblyr+1) :: &
208 173589150 : ikin ! permeability (m^2)
209 :
210 : integer (kind=int_kind) :: &
211 : k ! vertical biology layer index
212 :
213 : real (kind=dbl_kind) :: &
214 : surface_S , & ! salinity of ice above hin > hbr ! LCOV_EXCL_LINE
215 180306 : hinc_old ! mean ice thickness before current melt/growth (m)
216 :
217 : real (kind=dbl_kind), dimension (ntrcr+2) :: & ! nblyr+2)
218 : trtmp_s , & ! temporary, remapped tracers ! LCOV_EXCL_LINE
219 197840307 : trtmp_q ! temporary, remapped tracers
220 :
221 : real (kind=dbl_kind), dimension(nblyr+1) :: &
222 87696105 : 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 9121344243 : trtmp_s(:) = c0
234 9121344243 : trtmp_q(:) = c0
235 438077307 : iDin(:) = c0
236 :
237 : ! map Sin and qin (cice) profiles to bgc grid
238 86433963 : surface_S = min_salin
239 86433963 : hbr_old = min(hbr_old, maxhbr*hice_old)
240 86433963 : hinc_old = hice_old
241 :
242 : call remap_zbgc(nilyr, &
243 : nt_sice, & ! LCOV_EXCL_LINE
244 : trcrn, trtmp_s, & ! LCOV_EXCL_LINE
245 : 0, nblyr, & ! LCOV_EXCL_LINE
246 : hinc_old, hinc_old, & ! LCOV_EXCL_LINE
247 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
248 86433963 : bgrid(2:nblyr+1), surface_S )
249 86433963 : if (icepack_warnings_aborted(subname)) return
250 :
251 : call remap_zbgc(nilyr, &
252 : nt_qice, & ! LCOV_EXCL_LINE
253 : trcrn, trtmp_q, & ! LCOV_EXCL_LINE
254 : 0, nblyr, & ! LCOV_EXCL_LINE
255 : hinc_old, hinc_old, & ! LCOV_EXCL_LINE
256 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
257 86433963 : bgrid(2:nblyr+1), surface_S )
258 86433963 : if (icepack_warnings_aborted(subname)) return
259 :
260 351643344 : do k = 1, nblyr
261 265209381 : bqin (k+1) = min(c0, trtmp_q(nt_qice+k-1))
262 265209381 : bSin (k+1) = max(Smin, trtmp_s(nt_sice+k-1))
263 265209381 : bTin (k+1) = icepack_mushy_temperature_mush(bqin(k+1), bSin(k+1))
264 351643344 : bphin(k+1) = icepack_mushy_liquid_fraction (bTin(k+1), bSin(k+1))
265 : enddo ! k
266 :
267 86433963 : bSin (1) = bSin(2)
268 86433963 : bTin (1) = bTin(2)
269 86433963 : bphin(1) = bphin(2)
270 86433963 : bphin(nblyr+2) = c1
271 86433963 : bSin (nblyr+2) = sss
272 86433963 : bTin (nblyr+2) = sst
273 86433963 : bphin(nblyr+2) = c1
274 :
275 : !-----------------------------------------------------------------
276 : ! Define ice multiphase structure
277 : !-----------------------------------------------------------------
278 :
279 : call prepare_hbrine (nblyr, &
280 : bSin, bTin, iTin, & ! LCOV_EXCL_LINE
281 : brine_sal, brine_rho, & ! LCOV_EXCL_LINE
282 : ibrine_sal, ibrine_rho, & ! LCOV_EXCL_LINE
283 : sice_rho, & ! LCOV_EXCL_LINE
284 : bphin, iphin, & ! LCOV_EXCL_LINE
285 : kperm, bphi_min, & ! LCOV_EXCL_LINE
286 86433963 : igrid, sss)
287 86433963 : if (icepack_warnings_aborted(subname)) return
288 :
289 : call calculate_drho(nblyr, igrid, bgrid, &
290 86433963 : brine_rho, ibrine_rho, drho)
291 86433963 : if (icepack_warnings_aborted(subname)) return
292 :
293 351643344 : do k= 2, nblyr+1
294 265209381 : ikin(k) = k_o*iphin(k)**exp_h
295 265209381 : iDin(k) = iphin(k)*Dm/hbr_old**2
296 265209381 : if (hbr_old .GE. Ra_c) &
297 : iDin(k) = iDin(k) & ! LCOV_EXCL_LINE
298 345077537 : + 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 90745709 : subroutine prepare_hbrine (nblyr, &
306 : bSin, bTin, iTin, & ! LCOV_EXCL_LINE
307 : brine_sal, brine_rho, & ! LCOV_EXCL_LINE
308 : ibrine_sal, ibrine_rho, & ! LCOV_EXCL_LINE
309 : sice_rho, bphin, iphin,& ! LCOV_EXCL_LINE
310 : kperm, bphi_min, & ! LCOV_EXCL_LINE
311 90745709 : 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) :: & ! LCOV_EXCL_LINE
318 : bSin , & ! salinity of ice layers on bio grid (ppt) ! LCOV_EXCL_LINE
319 : bTin , & ! temperature of ice layers on bio grid for history (C) ! LCOV_EXCL_LINE
320 : i_grid ! biology grid interface points
321 :
322 : real (kind=dbl_kind), dimension (:), &
323 : intent(inout) :: & ! LCOV_EXCL_LINE
324 : brine_sal , & ! equilibrium brine salinity (ppt) ! LCOV_EXCL_LINE
325 : brine_rho , & ! internal brine density (kg/m^3) ! LCOV_EXCL_LINE
326 : ibrine_rho , & ! brine density on interface (kg/m^3) ! LCOV_EXCL_LINE
327 : ibrine_sal , & ! brine salinity on interface (ppt) ! LCOV_EXCL_LINE
328 : iphin , & ! porosity on interface ! LCOV_EXCL_LINE
329 : iTin , & ! Temperature on interface ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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 189491032 : kin ! permeability
346 :
347 : real (kind=dbl_kind) :: &
348 : k_min, ktemp, & ! LCOV_EXCL_LINE
349 1220076 : 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 90745709 : sice_rho = c0
361 :
362 476883021 : do k = 1, nblyr+1
363 :
364 386137312 : if (k == 1) then
365 90745709 : igrm = 0
366 : else
367 295391603 : igrm = i_grid(k) - i_grid(k-1)
368 : endif
369 :
370 9219690 : brine_sal(k) = a1*bTin(k) &
371 : + a2*bTin(k)**2 & ! LCOV_EXCL_LINE
372 386137312 : + a3*bTin(k)**3
373 386137312 : brine_rho(k) = b1 + b2*brine_sal(k)
374 9219690 : bphin (k) = max(puny, bSin(k)*rhosi &
375 386137312 : / (brine_sal(k)*brine_rho(k)))
376 386137312 : bphin (k) = min(c1, bphin(k))
377 386137312 : kin (k) = k_o*bphin(k)**exp_h
378 9219690 : sice_rho = sice_rho + (rhoi*(c1-bphin(k)) &
379 476883021 : + brine_rho(k)*bphin(k))*igrm
380 : enddo ! k
381 :
382 90745709 : brine_sal (nblyr+2) = sss
383 90745709 : brine_rho (nblyr+2) = rhow
384 90745709 : bphin (nblyr+2) = c1
385 90745709 : ibrine_sal(1) = brine_sal (2)
386 90745709 : ibrine_sal(nblyr+1) = brine_sal (nblyr+2)
387 90745709 : ibrine_rho(1) = brine_rho (2)
388 90745709 : ibrine_rho(nblyr+1) = brine_rho (nblyr+2)
389 90745709 : iTin (1) = bTin(2)
390 90745709 : iTin (nblyr+1) = bTin(nblyr+1)
391 90745709 : iphin (1) = bphin (2)
392 90745709 : iphin (nblyr+1) = bphin (nblyr+1)
393 476883021 : k_min = MINVAL(kin(2:nblyr+1))
394 90745709 : kperm = c0 ! initialize
395 90745709 : ktemp = c0
396 90745709 : 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 295391603 : do k = 2, nblyr
401 204645894 : if (k_min > c0) then
402 204645894 : ktemp = ktemp + c1/kin(k)
403 204645894 : kperm = k_min
404 : endif
405 :
406 204645894 : igrp = i_grid(k+1) - i_grid(k )
407 204645894 : igrm = i_grid(k ) - i_grid(k-1)
408 204645894 : rigr = c1 / (i_grid(k+1)-i_grid(k-1))
409 :
410 204645894 : ibrine_sal(k) = (brine_sal(k+1)*igrp + brine_sal(k)*igrm) * rigr
411 204645894 : ibrine_rho(k) = (brine_rho(k+1)*igrp + brine_rho(k)*igrm) * rigr
412 204645894 : iTin (k) = (bTin (k+1)*igrp + bTin (k)*igrm) * rigr
413 6779538 : iphin (k) = max(puny, &
414 204645894 : (bphin (k+1)*igrp + bphin (k)*igrm) * rigr)
415 295391603 : iphin (k) = min(c1, iphin (k))
416 : enddo ! k
417 :
418 90745709 : if (k_min > c0) then
419 90745709 : ktemp = ktemp + c1/kin(nblyr+1)
420 90745709 : kperm = real(nblyr,kind=dbl_kind)/ktemp
421 : endif
422 :
423 90745709 : 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 90745709 : subroutine update_hbrine (meltt, &
439 : melts, dt, & ! LCOV_EXCL_LINE
440 : hin, hsn, & ! LCOV_EXCL_LINE
441 : hin_old, hbr, & ! LCOV_EXCL_LINE
442 : hbr_old, & ! LCOV_EXCL_LINE
443 : fbri, & ! LCOV_EXCL_LINE
444 : dhS_top, dhS_bottom, & ! LCOV_EXCL_LINE
445 : dh_top_chl, dh_bot_chl, & ! LCOV_EXCL_LINE
446 : kperm, bphi_min, & ! LCOV_EXCL_LINE
447 : darcy_V, darcy_V_chl, & ! LCOV_EXCL_LINE
448 : bphin, aice0, & ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
456 : melts, & ! true snow melt over dt (m) ! LCOV_EXCL_LINE
457 : hin, & ! ice thickness (m) ! LCOV_EXCL_LINE
458 : hsn, & ! snow thickness (m) ! LCOV_EXCL_LINE
459 : hin_old, & ! past timestep ice thickness (m) ! LCOV_EXCL_LINE
460 : hbr_old, & ! previous timestep hbr ! LCOV_EXCL_LINE
461 : kperm, & ! avg ice permeability ! LCOV_EXCL_LINE
462 : bphin, & ! upper brine porosity ! LCOV_EXCL_LINE
463 : dhS_bottom, & ! change in bottom hbr initially before darcy flow ! LCOV_EXCL_LINE
464 : aice0 ! open water area fraction
465 :
466 : real (kind=dbl_kind), intent(inout):: &
467 : darcy_V , & ! Darcy velocity: m/s ! LCOV_EXCL_LINE
468 : darcy_V_chl, & ! Darcy velocity: m/s for bgc ! LCOV_EXCL_LINE
469 : dhS_top , & ! change in top hbr before darcy flow ! LCOV_EXCL_LINE
470 : dh_bot_chl , & ! change in bottom for algae ! LCOV_EXCL_LINE
471 : dh_top_chl , & ! change in bottom for algae ! LCOV_EXCL_LINE
472 : hbr , & ! thickness of brine (m) ! LCOV_EXCL_LINE
473 : fbri , & ! brine height ratio tracer (hbr/hin) ! LCOV_EXCL_LINE
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 : hbrmin , & ! thinS or hin ! LCOV_EXCL_LINE
483 : dhbr_hin , & ! hbr-hin ! LCOV_EXCL_LINE
484 : hbrocn , & ! brine height above sea level (m) hbr-h_ocn ! LCOV_EXCL_LINE
485 : dhbr , & ! change in brine surface ! LCOV_EXCL_LINE
486 : h_ocn , & ! new ocean surface from ice bottom (m) ! LCOV_EXCL_LINE
487 : darcy_coeff, & ! magnitude of the Darcy velocity/hbrocn (1/s) ! LCOV_EXCL_LINE
488 : hbrocn_new , & ! hbrocn after flushing ! LCOV_EXCL_LINE
489 : dhflood , & ! surface flooding by ocean ! LCOV_EXCL_LINE
490 : exp_arg , & ! temporary exp value ! LCOV_EXCL_LINE
491 1220076 : 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 90745709 : hbrocn = c0
500 90745709 : darcy_V = c0
501 90745709 : darcy_V_chl = c0
502 90745709 : hbrocn_new = c0
503 90745709 : h_ocn = rhosi/rhow*hin + rhos/rhow*hsn
504 90745709 : dh_direct = c0
505 :
506 90745709 : if (hbr_old > thinS .AND. hin_old > thinS .AND. hin > thinS ) then
507 88963691 : hbrmin = thinS
508 88963691 : dhS_top = -max(c0, min(hin_old-hbr_old, meltt)) * rhoi/rhow
509 88963691 : dhS_top = dhS_top - max(c0, melts) * rhos/rhow
510 88963691 : dh_top_chl = dhS_top
511 88963691 : dhbr = dhS_bottom - dhS_top
512 88963691 : hbr = max(puny, hbr_old+dhbr)
513 88963691 : hbrocn = hbr - h_ocn
514 88963691 : darcy_coeff = max(c0, kperm*gravit/(viscos*hbr_old))
515 :
516 88963691 : if (hbrocn > c0 .AND. hbr > thinS ) then
517 63799639 : bphi_min = bphin
518 63799639 : dhrunoff = -dhS_top*aice0
519 63799639 : hbrocn = max(c0,hbrocn - dhrunoff)
520 63799639 : exp_arg = darcy_coeff/bphi_min*dt
521 : ! tcraig avoids underflows
522 63799639 : if (exp_arg > exp_argmax) then
523 14737056 : hbrocn_new = c0
524 : else
525 49062583 : hbrocn_new = hbrocn*exp(-exp_arg)
526 : endif
527 63799639 : hbr = max(hbrmin, h_ocn + hbrocn_new)
528 63799639 : hbrocn_new = hbr-h_ocn
529 63799639 : darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn)
530 63799639 : darcy_V_chl= darcy_V
531 63799639 : dhS_top = dhS_top - darcy_V*dt/bphi_min + dhrunoff
532 63799639 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min + dhrunoff
533 63799639 : dh_direct = dhrunoff
534 25164052 : elseif (hbrocn < c0 .AND. hbr > thinS) then
535 25107603 : exp_arg = darcy_coeff/bphi_min*dt
536 : ! tcraig avoids underflows
537 25107603 : if (exp_arg > exp_argmax) then
538 8847544 : hbrocn_new = c0
539 : else
540 16260059 : hbrocn_new = hbrocn*exp(-exp_arg)
541 : endif
542 25107603 : dhflood = max(c0,hbrocn_new - hbrocn)*aice0
543 25107603 : hbr = max(hbrmin, h_ocn + hbrocn_new)
544 25107603 : darcy_V = -SIGN((hbrocn-hbrocn_new + dhflood)/dt*bphi_min, hbrocn)
545 25107603 : darcy_V_chl= darcy_V
546 25107603 : dhS_top = dhS_top - darcy_V*dt/bphi_min - dhflood
547 25107603 : dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min - dhflood
548 25107603 : dh_direct = -dhflood
549 : endif
550 :
551 88963691 : dh_bot_chl = dhS_bottom
552 :
553 : else ! very thin brine height
554 1782018 : hbrmin = min(thinS, hin)
555 1782018 : hbr = max(hbrmin, hbr_old+dhS_bottom-dhS_top)
556 1782018 : dhbr_hin = hbr - h_ocn
557 1782018 : if (abs(dhbr_hin) > dh_min) &
558 822107 : hbr = max(hbrmin, h_ocn + SIGN(dh_min,dhbr_hin))
559 1782018 : dhS_top = hbr_old - hbr + dhS_bottom
560 1782018 : dh_top_chl = dhS_top
561 1782018 : dh_bot_chl = dhS_bottom
562 : endif
563 :
564 90745709 : fbri = hbr/hin
565 :
566 90745709 : 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 4311746 : subroutine compute_microS (n_cat, nilyr, nblyr, &
579 : bgrid, cgrid, igrid, & ! LCOV_EXCL_LINE
580 : trcrn, hice_old, & ! LCOV_EXCL_LINE
581 : hbr_old, sss, sst, & ! LCOV_EXCL_LINE
582 : bTin, iTin, bphin, & ! LCOV_EXCL_LINE
583 : kperm, bphi_min, & ! LCOV_EXCL_LINE
584 : Rayleigh_criteria, firstice, & ! LCOV_EXCL_LINE
585 : bSin, brine_sal, & ! LCOV_EXCL_LINE
586 : brine_rho, iphin, ibrine_rho, & ! LCOV_EXCL_LINE
587 4311746 : ibrine_sal, sice_rho, sloss)
588 :
589 : integer (kind=int_kind), intent(in) :: &
590 : n_cat , & ! ice category ! LCOV_EXCL_LINE
591 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
605 : sss , & ! ocean salinity (ppt) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
617 : bphin , & ! Porosity of layers ! LCOV_EXCL_LINE
618 : brine_sal , & ! equilibrium brine salinity (ppt) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
629 : kperm , & ! average ice permeability (m^2) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
640 : ibrine_rho , & ! brine rho on interface ! LCOV_EXCL_LINE
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 : surface_S , & ! salinity of ice above hin > hbr ! LCOV_EXCL_LINE
650 1039770 : 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 : trtmp0 , & ! temporary, remapped tracers ! LCOV_EXCL_LINE
657 41743466 : trtmp ! temporary, remapped tracers
658 :
659 : real (kind=dbl_kind) :: &
660 1039770 : Tmlts ! melting temperature
661 :
662 : character(len=*),parameter :: subname='(compute_microS)'
663 :
664 : !-----------------------------------------------------------------
665 : ! Initialize
666 : !-----------------------------------------------------------------
667 :
668 4311746 : sloss = c0
669 43117460 : bTin(:) = c0
670 43117460 : bSin(:) = c0
671 :
672 150911110 : trtmp(:) = c0
673 4311746 : surface_S = min_salin
674 :
675 4311746 : 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 4311746 : if (firstice) then
694 :
695 153936 : do k = 1, nilyr
696 153936 : trcrn(nt_sice+k-1) = sss*salt_loss
697 : enddo
698 :
699 : call remap_zbgc(nilyr, &
700 : nt_sice, & ! LCOV_EXCL_LINE
701 : trcrn, trtmp, & ! LCOV_EXCL_LINE
702 : 0, nblyr, & ! LCOV_EXCL_LINE
703 : hinc_old, hinc_old, & ! LCOV_EXCL_LINE
704 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
705 19242 : bgrid(2:nblyr+1), surface_S )
706 19242 : if (icepack_warnings_aborted(subname)) return
707 :
708 153936 : do k = 1, nblyr
709 134694 : trcrn(nt_bgc_S+k-1) = max(min_salin,trtmp(nt_sice+k-1))
710 134694 : bSin(k+1) = max(min_salin,trcrn(nt_bgc_S+k-1))
711 134694 : if (trcrn(nt_bgc_S+k-1) < min_salin-puny) &
712 19242 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
713 : enddo ! k
714 :
715 19242 : bSin(1) = bSin(2)
716 19242 : bSin(nblyr+2) = sss
717 :
718 4292504 : elseif (hbr_old > maxhbr*hice_old) then
719 :
720 : call remap_zbgc(nblyr, &
721 : nt_bgc_S, & ! LCOV_EXCL_LINE
722 : trcrn, trtmp, & ! LCOV_EXCL_LINE
723 : 0, nblyr, & ! LCOV_EXCL_LINE
724 : hbr_old, & ! LCOV_EXCL_LINE
725 : maxhbr*hinc_old, & ! LCOV_EXCL_LINE
726 : bgrid(2:nblyr+1), & ! LCOV_EXCL_LINE
727 80 : bgrid(2:nblyr+1), surface_S )
728 80 : if (icepack_warnings_aborted(subname)) return
729 :
730 640 : do k = 1, nblyr
731 560 : bSin(k+1) = max(min_salin,trtmp(nt_bgc_S+k-1))
732 175 : sloss = sloss + rhosi*(hbr_old*trcrn(nt_bgc_S+k-1) &
733 735 : - maxhbr*hice_old*bSin(k+1))*(igrid(k+1)-igrid(k))
734 560 : trcrn(nt_bgc_S+k-1) = bSin(k+1)
735 560 : if (trcrn(nt_bgc_S+k-1) < min_salin-puny) &
736 80 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
737 : enddo ! k
738 :
739 80 : bSin(1) = bSin(2)
740 80 : bSin(nblyr+2) = sss
741 80 : hbr_old = maxhbr*hinc_old
742 :
743 : else ! old, thin ice
744 :
745 34339392 : do k = 1, nblyr
746 30046968 : trcrn(nt_bgc_S+k-1) = max(min_salin,trcrn(nt_bgc_S+k-1))
747 34339392 : bSin (k+1) = trcrn(nt_bgc_S+k-1)
748 : enddo ! k
749 :
750 4292424 : bSin (1) = bSin(2)
751 4292424 : bSin (nblyr+2) = sss
752 :
753 : endif ! ice type
754 :
755 : !-----------------------------------------------------------------
756 : ! sea ice temperature for bio grid
757 : !-----------------------------------------------------------------
758 :
759 34493968 : do k = 1, nilyr
760 30182222 : Tmlts = -trcrn(nt_sice+k-1)*depressT
761 34493968 : trtmp0(nt_qice+k-1) = calculate_Tin_from_qin(trcrn(nt_qice+k-1),Tmlts)
762 : enddo ! k
763 :
764 150911110 : trtmp(:) = c0
765 :
766 : ! CICE to Bio: remap temperatures
767 : call remap_zbgc (nilyr, nt_qice, &
768 : trtmp0(1:ntrcr), trtmp, & ! LCOV_EXCL_LINE
769 : 0, nblyr, & ! LCOV_EXCL_LINE
770 : hinc_old, hbr_old, & ! LCOV_EXCL_LINE
771 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
772 4311746 : bgrid(2:nblyr+1), surface_S )
773 4311746 : if (icepack_warnings_aborted(subname)) return
774 :
775 34493968 : do k = 1, nblyr
776 30182222 : Tmlts = -bSin(k+1) * depressT
777 34493968 : bTin (k+1) = min(Tmlts,trtmp(nt_qice+k-1))
778 : enddo !k
779 :
780 4311746 : Tmlts = -min_salin* depressT
781 4311746 : bTin (1) = min(Tmlts,(bTin(2) + trcrn(nt_Tsfc))*p5)
782 4311746 : Tmlts = -bSin(nblyr+2)* depressT
783 4311746 : bTin (nblyr+2) = sst
784 :
785 : !-----------------------------------------------------------------
786 : ! Define ice multiphase structure
787 : !-----------------------------------------------------------------
788 :
789 : call prepare_hbrine (nblyr, &
790 : bSin, bTin, iTin, & ! LCOV_EXCL_LINE
791 : brine_sal, brine_rho, & ! LCOV_EXCL_LINE
792 : ibrine_sal, ibrine_rho, & ! LCOV_EXCL_LINE
793 : sice_rho, & ! LCOV_EXCL_LINE
794 : bphin, iphin, & ! LCOV_EXCL_LINE
795 : kperm, bphi_min, & ! LCOV_EXCL_LINE
796 4311746 : igrid, sss)
797 4311746 : 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 272237127 : subroutine calculate_drho (nblyr, i_grid, b_grid, &
807 181491418 : 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 : rho_a , & ! average brine density above grid point (kg/m^3) ! LCOV_EXCL_LINE
837 189491032 : 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 : rho_b , & ! brine density above grid point (kg/m^3) ! LCOV_EXCL_LINE
841 190711108 : rho_2b ! brine density above and below grid points (kg/m^3)
842 :
843 : character(len=*),parameter :: subname='(calculate_drho)'
844 :
845 476883021 : rho_a (:) = c0
846 476883021 : rho_2a(:) = c0
847 476883021 : rho_b (:) = c0
848 476883021 : rho_2b(:) = c0
849 476883021 : drho (:) = c0 ! surface is snow or atmosphere
850 :
851 476883021 : 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 476883021 : if (k == 2) then
860 1220076 : rho_a(2) = (brine_rho(2)*b_grid(2) &
861 : + (ibrine_rho(2) + brine_rho(2)) & ! LCOV_EXCL_LINE
862 90745709 : * p5*(i_grid(2)-b_grid(2)) )/i_grid(2)
863 90745709 : rho_b(2) = brine_rho(2)
864 :
865 295391603 : elseif (k > 2 .AND. k < nblyr+1) then
866 22598460 : rho_a(k) = (rho_a(k-1)*i_grid(k-1) + (ibrine_rho(k-1) + brine_rho(k)) &
867 : * p5*(b_grid(k)-i_grid(k-1)) + (ibrine_rho(k ) + brine_rho(k)) & ! LCOV_EXCL_LINE
868 193136705 : * p5*(i_grid(k)-b_grid(k)))/i_grid(k)
869 170538245 : rho_b(k) = brine_rho(k)
870 : else
871 9399996 : rho_a(nblyr+1) = (rho_a(nblyr)*i_grid(nblyr) + (ibrine_rho(nblyr) + &
872 : brine_rho(nblyr+1))*p5*(b_grid(nblyr+1)-i_grid(nblyr)) + & ! LCOV_EXCL_LINE
873 141303351 : brine_rho(nblyr+1)*(i_grid(nblyr+1)-b_grid(nblyr+1)))/i_grid(nblyr+1)
874 124853358 : rho_a(1) = brine_rho(2) !for k == 1 use grid point value
875 124853358 : rho_b(nblyr+1) = brine_rho(nblyr+1)
876 124853358 : 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 476883021 : do k = 1, nblyr+1 !i_grid values
886 386137312 : if (k == 1) then
887 1220076 : rho_2a(1) = (rho_a(1)*b_grid(2) + p5*(brine_rho(2) + ibrine_rho(2)) &
888 90745709 : * (i_grid(2)-b_grid(2)))/i_grid(2)
889 90745709 : rho_2b(1) = brine_rho(2)
890 : else
891 295391603 : mstop = 2*(k-1) + 1
892 295391603 : if (mstop < nblyr+1) then
893 102322947 : rho_2a(k) = rho_a(mstop)
894 102322947 : mstart = 2
895 102322947 : mstop = 1
896 : else
897 193068656 : mstart = nblyr+2
898 193068656 : mstop = nblyr+3
899 : endif
900 :
901 681528915 : do mm = mstart,mstop
902 681528915 : rho_2a(k) =(rho_a(nblyr+1) + rhow*(c2*i_grid(k)-c1))*p5/i_grid(k)
903 : enddo
904 295391603 : rho_2b(k) = brine_rho(k+1)
905 : endif
906 9219690 : drho(k) = max(rho_b(k) - rho_2b(k),max(c0,c2*(rho_a(k)-rho_2a(k)), &
907 476883021 : c2*(brine_rho(k)-brine_rho(k+1))/real(nblyr,kind=dbl_kind)))
908 : enddo
909 :
910 90745709 : end subroutine calculate_drho
911 :
912 : !=======================================================================
913 : !autodocument_start icepack_init_hbrine
914 : ! Initialize brine height tracer
915 :
916 720 : subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
917 720 : icgrid, swgrid, nblyr, nilyr, phi_snow)
918 :
919 : integer (kind=int_kind), intent(in) :: &
920 : nilyr, & ! number of ice layers ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
934 : icgrid , & ! interface grid for CICE (shortwave variable) ! LCOV_EXCL_LINE
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 24 : zspace ! grid spacing for CICE vertical grid
946 :
947 : character(len=*),parameter :: subname='(icepack_init_hbrine)'
948 :
949 :
950 720 : 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 4908 : bgrid(:) = c0 ! zsalinity grid points
957 720 : bgrid(nblyr+2) = c1 ! bottom value
958 4188 : igrid(:) = c0 ! bgc interface grid points
959 720 : igrid(1) = c0 ! ice top
960 720 : igrid(nblyr+1) = c1 ! ice bottom
961 :
962 720 : zspace = c1/max(c1,(real(nblyr,kind=dbl_kind)))
963 3468 : do k = 2, nblyr+1
964 3468 : bgrid(k) = zspace*(real(k,kind=dbl_kind) - c1p5)
965 : enddo
966 :
967 2748 : do k = 2, nblyr
968 2748 : 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 720 : cgrid(1) = c0 ! CICE vertical grid top point
976 720 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
977 :
978 5760 : do k = 2, nilyr+1
979 5760 : 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 720 : icgrid(1) = c0
987 720 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
988 :
989 5760 : do k = 2, nilyr+1
990 5760 : 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 720 : zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
1000 720 : swgrid(1) = min(c1/60.0_dbl_kind, zspace/c2)
1001 720 : swgrid(2) = zspace/c2 !+ swgrid(1)
1002 5040 : do k = 3, nilyr+1
1003 5040 : swgrid(k) = zspace * (real(k,kind=dbl_kind)-c1p5)
1004 : enddo
1005 :
1006 720 : end subroutine icepack_init_hbrine
1007 :
1008 : !=======================================================================
1009 : !autodocument_start icepack_init_zsalinity
1010 : ! Initialize zSalinity
1011 :
1012 34365 : subroutine icepack_init_zsalinity(nblyr,ntrcr_o, Rayleigh_criteria, &
1013 34365 : Rayleigh_real, trcrn_bgc, nt_bgc_S, ncat, sss)
1014 :
1015 : integer (kind=int_kind), intent(in) :: &
1016 : nblyr, & ! number of biolayers ! LCOV_EXCL_LINE
1017 : ntrcr_o, & ! number of non bio tracers ! LCOV_EXCL_LINE
1018 : ncat , & ! number of categories ! LCOV_EXCL_LINE
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 34365 : if (nblyr .LE. 7) then
1043 34365 : dts_b = 300.0_dbl_kind
1044 : else
1045 0 : dts_b = 50.0_dbl_kind
1046 : endif
1047 :
1048 34365 : Rayleigh_criteria = .false. ! no ice initial condition
1049 34365 : Rayleigh_real = c0
1050 206190 : do n = 1,ncat
1051 1408965 : do k = 1,nblyr
1052 1374600 : trcrn_bgc(nt_bgc_S+k-1-ntrcr_o,n) = sss*salt_loss
1053 : enddo ! k
1054 : enddo ! n
1055 :
1056 34365 : end subroutine icepack_init_zsalinity
1057 :
1058 : !=======================================================================
1059 :
1060 : end module icepack_brine
1061 :
1062 : !=======================================================================
|
