Line data Source code
1 : !=======================================================================
2 : !
3 : ! Vertical salinity (trcrn(nt_bgc_S)) is solved on the bio grid (bgrid and igrid)
4 : ! with domain defined by the dynamic brine height (trcrn(nt_fbri) * vicen/aicen).
5 : ! The CICE Bitz and Lipscomb thermodynamics is solved on the cgrid with height
6 : ! vicen/aicen.
7 : ! Gravity drainage is parameterized as nonlinear advection
8 : ! Flushing is incorporated in the boundary changes and a darcy flow.
9 : ! (see Jeffery et al., JGR, 2011).
10 : !
11 : ! authors: Nicole Jeffery, LANL
12 : ! Elizabeth C. Hunke, LANL
13 : !
14 : module icepack_zsalinity
15 :
16 : use icepack_kinds
17 : use icepack_parameters, only: c0, c1, c2, p001, p5, puny, rhow, depressT, gravit
18 : use icepack_parameters, only: rhosi, min_salin, salt_loss
19 : use icepack_parameters, only: l_skS, grid_oS, l_sk
20 : use icepack_parameters, only: dts_b
21 : use icepack_tracers, only: nt_sice
22 : use icepack_zbgc_shared, only: remap_zbgc
23 : use icepack_zbgc_shared, only: Ra_c, k_o, viscos_dynamic, thinS, Dm, exp_h
24 : use icepack_warnings, only: warnstr, icepack_warnings_add
25 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
26 :
27 : use icepack_brine, only: calculate_drho
28 : use icepack_therm_shared, only: calculate_Tin_from_qin
29 :
30 : implicit none
31 :
32 : private
33 : public :: zsalinity
34 :
35 : real (kind=dbl_kind), parameter :: &
36 : max_salin = 200.0_dbl_kind, & !(ppt) maximum bulk salinity ! LCOV_EXCL_LINE
37 : lapidus_g = 0.3_dbl_kind ! constant for artificial
38 : ! viscosity/diffusion during growth
39 : ! lapidus_m = 0.007_dbl_kind ! constant for artificial diffusion during melt
40 :
41 : !=======================================================================
42 :
43 : contains
44 :
45 : !=======================================================================
46 :
47 4311746 : subroutine zsalinity (n_cat, dt, &
48 : nilyr, bgrid, & ! LCOV_EXCL_LINE
49 : cgrid, igrid, & ! LCOV_EXCL_LINE
50 : trcrn_S, trcrn_q, & ! LCOV_EXCL_LINE
51 : trcrn_Si, ntrcr, & ! LCOV_EXCL_LINE
52 : fbri, & ! LCOV_EXCL_LINE
53 : bSin, bTin, & ! LCOV_EXCL_LINE
54 : bphin, iphin, & ! LCOV_EXCL_LINE
55 : ikin, hbr_old, & ! LCOV_EXCL_LINE
56 : hbrin, hin, & ! LCOV_EXCL_LINE
57 : hin_old, iDin, & ! LCOV_EXCL_LINE
58 : darcy_V, brine_sal, & ! LCOV_EXCL_LINE
59 : brine_rho, ibrine_sal, & ! LCOV_EXCL_LINE
60 : ibrine_rho, dh_direct, & ! LCOV_EXCL_LINE
61 : Rayleigh_criteria, & ! LCOV_EXCL_LINE
62 : first_ice, sss, & ! LCOV_EXCL_LINE
63 : sst, dh_top, & ! LCOV_EXCL_LINE
64 : dh_bot, & ! LCOV_EXCL_LINE
65 : fzsal, & ! LCOV_EXCL_LINE
66 : fzsal_g, bphi_min, & ! LCOV_EXCL_LINE
67 : nblyr, vicen, & ! LCOV_EXCL_LINE
68 : aicen, zsal_tot)
69 :
70 : integer (kind=int_kind), intent(in) :: &
71 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
72 : nblyr , & ! number of bio layers ! LCOV_EXCL_LINE
73 : ntrcr , & ! number of tracers ! LCOV_EXCL_LINE
74 : n_cat ! category number
75 :
76 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
77 : bgrid ! biology nondimensional vertical grid points
78 :
79 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
80 : igrid ! biology vertical interface points
81 :
82 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
83 : cgrid ! CICE vertical coordinate
84 :
85 : real (kind=dbl_kind), intent(in) :: &
86 : sss , & ! ocean salinity (ppt) ! LCOV_EXCL_LINE
87 : sst , & ! ocean temperature (oC) ! LCOV_EXCL_LINE
88 : hin_old , & ! old ice thickness (m) ! LCOV_EXCL_LINE
89 : dh_top , & ! brine change in top and bottom for diagnostics (m) ! LCOV_EXCL_LINE
90 : dh_bot , & ! minimum porosity ! LCOV_EXCL_LINE
91 : darcy_V , & ! darcy velocity (m/s) ! LCOV_EXCL_LINE
92 : dt , & ! time step ! LCOV_EXCL_LINE
93 : fbri , & ! ratio of brine height to ice thickness ! LCOV_EXCL_LINE
94 : hbr_old , & ! old brine height (m) ! LCOV_EXCL_LINE
95 : hin , & ! new ice thickness (m) ! LCOV_EXCL_LINE
96 : hbrin , & ! new brine height (m) ! LCOV_EXCL_LINE
97 : vicen , & ! ice volume (m) ! LCOV_EXCL_LINE
98 : aicen , & ! ice area (m) ! LCOV_EXCL_LINE
99 : bphi_min , & ! ! LCOV_EXCL_LINE
100 : dh_direct ! flooded or runoff amount (m)
101 :
102 : real (kind=dbl_kind), intent(inout) :: &
103 : zsal_tot , & ! tot salinity (psu*rhosi*total vol ice) ! LCOV_EXCL_LINE
104 : fzsal , & ! total flux of salt out of ice over timestep(kg/m^2/s) ! LCOV_EXCL_LINE
105 : fzsal_g ! gravity drainage flux of salt over timestep(kg/m^2/s)
106 :
107 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
108 : bTin , & ! Ice Temperature ^oC (on bio grid) ! LCOV_EXCL_LINE
109 : bphin ! Ice porosity (on bio grid)
110 :
111 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
112 : bSin , & ! Ice salinity ppt (on bio grid) ! LCOV_EXCL_LINE
113 : brine_sal , & ! brine salinity (ppt) ! LCOV_EXCL_LINE
114 : brine_rho ! brine density (kg/m^3)
115 :
116 : real (kind=dbl_kind), dimension (nblyr), &
117 : intent(inout) :: & ! LCOV_EXCL_LINE
118 : trcrn_S ! salinity tracer ppt (on bio grid)
119 :
120 : real (kind=dbl_kind), dimension (nilyr), &
121 : intent(inout) :: & ! LCOV_EXCL_LINE
122 : trcrn_q , & ! enthalpy tracer ! LCOV_EXCL_LINE
123 : trcrn_Si ! salinity on CICE grid
124 :
125 : logical (kind=log_kind), intent(inout) :: &
126 : Rayleigh_criteria ! .true. if minimun ice thickness (Ra_c) was reached
127 :
128 : logical (kind=log_kind), intent(in) :: &
129 : first_ice ! for first category ice only .true.
130 : !initialized values should be used
131 :
132 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
133 : iDin , & ! Diffusivity on the igrid (1/s) ! LCOV_EXCL_LINE
134 : ikin ! permeability on the igrid
135 :
136 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
137 : iphin , & ! porosity on the igrid ! LCOV_EXCL_LINE
138 : ibrine_rho , & ! brine rho on interface ! LCOV_EXCL_LINE
139 : ibrine_sal ! brine sal on interface
140 :
141 : ! local variables
142 :
143 : real (kind=dbl_kind) :: &
144 : fzsaln , & ! category flux of salt out of ice over timestep(kg/m^2/s) ! LCOV_EXCL_LINE
145 : fzsaln_g , & ! category gravity drainage flux of salt over timestep(kg/m^2/s) ! LCOV_EXCL_LINE
146 1039770 : zsal_totn ! total salt content
147 :
148 : character(len=*),parameter :: subname='(zsalinity)'
149 :
150 : call solve_zsalinity (nilyr, nblyr, n_cat, dt, &
151 : bgrid, cgrid, igrid, & ! LCOV_EXCL_LINE
152 : trcrn_S, trcrn_q, & ! LCOV_EXCL_LINE
153 : trcrn_Si, ntrcr, & ! LCOV_EXCL_LINE
154 : bSin, bTin, & ! LCOV_EXCL_LINE
155 : bphin, iphin, & ! LCOV_EXCL_LINE
156 : ikin, hbr_old, & ! LCOV_EXCL_LINE
157 : hbrin, hin, & ! LCOV_EXCL_LINE
158 : hin_old, iDin, & ! LCOV_EXCL_LINE
159 : darcy_V, brine_sal, & ! LCOV_EXCL_LINE
160 : brine_rho, ibrine_sal, & ! LCOV_EXCL_LINE
161 : ibrine_rho, dh_direct, & ! LCOV_EXCL_LINE
162 : Rayleigh_criteria, & ! LCOV_EXCL_LINE
163 : first_ice, sss, & ! LCOV_EXCL_LINE
164 : sst, dh_top, & ! LCOV_EXCL_LINE
165 : dh_bot, & ! LCOV_EXCL_LINE
166 : fzsaln, & ! LCOV_EXCL_LINE
167 4311746 : fzsaln_g, bphi_min)
168 4311746 : if (icepack_warnings_aborted(subname)) return
169 :
170 4311746 : zsal_totn = c0
171 :
172 : call column_sum_zsal (zsal_totn, nblyr, &
173 : vicen, trcrn_S, & ! LCOV_EXCL_LINE
174 4311746 : fbri)
175 4311746 : if (icepack_warnings_aborted(subname)) return
176 :
177 : call merge_zsal_fluxes (aicen, &
178 : zsal_totn, zsal_tot, & ! LCOV_EXCL_LINE
179 : fzsal, fzsaln, & ! LCOV_EXCL_LINE
180 4311746 : fzsal_g, fzsaln_g)
181 4311746 : if (icepack_warnings_aborted(subname)) return
182 :
183 : end subroutine zsalinity
184 :
185 : !=======================================================================
186 : !
187 : ! update vertical salinity
188 : !
189 4311746 : subroutine solve_zsalinity (nilyr, nblyr, &
190 : n_cat, dt, & ! LCOV_EXCL_LINE
191 : bgrid, cgrid, igrid, & ! LCOV_EXCL_LINE
192 : trcrn_S, trcrn_q, & ! LCOV_EXCL_LINE
193 : trcrn_Si, ntrcr, & ! LCOV_EXCL_LINE
194 : bSin, bTin, & ! LCOV_EXCL_LINE
195 : bphin, iphin, & ! LCOV_EXCL_LINE
196 : ikin, hbr_old, & ! LCOV_EXCL_LINE
197 : hbrin, hin, & ! LCOV_EXCL_LINE
198 : hin_old, iDin, & ! LCOV_EXCL_LINE
199 : darcy_V, brine_sal, & ! LCOV_EXCL_LINE
200 : brine_rho, ibrine_sal, & ! LCOV_EXCL_LINE
201 : ibrine_rho, dh_direct, & ! LCOV_EXCL_LINE
202 : Rayleigh_criteria, & ! LCOV_EXCL_LINE
203 : first_ice, sss, & ! LCOV_EXCL_LINE
204 : sst, dh_top, & ! LCOV_EXCL_LINE
205 : dh_bot, & ! LCOV_EXCL_LINE
206 : fzsaln, & ! LCOV_EXCL_LINE
207 : fzsaln_g, bphi_min)
208 :
209 : integer (kind=int_kind), intent(in) :: &
210 : nilyr, & ! number of ice layers ! LCOV_EXCL_LINE
211 : nblyr, & ! number of bio layers ! LCOV_EXCL_LINE
212 : ntrcr, & ! number of tracers ! LCOV_EXCL_LINE
213 : n_cat ! category number
214 :
215 : real (kind=dbl_kind), intent(in) :: &
216 : dt ! time step
217 :
218 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
219 : bgrid ! biology nondimensional vertical grid points
220 :
221 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
222 : igrid ! biology vertical interface points
223 :
224 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
225 : cgrid ! CICE vertical coordinate
226 :
227 : real (kind=dbl_kind), intent(in) :: &
228 : sss , & ! ocean salinity (ppt) ! LCOV_EXCL_LINE
229 : sst , & ! ocean temperature (oC) ! LCOV_EXCL_LINE
230 : hin_old , & ! old ice thickness (m) ! LCOV_EXCL_LINE
231 : dh_top , & ! brine change in top and bottom for diagnostics (m) ! LCOV_EXCL_LINE
232 : dh_bot , & ! LCOV_EXCL_LINE
233 : darcy_V
234 :
235 : real (kind=dbl_kind), intent(in) :: &
236 : hbr_old , & ! old brine height (m) ! LCOV_EXCL_LINE
237 : hin , & ! new ice thickness (m) ! LCOV_EXCL_LINE
238 : hbrin , & ! new brine height (m) ! LCOV_EXCL_LINE
239 : bphi_min , & ! ! LCOV_EXCL_LINE
240 : dh_direct ! flooded or runoff amount (m)
241 :
242 : real (kind=dbl_kind), intent(out) :: &
243 : fzsaln , & ! total flux of salt out of ice over timestep(kg/m^2/s) ! LCOV_EXCL_LINE
244 : fzsaln_g ! gravity drainage flux of salt over timestep(kg/m^2/s)
245 :
246 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
247 : bTin , & ! Ice Temperature ^oC (on bio grid) ! LCOV_EXCL_LINE
248 : bphin ! Ice porosity (on bio grid)
249 :
250 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
251 : bSin , & ! Ice salinity ppt (on bio grid) ! LCOV_EXCL_LINE
252 : brine_sal , & ! brine salinity (ppt) ! LCOV_EXCL_LINE
253 : brine_rho ! brine density (kg/m^3)
254 :
255 : real (kind=dbl_kind), dimension (nblyr), &
256 : intent(inout) :: & ! LCOV_EXCL_LINE
257 : trcrn_S ! salinity tracer ppt (on bio grid)
258 :
259 : real (kind=dbl_kind), dimension (nilyr), &
260 : intent(inout) :: & ! LCOV_EXCL_LINE
261 : trcrn_q , & ! enthalpy tracer ! LCOV_EXCL_LINE
262 : trcrn_Si ! salinity on CICE grid
263 :
264 : logical (kind=log_kind), intent(inout) :: &
265 : Rayleigh_criteria ! .true. if minimun ice thickness (Ra_c) was reached
266 :
267 : logical (kind=log_kind), intent(in) :: &
268 : first_ice ! for first category ice only .true.
269 : !initialized values should be used
270 :
271 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
272 : iDin , & ! Diffusivity on the igrid (1/s) ! LCOV_EXCL_LINE
273 : ikin ! permeability on the igrid
274 :
275 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
276 : iphin , & ! porosity on the igrid ! LCOV_EXCL_LINE
277 : ibrine_rho , & ! brine rho on interface ! LCOV_EXCL_LINE
278 : ibrine_sal ! brine sal on interface
279 :
280 : ! local variables
281 :
282 : integer (kind=int_kind) :: &
283 : k, nint ! vertical biology layer index
284 :
285 : real (kind=dbl_kind) :: &
286 1039770 : surface_S ! salinity of ice above hin > hbrin
287 :
288 : real (kind=dbl_kind), dimension(2) :: &
289 3119310 : S_bot
290 :
291 : real (kind=dbl_kind) :: &
292 : Tmlts , & ! melting temperature ! LCOV_EXCL_LINE
293 1039770 : dts ! local timestep (s)
294 :
295 : logical (kind=log_kind) :: &
296 : Rayleigh
297 :
298 : real (kind=dbl_kind):: &
299 1039770 : Ttemp ! initial temp profile on the CICE grid
300 :
301 : real (kind=dbl_kind), dimension (ntrcr+2) :: &
302 : trtmp0 , & ! temporary, remapped tracers !need extra ! LCOV_EXCL_LINE
303 41743466 : trtmp ! temporary, remapped tracers !
304 :
305 : logical (kind=log_kind) :: &
306 : cflag
307 :
308 : character(len=*),parameter :: subname='(solve_zsalinity)'
309 :
310 : !-----------------------------------------------------------------
311 : ! Initialize
312 : !-----------------------------------------------------------------
313 :
314 4311746 : dts = dts_b
315 4311746 : nint = max(1,INT(dt/dts))
316 4311746 : dts = dt/nint
317 :
318 : !----------------------------------------------------------------
319 : ! Update boundary conditions
320 : !----------------------------------------------------------------
321 :
322 4311746 : surface_S = min_salin
323 :
324 4311746 : Rayleigh = .true.
325 4311746 : if (n_cat == 1 .AND. hbr_old < Ra_c) then
326 59446 : Rayleigh = Rayleigh_criteria ! only category 1 ice can be false
327 : endif
328 :
329 4311746 : if (dh_bot + darcy_V*dt > c0) then
330 :
331 3182640 : bSin (nblyr+2) = sss
332 3182640 : bTin (nblyr+2) = sst
333 3182640 : brine_sal(nblyr+2) = sss
334 3182640 : brine_rho(nblyr+2) = rhow
335 3182640 : bphin (nblyr+2) = c1
336 3182640 : S_bot (1) = c0
337 3182640 : S_bot (2) = c1
338 :
339 : ! bottom melt
340 : else
341 1129106 : bSin (nblyr+2) = bSin(nblyr+1)
342 1129106 : Tmlts = -bSin(nblyr+2)* depressT
343 1129106 : bTin (nblyr+2) = bTin(nblyr+1)
344 1129106 : bphin(nblyr+2) = iphin(nblyr+1)
345 1129106 : S_bot(1) = c1
346 1129106 : S_bot(2) = c0
347 : endif
348 :
349 4311746 : if (abs(dh_top) > puny .AND. abs(darcy_V) > puny) then
350 1002510 : bSin(1) = max(min_salin,-(brine_rho(2)*brine_sal(2)/rhosi &
351 : * darcy_V*dt - (dh_top + darcy_V*dt/bphi_min - dh_direct)*min_salin & ! LCOV_EXCL_LINE
352 4130244 : + max(c0,-dh_direct) * sss )/dh_top)
353 4130244 : brine_sal(1) = brine_sal(2)
354 4130244 : brine_rho(1) = brine_rho(2)
355 4130244 : bphin(1) = bphi_min
356 : else
357 181502 : bSin(1) = min_salin
358 : endif
359 :
360 : !-----------------------------------------------------------------
361 : ! Solve for S using CICE T. If solve_zsal = .true., then couple back
362 : ! to the thermodynamics
363 : !-----------------------------------------------------------------
364 :
365 : call solve_S_dt (cflag, nblyr, &
366 : nint , dts , & ! LCOV_EXCL_LINE
367 : bSin , bTin , & ! LCOV_EXCL_LINE
368 : bphin , iphin , & ! LCOV_EXCL_LINE
369 : igrid , bgrid , & ! LCOV_EXCL_LINE
370 : ikin , & ! LCOV_EXCL_LINE
371 : hbr_old , hbrin , & ! LCOV_EXCL_LINE
372 : hin_old , & ! LCOV_EXCL_LINE
373 : iDin , darcy_V , & ! LCOV_EXCL_LINE
374 : brine_sal , Rayleigh , & ! LCOV_EXCL_LINE
375 : first_ice , sss , & ! LCOV_EXCL_LINE
376 : dt , dh_top , & ! LCOV_EXCL_LINE
377 : dh_bot , brine_rho , & ! LCOV_EXCL_LINE
378 : ibrine_sal , ibrine_rho , & ! LCOV_EXCL_LINE
379 : fzsaln , fzsaln_g , & ! LCOV_EXCL_LINE
380 4311746 : S_bot )
381 4311746 : if (icepack_warnings_aborted(subname)) return
382 :
383 4311746 : if (n_cat == 1) Rayleigh_criteria = Rayleigh
384 :
385 150911110 : trtmp0(:) = c0
386 150911110 : trtmp (:) = c0
387 :
388 34493968 : do k = 1,nblyr ! back to bulk quantity
389 30182222 : trcrn_S(k) = bSin(k+1)
390 34493968 : trtmp0(nt_sice+k-1) = trcrn_S(k)
391 : enddo ! k
392 :
393 : call remap_zbgc (nilyr, &
394 : nt_sice, & ! LCOV_EXCL_LINE
395 : trtmp0(1:ntrcr), & ! LCOV_EXCL_LINE
396 : trtmp, & ! LCOV_EXCL_LINE
397 : 1, nblyr, & ! LCOV_EXCL_LINE
398 : hin, hbrin, & ! LCOV_EXCL_LINE
399 : cgrid(2:nilyr+1), & ! LCOV_EXCL_LINE
400 : bgrid(2:nblyr+1), & ! LCOV_EXCL_LINE
401 4311746 : surface_S )
402 4311746 : if (icepack_warnings_aborted(subname)) return
403 :
404 34493968 : do k = 1, nilyr
405 30182222 : Tmlts = -trcrn_Si(k)*depressT
406 : Ttemp = min(-(min_salin+puny)*depressT, &
407 30182222 : calculate_Tin_from_qin(trcrn_q(k),Tmlts))
408 7278390 : trcrn_Si(k) = min(-Ttemp/depressT, max(min_salin, &
409 30182222 : trtmp(nt_sice+k-1)))
410 34493968 : Tmlts = - trcrn_Si(k)*depressT
411 : ! if (cflag) trcrn_q(k) = calculate_qin_from_Sin(Ttemp,Tmlts)
412 : enddo ! k
413 :
414 : end subroutine solve_zsalinity
415 :
416 : !=======================================================================
417 : !
418 : ! solves salt continuity explicitly using
419 : ! Lax-Wendroff-type scheme (MacCormack)
420 : ! (Mendez-Nunez and Carroll, Monthly Weather Review, 1993)
421 : !
422 : ! authors Nicole Jeffery, LANL
423 : !
424 4311746 : subroutine solve_S_dt (cflag, nblyr, nint, &
425 : dts, bSin, bTin, & ! LCOV_EXCL_LINE
426 : bphin, iphin, igrid, & ! LCOV_EXCL_LINE
427 : bgrid, ikin, hbri_old, & ! LCOV_EXCL_LINE
428 : hbrin, hice_old, & ! LCOV_EXCL_LINE
429 : iDin, darcy_V, & ! LCOV_EXCL_LINE
430 : brine_sal, Rayleigh, & ! LCOV_EXCL_LINE
431 : first_ice, sss, & ! LCOV_EXCL_LINE
432 : dt, dht, & ! LCOV_EXCL_LINE
433 : dhb, brine_rho, & ! LCOV_EXCL_LINE
434 : ibrine_sal, ibrine_rho, & ! LCOV_EXCL_LINE
435 : fzsaln, fzsaln_g, & ! LCOV_EXCL_LINE
436 : S_bot )
437 :
438 : integer (kind=int_kind), intent(in) :: &
439 : nblyr , & ! number of bio layers ! LCOV_EXCL_LINE
440 : nint ! number of interations
441 :
442 : logical (kind=log_kind), intent(out) :: &
443 : cflag ! thin or not
444 :
445 : real (kind=dbl_kind), intent(in) :: &
446 : dt , & ! timestep (s) ! LCOV_EXCL_LINE
447 : dts , & ! local timestep (s) ! LCOV_EXCL_LINE
448 : sss , & ! sea surface salinity ! LCOV_EXCL_LINE
449 : dht , & ! change in the ice top (positive for melting) ! LCOV_EXCL_LINE
450 : dhb , & ! change in the ice bottom (positive for freezing) ! LCOV_EXCL_LINE
451 : hice_old , & ! old ice thickness (m) ! LCOV_EXCL_LINE
452 : hbri_old , & ! brine thickness (m) ! LCOV_EXCL_LINE
453 : hbrin , & ! new brine thickness (m) ! LCOV_EXCL_LINE
454 : darcy_V ! Darcy velocity due to a pressure head (m/s) or melt
455 :
456 : real (kind=dbl_kind), intent(out) :: &
457 : fzsaln , & ! salt flux +ive to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
458 : fzsaln_g ! gravity drainage salt flux +ive to ocean (kg/m^2/s)
459 :
460 : logical (kind=log_kind), intent(inout) :: &
461 : Rayleigh ! if .true. convection is allowed; if .false. not yet
462 :
463 : logical (kind=log_kind), intent(in) :: &
464 : first_ice
465 :
466 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
467 : brine_sal , & ! Internal brine salinity (ppt) ! LCOV_EXCL_LINE
468 : brine_rho , & ! Internal brine density (kg/m^3) ! LCOV_EXCL_LINE
469 : bgrid , & ! biology nondimensional grid layer points ! LCOV_EXCL_LINE
470 : bTin ! Temperature of ice layers on bio grid for history (C)
471 :
472 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
473 : bphin , & ! Porosity of layers ! LCOV_EXCL_LINE
474 : bSin ! Bulk Salinity (ppt) contains previous timestep
475 : ! and ocean ss
476 :
477 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
478 : ibrine_rho , & ! brine rho on interface ! LCOV_EXCL_LINE
479 : ibrine_sal , & ! brine sal on interface ! LCOV_EXCL_LINE
480 : igrid ! biology grid interface points
481 :
482 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
483 : iphin ! Porosity of layers on interface
484 :
485 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
486 : iDin , & ! Diffusivity on the igrid (1/s) with minimum bphi condition ! LCOV_EXCL_LINE
487 : ikin ! permeability on interface
488 :
489 : real (kind=dbl_kind), dimension (2), intent(in) :: &
490 : S_bot
491 :
492 : ! local variables
493 :
494 : integer (kind=int_kind) :: &
495 : k, m ! vertical biology layer index
496 :
497 : real (kind=dbl_kind), dimension (nblyr+1) :: &
498 : iDin_p , & ! Diffusivity on the igrid (1/s)/bphi^3 ! LCOV_EXCL_LINE
499 : dSbdx , & ! gradient of brine rho on grid ! LCOV_EXCL_LINE
500 : drho , & ! brine difference rho_a-rho_b (kg/m^3) ! LCOV_EXCL_LINE
501 : ! Ci_s , & ! ! LCOV_EXCL_LINE
502 : Ui_s , & ! interface function ! LCOV_EXCL_LINE
503 : ! Vi_s , & ! for conservation check ! LCOV_EXCL_LINE
504 15901882 : ivel
505 :
506 : real (kind=dbl_kind), dimension (nblyr+2) :: &
507 : Din_p , & ! Diffusivity on the igrid (1/s)/bphi^3 ! LCOV_EXCL_LINE
508 : Sintemp , & ! initial salinity ! LCOV_EXCL_LINE
509 : pre_sin , & ! estimate of salinity of layers ! LCOV_EXCL_LINE
510 : pre_sinb , & ! estimate of salinity of layers ! LCOV_EXCL_LINE
511 : bgrid_temp , & ! biology nondimensional grid layer points ! LCOV_EXCL_LINE
512 : ! with boundary values
513 33883304 : Q_s, C_s , & ! Functions in continuity equation
514 50824956 : V_s, U_s, F_s
515 :
516 : real (kind=dbl_kind) :: &
517 : dh , & ! (m) change in hbrine over dts ! LCOV_EXCL_LINE
518 : dbgrid , & ! ratio of grid space to spacing across boundary ! LCOV_EXCL_LINE
519 : ! i.e. 1/nilyr/(dbgrid(2)-dbgrid(1))
520 1039770 : lapidus , & ! artificial viscosity: use lapidus_g for growth
521 : Ssum_old,Ssum_new, & ! depth integrated salt before and after timestep ! LCOV_EXCL_LINE
522 : fluxcorr, & ! flux correction to prevent S < min_salin ! LCOV_EXCL_LINE
523 : Ssum_corr, & ! numerical boundary flux correction ! LCOV_EXCL_LINE
524 : fluxb , & ! bottom, top and total boundary flux (g/kg/m^2) ! LCOV_EXCL_LINE
525 : fluxg , & ! bottom, top and total gravity drainage flux (g/kg/m^2) ! LCOV_EXCL_LINE
526 : fluxm , & ! bottom, top and total molecular diffusion flux (g/kg/m^2) ! LCOV_EXCL_LINE
527 : sum_old,sum_new , & ! integrated salinity at t and t+dt ! LCOV_EXCL_LINE
528 : dh_dt, dS_dt , & ! LCOV_EXCL_LINE
529 1039770 : Ssum_tmp
530 :
531 : real (kind=dbl_kind), dimension (nblyr) :: &
532 : vel , & ! advective velocity times dt (m) ! LCOV_EXCL_LINE
533 : lapidus_diff , & ! lapidus term and ! LCOV_EXCL_LINE
534 : flux_corr , & ! LCOV_EXCL_LINE
535 : lapA , & ! LCOV_EXCL_LINE
536 13669676 : lapB
537 :
538 : logical (kind=log_kind) :: &
539 : test_conservation ! test that salt change is balanced by fluxes
540 :
541 : character(len=*),parameter :: subname='(solve_S_dt)'
542 :
543 : !-----------------------------------------------------------------
544 : ! Initialize
545 : !-----------------------------------------------------------------
546 :
547 4311746 : cflag = .false.
548 4311746 : test_conservation = .false.
549 38805714 : iDin_p(:) = c0
550 43117460 : Din_p(:) = c0
551 34493968 : lapA(:) = c1
552 34493968 : lapB(:) = c1
553 4311746 : lapA(nblyr) = c0
554 4311746 : lapB(1) = c0
555 43117460 : V_s(:) = c0
556 43117460 : U_s(:) = c0
557 43117460 : Q_s(:) = c0
558 43117460 : C_s(:) = c0
559 : ! Ci_s(:) = c0
560 43117460 : F_s(:) = c0
561 38805714 : Ui_s(:) = c0
562 : ! Vi_s(:) = c0
563 38805714 : ivel(:) = c0
564 34493968 : vel(:) = c0
565 4311746 : dh = c0
566 4311746 : dbgrid = c2
567 4311746 : fzsaln = c0
568 4311746 : fzsaln_g = c0
569 :
570 : !-----------------------------------------------------------------
571 : ! Find brine density gradient for gravity drainage parameterization
572 : !-----------------------------------------------------------------
573 :
574 : call calculate_drho(nblyr, igrid, bgrid,&
575 4311746 : brine_rho, ibrine_rho, drho)
576 4311746 : if (icepack_warnings_aborted(subname)) return
577 :
578 : !-----------------------------------------------------------------
579 : ! Calculate bphi diffusivity on the grid points
580 : ! rhosi = 919-974 kg/m^2 set in bio_in
581 : ! rhow = 1026.0 density of sea water: uses kinematic viscosity (m^2/s) in Q18
582 : ! dynamic viscosity divided by density = kinematic.
583 : !-----------------------------------------------------------------
584 :
585 34493968 : do k = 2, nblyr+1
586 30182222 : iDin_p(k) = k_o*gravit*l_skS/viscos_dynamic*drho(k)/(hbri_old**2)
587 7278390 : Din_p(k) = (iDin_p(k)*(igrid(k)-bgrid(k)) &
588 34493968 : + iDin_p(k-1)*(bgrid(k)-igrid(k-1)))/(igrid(k)-igrid(k-1))
589 : enddo !k
590 :
591 : !-----------------------------------------------------------------
592 : ! Critical Ra_c value is only for the onset of convection in thinS
593 : ! ice and not throughout, therefore I need a flag to indicate the
594 : ! Ra_c was reached for a particular ice block.
595 : ! Using a thickness minimum (Ra_c) for simplicity.
596 : !-----------------------------------------------------------------
597 :
598 43117460 : bgrid_temp(:) = bgrid(:)
599 4311746 : Din_p(nblyr+2) = iDin_p(nblyr+1)
600 4311746 : if (.NOT. Rayleigh .AND. hbrin < Ra_c) then
601 51280 : Din_p(:) = c0
602 46152 : iDin_p(:) = c0
603 : else
604 4306618 : Rayleigh = .true.
605 : endif
606 :
607 : if (hbri_old > thinS .AND. hbrin > thinS .and. &
608 4311746 : hice_old > thinS .AND. .NOT. first_ice) then
609 :
610 4227698 : cflag = .true.
611 :
612 4227698 : bgrid_temp(1) = c0
613 4227698 : bgrid_temp(nblyr+2) = c1
614 4227698 : dbgrid = igrid(2)/(bgrid_temp(2)-bgrid_temp(1))
615 :
616 : !-----------------------------------
617 : ! surface boundary terms
618 : !-----------------------------------
619 :
620 4227698 : lapidus = lapidus_g/real(nblyr,kind=dbl_kind)**2
621 4227698 : ivel(1) = dht/hbri_old
622 4227698 : U_s (1) = ivel(1)/dt*dts
623 4227698 : Ui_s(1) = U_s(1)
624 : ! Ci_s(1) = c0
625 4227698 : F_s (1) = brine_rho(2)*brine_sal(2)/rhosi*darcy_V*dts/hbri_old/bSin(1)
626 :
627 : !-----------------------------------
628 : ! bottom boundary terms
629 : !-----------------------------------
630 :
631 4227698 : ivel(nblyr+1) = dhb/hbri_old
632 4227698 : Ui_s(nblyr+1) = ivel(nblyr+1)/dt*dts
633 3046935 : dSbdx(nblyr) = (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) &
634 : - ibrine_sal(nblyr)*ibrine_rho(nblyr)) & ! LCOV_EXCL_LINE
635 9305923 : / (igrid(nblyr+1)-igrid(nblyr))
636 3046935 : C_s(nblyr+1) = Dm/brine_sal(nblyr+1)/brine_rho(nblyr+1)*dts/hbri_old**2 &
637 : * (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) & ! LCOV_EXCL_LINE
638 : - ibrine_sal(nblyr)*ibrine_rho(nblyr)) & ! LCOV_EXCL_LINE
639 11337213 : / (igrid(nblyr+1)-igrid(nblyr))
640 4227698 : F_s(nblyr+1) = darcy_V*dts/hbri_old/bphin(nblyr+1)
641 4227698 : F_s(nblyr+2) = darcy_V*dts/hbri_old/bphin(nblyr+2)
642 4227698 : vel(nblyr) =(bgrid(nblyr+1)*(dhb) -(bgrid(nblyr+1) - c1)*dht)/hbri_old
643 4227698 : U_s(nblyr+1) = vel(nblyr)/dt*dts
644 2031290 : V_s(nblyr+1) = Din_p(nblyr+1)/rhosi &
645 : * (rhosi/brine_sal(nblyr+1)/brine_rho(nblyr+1))**exp_h & ! LCOV_EXCL_LINE
646 8290278 : * dts*dSbdx(nblyr)
647 3046935 : dSbdx(nblyr+1) = (brine_sal(nblyr+2)*brine_rho(nblyr+2) &
648 : - brine_sal(nblyr+1)*brine_rho(nblyr+1)) & ! LCOV_EXCL_LINE
649 9305923 : / (bgrid(nblyr+2)-bgrid(nblyr+1)+ grid_oS/hbri_old )
650 3046935 : C_s( nblyr+2) = Dm/brine_sal(nblyr+2)/brine_rho(nblyr+2)*dts/hbri_old**2 &
651 : * (brine_sal(nblyr+2)*brine_rho(nblyr+2) & ! LCOV_EXCL_LINE
652 : - brine_sal(nblyr+1)*brine_rho(nblyr+1)) & ! LCOV_EXCL_LINE
653 11337213 : / (bgrid(nblyr+2)-bgrid(nblyr+1) + grid_oS/hbri_old )
654 4227698 : U_s(nblyr+2) = ivel(nblyr+1)/dt*dts
655 2031290 : V_s(nblyr+2) = Din_p(nblyr+2)/rhosi &
656 : * (bphin(nblyr+1)/bSin(nblyr+2))**exp_h & ! LCOV_EXCL_LINE
657 8290278 : * dts*dSbdx(nblyr+1)
658 : ! Ci_s(nblyr+1) = C_s(nblyr+2)
659 : ! Vi_s(nblyr+1) = V_s(nblyr+2)
660 4227698 : dh = (dhb-dht)/dt*dts
661 :
662 29593886 : do k = 2, nblyr
663 25366188 : ivel(k) = (igrid(k)*dhb - (igrid(k)-c1)*dht)/hbri_old
664 25366188 : Ui_s(k) = ivel(k)/dt*dts
665 : ! Vi_s(k) = iDin_p(k)/rhosi &
666 : ! *(rhosi/ibrine_rho(k)/ibrine_sal(k))**exp_h*dts & ! LCOV_EXCL_LINE
667 : ! * (brine_sal(k+1)*brine_rho(k+1) & ! LCOV_EXCL_LINE
668 : ! - brine_sal(k)*brine_rho(k)) & ! LCOV_EXCL_LINE
669 : ! / (bgrid(k+1)-bgrid(k))
670 6093870 : dSbdx(k-1) = (ibrine_sal(k)*ibrine_rho(k) &
671 31460058 : - ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1))
672 25366188 : F_s(k) = darcy_V*dts/hbri_old/bphin(k)
673 6093870 : C_s(k) = Dm/brine_sal(k)/brine_rho(k)*dts/hbri_old**2 &
674 : * (ibrine_sal(k)*ibrine_rho(k) & ! LCOV_EXCL_LINE
675 31460058 : - ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1))
676 : ! Ci_s(k) = Dm/ibrine_sal(k)/ibrine_rho(k)*dts/hbri_old**2 &
677 : ! * (brine_sal(k+1)*brine_rho(k+1) & ! LCOV_EXCL_LINE
678 : ! - brine_sal(k)*brine_rho(k))/(bgrid(k+1)-bgrid(k))
679 25366188 : vel(k-1) = (bgrid(k)*(dhb) - (bgrid(k) - c1)* dht)/hbri_old
680 25366188 : U_s(k) = vel(k-1)/dt*dts
681 6093870 : V_s(k) = Din_p(k)/rhosi &
682 25366188 : * (rhosi/brine_sal(k)/brine_rho(k))**exp_h*dts*dSbdx(k-1)
683 25366188 : C_s(2) = c0
684 29593886 : V_s(2) = c0
685 : enddo !k
686 :
687 : !-----------------------------------------------------------------
688 : ! Solve
689 : !-----------------------------------------------------------------
690 :
691 59187772 : do m = 1, nint
692 :
693 507323760 : Sintemp(:) = bSin(:)
694 507323760 : pre_sin(:) = bSin(:)
695 507323760 : pre_sinb(:) = bSin(:)
696 50732376 : Ssum_old = bSin(nblyr+1)*(igrid(nblyr+1)-igrid(nblyr))
697 :
698 : ! forward-difference
699 :
700 355126632 : do k = 2, nblyr
701 304394256 : Ssum_old = Ssum_old + bSin(k)*(igrid(k)-igrid(k-1))
702 :
703 146252880 : pre_sin(k) =bSin(k) + (Ui_s(k)*(bSin(k+1) - bSin(k)) + &
704 : V_s(k+1)*bSin(k+1)**3 - V_s(k)*bSin(k)**3 + & ! LCOV_EXCL_LINE
705 : (C_s(k+1)+F_s(k+1))*bSin(k+1)-& ! LCOV_EXCL_LINE
706 867011712 : (C_s(k)+F_s(k))*bSin(k))/(bgrid_temp(k+1)-bgrid_temp(k))
707 : enddo !k
708 :
709 48750960 : pre_sin(nblyr+1) = bSin(nblyr+1) + (Ui_s(nblyr+1)*(bSin(nblyr+2) - &
710 : bSin(nblyr+1)) + V_s(nblyr+2)*bSin(nblyr+2)**3 - & ! LCOV_EXCL_LINE
711 : V_s(nblyr+1)*bSin(nblyr+1)**3+ (C_s(nblyr+2)+F_s(nblyr+2))*& ! LCOV_EXCL_LINE
712 : bSin(nblyr+2)-(C_s(nblyr+1)+F_s(nblyr+1))*bSin(nblyr+1) )/& ! LCOV_EXCL_LINE
713 233548476 : (bgrid_temp(nblyr+2)- bgrid_temp(nblyr+1))
714 :
715 : ! backward-difference
716 :
717 12187740 : pre_sinb(2) = p5*(bSin(2) + pre_sin(2) + (Ui_s(1)&
718 : *(pre_sin(2) -pre_sin(1)) + & ! LCOV_EXCL_LINE
719 : V_s(2)*pre_sin(2)**3 - & ! LCOV_EXCL_LINE
720 : V_s(1)*pre_sin(1)**3 + (C_s(2)+F_s(2))*pre_sin(2)-& ! LCOV_EXCL_LINE
721 75107856 : (C_s(1)+F_s(1))*pre_sin(1) )/(bgrid_temp(2)-bgrid_temp(1)) )
722 :
723 355126632 : do k = nblyr+1, 3, -1 !nblyr+1
724 146252880 : pre_sinb(k) =p5*(bSin(k) + pre_sin(k) + (Ui_s(k-1)&
725 : *(pre_sin(k) - pre_sin(k-1)) + & ! LCOV_EXCL_LINE
726 : V_s(k)*pre_sin(k)**3 - & ! LCOV_EXCL_LINE
727 : V_s(k-1)*pre_sin(k-1)**3 + (C_s(k)+F_s(k))*pre_sin(k)-& ! LCOV_EXCL_LINE
728 670026456 : (C_s(k-1)+F_s(k-1))*pre_sin(k-1))/(bgrid_temp(k)-bgrid_temp(k-1)) )
729 :
730 355126632 : Q_s(k) = V_s(k)*pre_sin(k)**2 + U_s(k) + C_s(k) + F_s(k)
731 : enddo !k
732 :
733 50732376 : Q_s(2) = V_s(2)*pre_sin(2)**2 + U_s(2) + C_s(2) + F_s(2)
734 50732376 : Q_s(1) = V_s(1)*pre_sin(2)**2 + Ui_s(1) + C_s(1)+ F_s(1)
735 36563220 : Q_s(nblyr+2) = V_s(nblyr+2)*pre_sin(nblyr+1)**2 + &
736 87295596 : Ui_s(nblyr+1) + C_s(nblyr+2) + F_s(nblyr+2)
737 :
738 : !-----------------------------------------------------------------
739 : ! Add artificial viscosity [Lapidus,1967] [Lohner et al, 1985]
740 : ! * more for melting ice
741 : !-----------------------------------------------------------------
742 :
743 405859008 : lapidus_diff(:) = c0
744 405859008 : flux_corr(:) = c0
745 50732376 : Ssum_new = c0
746 50732376 : Ssum_corr = c0
747 50732376 : fluxcorr = c0
748 50732376 : fluxg = c0
749 50732376 : fluxb = c0
750 50732376 : fluxm = c0
751 :
752 405859008 : do k = 2, nblyr+1
753 :
754 85314180 : lapidus_diff(k-1) = lapidus/& ! lapidus/real(nblyr,kind=dbl_kind)**2/&
755 : (igrid(k)-igrid(k-1))* & ! LCOV_EXCL_LINE
756 : ( lapA(k-1)*ABS(Q_s(k+1)-Q_s(k))*(abs(pre_sinb(k+1))-abs(pre_sinb(k)))/& ! LCOV_EXCL_LINE
757 : (bgrid_temp(k+1)-bgrid_temp(k) )**2 - & ! LCOV_EXCL_LINE
758 : lapB(k-1)*ABS(Q_s(k)-Q_s(k-1))*(abs(pre_sinb(k))-abs(pre_sinb(k-1)))/& ! LCOV_EXCL_LINE
759 1122954252 : (bgrid_temp(k)-bgrid_temp(k-1))**2)
760 :
761 355126632 : bSin(k) = pre_sinb(k) + lapidus_diff(k-1)
762 :
763 355126632 : if (bSin(k) < min_salin) then
764 122 : flux_corr(k-1) = min_salin - bSin(k) ! flux into the ice
765 122 : bSin(k) = min_salin
766 355126510 : elseif (bSin(k) > -bTin(k)/depressT) then
767 21678 : flux_corr(k-1) = bSin(k)+bTin(k)/depressT ! flux into the ice
768 21678 : bSin(k) = -bTin(k)/depressT
769 355104832 : elseif (bSin(k) > max_salin) then
770 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
771 0 : call icepack_warnings_add(subname//' bSin(k) > max_salin')
772 : endif
773 :
774 379502112 : if (k == nblyr+1) bSin(nblyr+2) = S_bot(1)*bSin(nblyr+1) &
775 75107856 : + S_bot(2)*bSin(nblyr+2)
776 :
777 355126632 : Ssum_new = Ssum_new + bSin(k)*(igrid(k)-igrid(k-1))
778 85314180 : fluxcorr = fluxcorr + (flux_corr(k-1) &
779 491173188 : + lapidus_diff(k-1))*(igrid(k)-igrid(k-1))
780 :
781 : enddo !k
782 :
783 50732376 : Ssum_tmp = Ssum_old
784 :
785 : call calc_salt_fluxes (nint, nblyr, &
786 : Ui_s, dh,dbgrid,hbri_old,Sintemp, & ! LCOV_EXCL_LINE
787 : pre_sin, fluxb,fluxg,fluxm,V_s, & ! LCOV_EXCL_LINE
788 : C_s, F_s, Ssum_corr,fzsaln_g,fzsaln, & ! LCOV_EXCL_LINE
789 50732376 : Ssum_tmp,dts, Ssum_new)
790 50732376 : if (icepack_warnings_aborted(subname)) return
791 :
792 54960074 : if (test_conservation) then
793 : call check_conserve_salt(nint, m, dt, &
794 : Ssum_tmp, Ssum_new, Ssum_corr,& ! LCOV_EXCL_LINE
795 : fluxcorr, fluxb, fluxg, fluxm, & ! LCOV_EXCL_LINE
796 0 : hbrin, hbri_old)
797 0 : call icepack_warnings_add(subname//' check_conserve_salt fails')
798 0 : if (icepack_warnings_aborted(subname)) return
799 : endif ! test_conservation
800 :
801 : enddo !m
802 :
803 : else ! add/melt ice only
804 :
805 84048 : sum_old = c0
806 84048 : sum_new = c0
807 84048 : dh_dt = hbrin-hbri_old
808 84048 : dS_dt = c0
809 84048 : if (dh_dt > c0) then
810 2684 : dS_dt = sss*dh_dt*salt_loss
811 21472 : do k = 2, nblyr+1
812 21472 : bSin(k) = max(min_salin,(bSin(k)*hbri_old + dS_dt)/hbrin)
813 : enddo !k
814 : else
815 650912 : do k = 2, nblyr+1
816 569548 : sum_old = sum_old + bSin(k)*hbri_old*(igrid(k)-igrid(k-1))
817 569548 : bSin(k) = max(min_salin,bSin(k) * (c1-abs(dh_dt)/hbri_old))
818 650912 : sum_new = sum_new + bSin(k)*hbrin*(igrid(k)-igrid(k-1))
819 : enddo !k
820 : endif
821 84048 : fzsaln = rhosi*(sum_old-sum_new + dS_dt)*p001/dt !kg/m^2/s
822 84048 : fzsaln_g = c0
823 :
824 : endif ! (hbri_old > thinS .AND. hbrin > thinS &
825 : ! .and. hice_old > thinS .AND. .NOT. first_ice)
826 :
827 : !-----------------------------------------------------------------
828 : ! Move this to bgc calculation if using tr_salinity
829 : ! Calculate bphin, iphin, ikin, iDin and iDin_N
830 : !-----------------------------------------------------------------
831 :
832 38805714 : iDin(:) = c0
833 38805714 : iphin(:) = c1
834 38805714 : ikin(:) = c0
835 :
836 38805714 : do k = 1, nblyr+1
837 41772358 : if (k < nblyr+1) bphin(k+1) = min(c1,max(puny, &
838 37460612 : bSin(k+1)*rhosi/(brine_sal(k+1)*brine_rho(k+1))))
839 34493968 : if (k == 1) then
840 4311746 : bphin(k) = min(c1,max(puny, bSin(k)*rhosi/(brine_sal(k)*brine_rho(k))))
841 4311746 : iphin(k) = bphin(2)
842 30182222 : elseif (k == nblyr+1) then
843 4311746 : iphin(nblyr+1) = bphin(nblyr+1)
844 : else
845 18715860 : iphin(k) = min(c1, max(c0,(bphin(k+1) - bphin(k))/(bgrid(k+1) &
846 38347716 : - bgrid(k))*(igrid(k)-bgrid(k)) + bphin(k)))
847 : endif
848 38805714 : ikin(k) = k_o*iphin(k)**exp_h
849 : enddo !k
850 :
851 4311746 : if (cflag) then
852 :
853 33821584 : do k = 2, nblyr+1
854 29593886 : iDin(k) = iphin(k)*Dm/hbri_old**2
855 29593886 : if (Rayleigh .AND. hbrin .GE. Ra_c) &
856 : iDin(k) = iDin(k) + l_sk*ikin(k)*gravit/viscos_dynamic & ! LCOV_EXCL_LINE
857 33821584 : * drho(k)/hbri_old**2
858 : enddo !k
859 : else ! .not. cflag
860 672384 : do k = 2, nblyr+1
861 672384 : iDin(k) = iphin(k)*Dm/hbri_old**2
862 : enddo !k
863 : endif
864 :
865 : end subroutine solve_S_dt
866 :
867 : !=======================================================================
868 : !
869 : ! Calculate salt fluxes
870 : !
871 50732376 : subroutine calc_salt_fluxes (mint, nblyr, &
872 : Ui_s,dh,dbgrid,hbri_old,Sintemp,pre_sin,& ! LCOV_EXCL_LINE
873 : fluxb,fluxg,fluxm,V_s,& ! LCOV_EXCL_LINE
874 : C_s,F_s,Ssum_corr,fzsaln_g,fzsaln,Ssum_old, & ! LCOV_EXCL_LINE
875 : ! fluxcorr,dts, Ssum_new)
876 : dts, Ssum_new)
877 :
878 : integer(kind=int_kind), intent(in) :: &
879 : nblyr, & ! number of bio layers ! LCOV_EXCL_LINE
880 : mint ! current iteration
881 :
882 : real (kind=dbl_kind), intent(in) :: &
883 : dts , & ! halodynamic timesteps (s) ! LCOV_EXCL_LINE
884 : ! hbrin , & ! new brine height after all iterations (m) ! LCOV_EXCL_LINE
885 : dh , & ! (m) change in hbrine over dts ! LCOV_EXCL_LINE
886 : dbgrid , & ! ratio of grid space to spacing across boundary ! LCOV_EXCL_LINE
887 : ! ie. 1/nilyr/(dbgrid(2)-dbgrid(1))
888 : ! fluxcorr , & ! flux correction to ensure S >= min_salin
889 : hbri_old ! initial brine height (m)
890 :
891 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
892 : Ui_s ! interface function
893 :
894 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
895 : Sintemp , & ! initial salinity ! LCOV_EXCL_LINE
896 : pre_sin , & ! estimate of salinity of layers ! LCOV_EXCL_LINE
897 : C_s , & ! Functions in continuity equation ! LCOV_EXCL_LINE
898 : F_s , & ! LCOV_EXCL_LINE
899 : V_s
900 :
901 : real (kind=dbl_kind), intent(in) :: &
902 : Ssum_old , & ! initial integrated salt content (ppt)/h ! LCOV_EXCL_LINE
903 : Ssum_new ! next integrated salt content(ppt)/h
904 :
905 : real (kind=dbl_kind), intent(inout) :: &
906 : fzsaln , & ! total salt flux out of ice over timestep(kg/m^2/s) ! LCOV_EXCL_LINE
907 : fzsaln_g , & ! gravity drainage flux of salt over timestep(kg/m^2/s) ! LCOV_EXCL_LINE
908 : Ssum_corr, & ! boundary flux correction due to numerics ! LCOV_EXCL_LINE
909 : fluxb , & ! total boundary salt flux into the ice (+ into ice) ! LCOV_EXCL_LINE
910 : fluxg , & ! total gravity drainage salt flux into the ice (+ into ice) ! LCOV_EXCL_LINE
911 : fluxm ! total molecular diffusive salt flux into the ice (+ into ice)
912 :
913 : ! local variables
914 :
915 : real (kind=dbl_kind) :: &
916 : Ssum_corr_flux , & ! numerical boundary flux correction ! LCOV_EXCL_LINE
917 : fluxb_b, fluxb_t, & ! bottom, top and total boundary flux (g/kg/m^2) ! LCOV_EXCL_LINE
918 : fluxg_b, fluxg_t, & ! bottom, top and total gravity drainage flux (g/kg/m^2) ! LCOV_EXCL_LINE
919 12187740 : fluxm_b, fluxm_t ! bottom, top and total molecular diffusion flux (g/kg/m^2)
920 :
921 24375480 : real (kind=dbl_kind) :: hin_old, hin_next, dhtmp !, dh
922 :
923 : character(len=*),parameter :: subname='(calc_salt_fluxes)'
924 :
925 50732376 : dhtmp = c1-dh/hbri_old
926 50732376 : hin_next = hbri_old + real(mint,kind=dbl_kind)*dh
927 50732376 : hin_old = hbri_old + (real(mint,kind=dbl_kind)-c1)*dh
928 :
929 : !-----------------------------------------------------------------
930 : ! boundary fluxes (positive into the ice)
931 : !---------------------------------------------
932 : ! without higher order numerics corrections
933 : ! fluxb = (Ui_s(nblyr+1) + F_s(nblyr+2))*Sintemp(nblyr+2) - (Ui_s(1) + F_s(1))*Sintemp(1)
934 : !-----------------------------------------------------------------
935 :
936 24375480 : fluxb_b = p5* Ui_s(nblyr+1) * (dhtmp*Sintemp(nblyr+2)*dbgrid &
937 : + pre_sin(nblyr+1) & ! LCOV_EXCL_LINE
938 : + dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)) & ! LCOV_EXCL_LINE
939 : + p5*(F_s(nblyr+2) * dhtmp*Sintemp(nblyr+2)*dbgrid & ! LCOV_EXCL_LINE
940 : + F_s(nblyr+1) * (pre_sin(nblyr+1) & ! LCOV_EXCL_LINE
941 148234296 : + dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)))
942 :
943 12187740 : fluxb_t = -p5*Ui_s(1)*(pre_sin(1)*dbgrid + &
944 : dhtmp*Sintemp(2) - & ! LCOV_EXCL_LINE
945 : (dbgrid-c1)*pre_sin(2)) - & ! LCOV_EXCL_LINE
946 : p5*(dbgrid*F_s(1)*pre_sin(1) + & ! LCOV_EXCL_LINE
947 : F_s(2)*(dhtmp*Sintemp(2) & ! LCOV_EXCL_LINE
948 50732376 : +(c1-dbgrid)*pre_sin(2)))
949 :
950 50732376 : fluxb = fluxb_b + fluxb_t
951 :
952 : !-----------------------------------------------------------------
953 : ! gravity drainage fluxes (positive into the ice)
954 : ! without higher order numerics corrections
955 : ! fluxg = V_s(nblyr+2)*Sintemp(nblyr+1)**3
956 : !-----------------------------------------------------------------
957 :
958 : fluxg_b = p5*(dhtmp* dbgrid* &
959 : V_s(nblyr+2)*Sintemp(nblyr+1)**3 + & ! LCOV_EXCL_LINE
960 : V_s(nblyr+1)*(pre_sin(nblyr+1)**3 - & ! LCOV_EXCL_LINE
961 : dhtmp*(dbgrid - c1)* & ! LCOV_EXCL_LINE
962 99483336 : Sintemp(nblyr+1)**3))
963 :
964 12187740 : fluxg_t = -p5*(dbgrid*V_s(1)*pre_sin(1)**3 + &
965 : V_s(2)*(dhtmp*Sintemp(2)**3- & ! LCOV_EXCL_LINE
966 62920116 : (dbgrid-c1)*pre_sin(2)**3))
967 :
968 50732376 : fluxg = fluxg_b + fluxg_t
969 :
970 : !-----------------------------------------------------------------
971 : ! diffusion fluxes (positive into the ice)
972 : ! without higher order numerics corrections
973 : ! fluxm = C_s(nblyr+2)*Sintemp(nblyr+2)
974 : !-----------------------------------------------------------------
975 :
976 24375480 : fluxm_b = p5*(dhtmp*C_s(nblyr+2)* Sintemp(nblyr+2)*dbgrid &
977 : + C_s(nblyr+1)*(pre_sin(nblyr+1) & ! LCOV_EXCL_LINE
978 99483336 : + dhtmp * Sintemp(nblyr+1)*(c1-dbgrid)))
979 :
980 12187740 : fluxm_t = -p5 * (C_s(1) * pre_sin(1)*dbgrid &
981 50732376 : + C_s(2) * (pre_sin(2)*(c1-dbgrid) + dhtmp*Sintemp(2)))
982 :
983 50732376 : fluxm = fluxm_b + fluxm_t
984 :
985 50732376 : Ssum_corr = (-dh/hbri_old + p5*(dh/hbri_old)**2)*Ssum_old
986 50732376 : Ssum_corr_flux = dh*Ssum_old/hin_next + Ssum_corr
987 50732376 : Ssum_corr = Ssum_corr_flux
988 :
989 : fzsaln_g = fzsaln_g - hin_next * fluxg_b &
990 50732376 : * rhosi*p001/dts
991 :
992 : !approximate fluxes
993 : !fzsaln = fzsaln - hin_next * (fluxg &
994 : ! + fluxb + fluxm + fluxcorr + Ssum_corr_flux) & ! LCOV_EXCL_LINE
995 : ! * rhosi*p001/dts
996 :
997 : fzsaln = fzsaln + (Ssum_old*hin_old - Ssum_new*hin_next) &
998 50732376 : * rhosi*p001/dts ! positive into the ocean
999 :
1000 50732376 : end subroutine calc_salt_fluxes
1001 :
1002 : !=======================================================================
1003 : !
1004 : ! Test salt conservation: flux conservative form d(hSin)/dt = -dF(x,Sin)/dx
1005 : !
1006 0 : subroutine check_conserve_salt (mmax, mint, dt, &
1007 : Ssum_old, Ssum_new, Ssum_corr, & ! LCOV_EXCL_LINE
1008 : fluxcorr, fluxb, fluxg, fluxm, & ! LCOV_EXCL_LINE
1009 : hbrin, hbri_old)
1010 :
1011 : integer(kind=int_kind), intent(in) :: &
1012 : mint , & ! current iteration ! LCOV_EXCL_LINE
1013 : mmax ! maximum number of iterations
1014 :
1015 : real (kind=dbl_kind), intent(in) :: &
1016 : dt , & ! thermodynamic and halodynamic timesteps (s) ! LCOV_EXCL_LINE
1017 : hbrin , & ! (m) final brine height ! LCOV_EXCL_LINE
1018 : hbri_old , & ! (m) initial brine height ! LCOV_EXCL_LINE
1019 : Ssum_old , & ! initial integrated salt content ! LCOV_EXCL_LINE
1020 : Ssum_new , & ! final integrated salt content ! LCOV_EXCL_LINE
1021 : fluxcorr , & ! flux correction to ensure S >= min_salin ! LCOV_EXCL_LINE
1022 : Ssum_corr , & ! boundary flux correction due to numerics ! LCOV_EXCL_LINE
1023 : fluxb , & ! total boundary salt flux into the ice (+ into ice) ! LCOV_EXCL_LINE
1024 : fluxg , & ! total gravity drainage salt flux into the ice (+ into ice) ! LCOV_EXCL_LINE
1025 : fluxm !
1026 :
1027 : ! local variables
1028 :
1029 : real (kind=dbl_kind):: &
1030 : diff2 , & ! ! LCOV_EXCL_LINE
1031 : dsum_flux , & ! salt change in kg/m^2/s ! LCOV_EXCL_LINE
1032 : flux_tot , & ! fluxg + fluxb ! LCOV_EXCL_LINE
1033 : order , & ! ! LCOV_EXCL_LINE
1034 0 : dh
1035 :
1036 : real (kind=dbl_kind), parameter :: &
1037 : accuracy = 1.0e-7_dbl_kind ! g/kg/m^2/s difference between boundary fluxes
1038 :
1039 : character(len=*),parameter :: subname='(check_conserve_salt)'
1040 :
1041 0 : dh = (hbrin-hbri_old)/real(mmax,kind=dbl_kind)
1042 :
1043 : flux_tot = (fluxb + fluxg + fluxm + fluxcorr + Ssum_corr)*&
1044 0 : (hbri_old + (real(mint,kind=dbl_kind))*dh)/dt
1045 : dsum_flux =(Ssum_new*(hbri_old + (real(mint,kind=dbl_kind))*dh) - &
1046 : Ssum_old*(hbri_old + (real(mint,kind=dbl_kind)-c1)* & ! LCOV_EXCL_LINE
1047 0 : dh) )/dt
1048 0 : order = abs(dh/min(hbri_old,hbrin))
1049 0 : if (abs(dsum_flux) > accuracy) then
1050 0 : diff2 = abs(dsum_flux - flux_tot)
1051 0 : if (diff2 > puny .AND. diff2 > order ) then
1052 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1053 0 : write(warnstr,*) subname, 'Poor salt conservation: check_conserve_salt'
1054 0 : call icepack_warnings_add(warnstr)
1055 0 : write(warnstr,*) subname, 'mint:', mint
1056 0 : call icepack_warnings_add(warnstr)
1057 0 : write(warnstr,*) subname, 'Ssum_corr',Ssum_corr
1058 0 : call icepack_warnings_add(warnstr)
1059 0 : write(warnstr,*) subname, 'fluxb,fluxg,fluxm,flux_tot,fluxcorr:'
1060 0 : call icepack_warnings_add(warnstr)
1061 0 : write(warnstr,*) subname, fluxb,fluxg,fluxm,flux_tot,fluxcorr
1062 0 : call icepack_warnings_add(warnstr)
1063 0 : write(warnstr,*) subname, 'fluxg,',fluxg
1064 0 : call icepack_warnings_add(warnstr)
1065 0 : write(warnstr,*) subname, 'dsum_flux,',dsum_flux
1066 0 : call icepack_warnings_add(warnstr)
1067 0 : write(warnstr,*) subname, 'Ssum_new,Ssum_old,hbri_old,dh:'
1068 0 : call icepack_warnings_add(warnstr)
1069 0 : write(warnstr,*) subname, Ssum_new,Ssum_old,hbri_old,dh
1070 0 : call icepack_warnings_add(warnstr)
1071 0 : write(warnstr,*) subname, 'diff2,order,puny',diff2,order,puny
1072 0 : call icepack_warnings_add(warnstr)
1073 : endif
1074 : endif
1075 :
1076 0 : end subroutine check_conserve_salt
1077 :
1078 : !=======================================================================
1079 : !
1080 : ! Aggregate flux information from all ice thickness categories
1081 : !
1082 4311746 : subroutine merge_zsal_fluxes(aicenS, &
1083 : zsal_totn, zsal_tot, & ! LCOV_EXCL_LINE
1084 : fzsal, fzsaln, & ! LCOV_EXCL_LINE
1085 : fzsal_g, fzsaln_g)
1086 :
1087 : ! single category fluxes
1088 : real (kind=dbl_kind), intent(in):: &
1089 : aicenS , & ! concentration of ice ! LCOV_EXCL_LINE
1090 : fzsaln , & ! salt flux (kg/m**2/s) ! LCOV_EXCL_LINE
1091 : fzsaln_g ! gravity drainage salt flux (kg/m**2/s)
1092 :
1093 : real (kind=dbl_kind), intent(in):: &
1094 : zsal_totn ! tot salinity in category (psu*volume*rhosi)
1095 :
1096 : real (kind=dbl_kind), intent(inout):: &
1097 : zsal_tot, & ! tot salinity (psu*rhosi*total vol ice) ! LCOV_EXCL_LINE
1098 : fzsal , & ! salt flux (kg/m**2/s) ! LCOV_EXCL_LINE
1099 : fzsal_g ! gravity drainage salt flux (kg/m**2/s)
1100 :
1101 : character(len=*),parameter :: subname='(merge_zsal_fluxes)'
1102 :
1103 : !-----------------------------------------------------------------
1104 : ! Merge fluxes
1105 : !-----------------------------------------------------------------
1106 :
1107 4311746 : zsal_tot = zsal_tot + zsal_totn ! already *aicenS
1108 :
1109 : ! ocean tot and gravity drainage salt fluxes
1110 4311746 : fzsal = fzsal + fzsaln * aicenS
1111 4311746 : fzsal_g = fzsal_g + fzsaln_g * aicenS
1112 :
1113 4311746 : end subroutine merge_zsal_fluxes
1114 :
1115 : !==========================================================================
1116 : !
1117 : ! For each grid cell, sum field over all ice layers. "Net" refers to the column
1118 : ! integration while "avg" is normalized by the ice thickness
1119 :
1120 4311746 : subroutine column_sum_zsal (zsal_totn, nblyr, &
1121 4311746 : vicenS, trcrn_S, fbri)
1122 :
1123 : integer (kind=int_kind), intent(in) :: &
1124 : nblyr ! number of layers
1125 :
1126 : real (kind=dbl_kind), intent(in):: &
1127 : vicenS , & ! volume of ice (m) ! LCOV_EXCL_LINE
1128 : fbri ! brine height to ice thickness ratio
1129 :
1130 : real (kind=dbl_kind), dimension (nblyr), intent(in) :: &
1131 : trcrn_S ! input field
1132 :
1133 : real (kind=dbl_kind), intent(inout) :: &
1134 : zsal_totn ! avg salinity (psu*rhosi*vol of ice)
1135 :
1136 : ! local variables
1137 :
1138 : integer (kind=int_kind) :: &
1139 : k ! layer index
1140 :
1141 : character(len=*),parameter :: subname='(column_sum_zsal)'
1142 :
1143 34493968 : do k = 1, nblyr
1144 : zsal_totn = zsal_totn &
1145 : + rhosi * trcrn_S(k) & ! LCOV_EXCL_LINE
1146 : * fbri & ! LCOV_EXCL_LINE
1147 34493968 : * vicenS/real(nblyr,kind=dbl_kind)
1148 : enddo ! k
1149 :
1150 4311746 : end subroutine column_sum_zsal
1151 :
1152 : !=======================================================================
1153 :
1154 : end module icepack_zsalinity
1155 :
1156 : !=======================================================================
|
