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
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 0 : subroutine zsalinity (n_cat, dt, &
48 0 : nilyr, bgrid, &
49 0 : cgrid, igrid, &
50 0 : trcrn_S, trcrn_q, &
51 0 : trcrn_Si, ntrcr, &
52 : fbri, &
53 0 : bSin, bTin, &
54 0 : bphin, iphin, &
55 0 : ikin, hbr_old, &
56 : hbrin, hin, &
57 0 : hin_old, iDin, &
58 0 : darcy_V, brine_sal, &
59 0 : brine_rho, ibrine_sal, &
60 0 : ibrine_rho, dh_direct, &
61 : Rayleigh_criteria, &
62 : first_ice, sss, &
63 : sst, dh_top, &
64 : dh_bot, &
65 : fzsal, &
66 : fzsal_g, bphi_min, &
67 : nblyr, vicen, &
68 : aicen, zsal_tot)
69 :
70 : integer (kind=int_kind), intent(in) :: &
71 : nilyr , & ! number of ice layers
72 : nblyr , & ! number of bio layers
73 : ntrcr , & ! number of tracers
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)
87 : sst , & ! ocean temperature (oC)
88 : hin_old , & ! old ice thickness (m)
89 : dh_top , & ! brine change in top and bottom for diagnostics (m)
90 : dh_bot , & ! minimum porosity
91 : darcy_V , & ! darcy velocity (m/s)
92 : dt , & ! time step
93 : fbri , & ! ratio of brine height to ice thickness
94 : hbr_old , & ! old brine height (m)
95 : hin , & ! new ice thickness (m)
96 : hbrin , & ! new brine height (m)
97 : vicen , & ! ice volume (m)
98 : aicen , & ! ice area (m)
99 : bphi_min , & !
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)
104 : fzsal , & ! total flux of salt out of ice over timestep(kg/m^2/s)
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)
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)
113 : brine_sal , & ! brine salinity (ppt)
114 : brine_rho ! brine density (kg/m^3)
115 :
116 : real (kind=dbl_kind), dimension (nblyr), &
117 : intent(inout) :: &
118 : trcrn_S ! salinity tracer ppt (on bio grid)
119 :
120 : real (kind=dbl_kind), dimension (nilyr), &
121 : intent(inout) :: &
122 : trcrn_q , & ! enthalpy tracer
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)
134 : ikin ! permeability on the igrid
135 :
136 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
137 : iphin , & ! porosity on the igrid
138 : ibrine_rho , & ! brine rho on interface
139 : ibrine_sal ! brine sal on interface
140 :
141 : ! local variables
142 :
143 : real (kind=dbl_kind) :: &
144 0 : fzsaln , & ! category flux of salt out of ice over timestep(kg/m^2/s)
145 0 : fzsaln_g , & ! category gravity drainage flux of salt over timestep(kg/m^2/s)
146 0 : 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, &
152 : trcrn_S, trcrn_q, &
153 : trcrn_Si, ntrcr, &
154 : bSin, bTin, &
155 : bphin, iphin, &
156 : ikin, hbr_old, &
157 : hbrin, hin, &
158 : hin_old, iDin, &
159 : darcy_V, brine_sal, &
160 : brine_rho, ibrine_sal, &
161 : ibrine_rho, dh_direct, &
162 : Rayleigh_criteria, &
163 : first_ice, sss, &
164 : sst, dh_top, &
165 : dh_bot, &
166 : fzsaln, &
167 0 : fzsaln_g, bphi_min)
168 0 : if (icepack_warnings_aborted(subname)) return
169 :
170 0 : zsal_totn = c0
171 :
172 : call column_sum_zsal (zsal_totn, nblyr, &
173 : vicen, trcrn_S, &
174 0 : fbri)
175 0 : if (icepack_warnings_aborted(subname)) return
176 :
177 : call merge_zsal_fluxes (aicen, &
178 : zsal_totn, zsal_tot, &
179 : fzsal, fzsaln, &
180 0 : fzsal_g, fzsaln_g)
181 0 : if (icepack_warnings_aborted(subname)) return
182 :
183 : end subroutine zsalinity
184 :
185 : !=======================================================================
186 : !
187 : ! update vertical salinity
188 : !
189 0 : subroutine solve_zsalinity (nilyr, nblyr, &
190 : n_cat, dt, &
191 0 : bgrid, cgrid, igrid, &
192 0 : trcrn_S, trcrn_q, &
193 0 : trcrn_Si, ntrcr, &
194 0 : bSin, bTin, &
195 0 : bphin, iphin, &
196 0 : ikin, hbr_old, &
197 : hbrin, hin, &
198 0 : hin_old, iDin, &
199 0 : darcy_V, brine_sal, &
200 0 : brine_rho, ibrine_sal, &
201 0 : ibrine_rho, dh_direct, &
202 : Rayleigh_criteria, &
203 : first_ice, sss, &
204 : sst, dh_top, &
205 : dh_bot, &
206 : fzsaln, &
207 : fzsaln_g, bphi_min)
208 :
209 : integer (kind=int_kind), intent(in) :: &
210 : nilyr, & ! number of ice layers
211 : nblyr, & ! number of bio layers
212 : ntrcr, & ! number of tracers
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)
229 : sst , & ! ocean temperature (oC)
230 : hin_old , & ! old ice thickness (m)
231 : dh_top , & ! brine change in top and bottom for diagnostics (m)
232 : dh_bot , &
233 : darcy_V
234 :
235 : real (kind=dbl_kind), intent(in) :: &
236 : hbr_old , & ! old brine height (m)
237 : hin , & ! new ice thickness (m)
238 : hbrin , & ! new brine height (m)
239 : bphi_min , & !
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)
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)
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)
252 : brine_sal , & ! brine salinity (ppt)
253 : brine_rho ! brine density (kg/m^3)
254 :
255 : real (kind=dbl_kind), dimension (nblyr), &
256 : intent(inout) :: &
257 : trcrn_S ! salinity tracer ppt (on bio grid)
258 :
259 : real (kind=dbl_kind), dimension (nilyr), &
260 : intent(inout) :: &
261 : trcrn_q , & ! enthalpy tracer
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)
273 : ikin ! permeability on the igrid
274 :
275 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
276 : iphin , & ! porosity on the igrid
277 : ibrine_rho , & ! brine rho on interface
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 0 : surface_S ! salinity of ice above hin > hbrin
287 :
288 : real (kind=dbl_kind), dimension(2) :: &
289 0 : S_bot
290 :
291 : real (kind=dbl_kind) :: &
292 0 : Tmlts , & ! melting temperature
293 0 : dts ! local timestep (s)
294 :
295 : logical (kind=log_kind) :: &
296 : Rayleigh
297 :
298 : real (kind=dbl_kind):: &
299 0 : Ttemp ! initial temp profile on the CICE grid
300 :
301 : real (kind=dbl_kind), dimension (ntrcr+2) :: &
302 0 : trtmp0 , & ! temporary, remapped tracers !need extra
303 0 : 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 0 : dts = dts_b
315 0 : nint = max(1,INT(dt/dts))
316 0 : dts = dt/nint
317 :
318 : !----------------------------------------------------------------
319 : ! Update boundary conditions
320 : !----------------------------------------------------------------
321 :
322 0 : surface_S = min_salin
323 :
324 0 : Rayleigh = .true.
325 0 : if (n_cat == 1 .AND. hbr_old < Ra_c) then
326 0 : Rayleigh = Rayleigh_criteria ! only category 1 ice can be false
327 : endif
328 :
329 0 : if (dh_bot + darcy_V*dt > c0) then
330 :
331 0 : bSin (nblyr+2) = sss
332 0 : bTin (nblyr+2) = sst
333 0 : brine_sal(nblyr+2) = sss
334 0 : brine_rho(nblyr+2) = rhow
335 0 : bphin (nblyr+2) = c1
336 0 : S_bot (1) = c0
337 0 : S_bot (2) = c1
338 :
339 : ! bottom melt
340 : else
341 0 : bSin (nblyr+2) = bSin(nblyr+1)
342 0 : Tmlts = -bSin(nblyr+2)* depressT
343 0 : bTin (nblyr+2) = bTin(nblyr+1)
344 0 : bphin(nblyr+2) = iphin(nblyr+1)
345 0 : S_bot(1) = c1
346 0 : S_bot(2) = c0
347 : endif
348 :
349 0 : if (abs(dh_top) > puny .AND. abs(darcy_V) > puny) then
350 0 : 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 &
352 0 : + max(c0,-dh_direct) * sss )/dh_top)
353 0 : brine_sal(1) = brine_sal(2)
354 0 : brine_rho(1) = brine_rho(2)
355 0 : bphin(1) = bphi_min
356 : else
357 0 : 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 , &
367 : bSin , bTin , &
368 : bphin , iphin , &
369 : igrid , bgrid , &
370 : ikin , &
371 : hbr_old , hbrin , &
372 : hin_old , &
373 : iDin , darcy_V , &
374 : brine_sal , Rayleigh , &
375 : first_ice , sss , &
376 : dt , dh_top , &
377 : dh_bot , brine_rho , &
378 : ibrine_sal , ibrine_rho , &
379 : fzsaln , fzsaln_g , &
380 0 : S_bot )
381 0 : if (icepack_warnings_aborted(subname)) return
382 :
383 0 : if (n_cat == 1) Rayleigh_criteria = Rayleigh
384 :
385 0 : trtmp0(:) = c0
386 0 : trtmp (:) = c0
387 :
388 0 : do k = 1,nblyr ! back to bulk quantity
389 0 : trcrn_S(k) = bSin(k+1)
390 0 : trtmp0(nt_sice+k-1) = trcrn_S(k)
391 : enddo ! k
392 :
393 : call remap_zbgc (nilyr, &
394 : nt_sice, &
395 0 : trtmp0(1:ntrcr), &
396 0 : trtmp, &
397 : 1, nblyr, &
398 : hin, hbrin, &
399 0 : cgrid(2:nilyr+1), &
400 0 : bgrid(2:nblyr+1), &
401 0 : surface_S )
402 0 : if (icepack_warnings_aborted(subname)) return
403 :
404 0 : do k = 1, nilyr
405 0 : Tmlts = -trcrn_Si(k)*depressT
406 : Ttemp = min(-(min_salin+puny)*depressT, &
407 0 : calculate_Tin_from_qin(trcrn_q(k),Tmlts))
408 0 : trcrn_Si(k) = min(-Ttemp/depressT, max(min_salin, &
409 0 : trtmp(nt_sice+k-1)))
410 0 : 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 0 : subroutine solve_S_dt (cflag, nblyr, nint, &
425 0 : dts, bSin, bTin, &
426 0 : bphin, iphin, igrid, &
427 0 : bgrid, ikin, hbri_old, &
428 : hbrin, hice_old, &
429 0 : iDin, darcy_V, &
430 0 : brine_sal, Rayleigh, &
431 : first_ice, sss, &
432 : dt, dht, &
433 0 : dhb, brine_rho, &
434 0 : ibrine_sal, ibrine_rho, &
435 : fzsaln, fzsaln_g, &
436 : S_bot )
437 :
438 : integer (kind=int_kind), intent(in) :: &
439 : nblyr , & ! number of bio layers
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)
447 : dts , & ! local timestep (s)
448 : sss , & ! sea surface salinity
449 : dht , & ! change in the ice top (positive for melting)
450 : dhb , & ! change in the ice bottom (positive for freezing)
451 : hice_old , & ! old ice thickness (m)
452 : hbri_old , & ! brine thickness (m)
453 : hbrin , & ! new brine thickness (m)
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)
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)
468 : brine_rho , & ! Internal brine density (kg/m^3)
469 : bgrid , & ! biology nondimensional grid layer points
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
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
479 : ibrine_sal , & ! brine sal on interface
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
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 0 : iDin_p , & ! Diffusivity on the igrid (1/s)/bphi^3
499 0 : dSbdx , & ! gradient of brine rho on grid
500 0 : drho , & ! brine difference rho_a-rho_b (kg/m^3)
501 : ! Ci_s , & !
502 0 : Ui_s , & ! interface function
503 : ! Vi_s , & ! for conservation check
504 0 : ivel
505 :
506 : real (kind=dbl_kind), dimension (nblyr+2) :: &
507 0 : Din_p , & ! Diffusivity on the igrid (1/s)/bphi^3
508 0 : Sintemp , & ! initial salinity
509 0 : pre_sin , & ! estimate of salinity of layers
510 0 : pre_sinb , & ! estimate of salinity of layers
511 0 : bgrid_temp , & ! biology nondimensional grid layer points
512 : ! with boundary values
513 0 : Q_s, C_s , & ! Functions in continuity equation
514 0 : V_s, U_s, F_s
515 :
516 : real (kind=dbl_kind) :: &
517 0 : dh , & ! (m) change in hbrine over dts
518 0 : dbgrid , & ! ratio of grid space to spacing across boundary
519 : ! i.e. 1/nilyr/(dbgrid(2)-dbgrid(1))
520 0 : lapidus , & ! artificial viscosity: use lapidus_g for growth
521 0 : Ssum_old,Ssum_new, & ! depth integrated salt before and after timestep
522 0 : fluxcorr, & ! flux correction to prevent S < min_salin
523 0 : Ssum_corr, & ! numerical boundary flux correction
524 0 : fluxb , & ! bottom, top and total boundary flux (g/kg/m^2)
525 0 : fluxg , & ! bottom, top and total gravity drainage flux (g/kg/m^2)
526 0 : fluxm , & ! bottom, top and total molecular diffusion flux (g/kg/m^2)
527 0 : sum_old,sum_new , & ! integrated salinity at t and t+dt
528 0 : dh_dt, dS_dt , &
529 0 : Ssum_tmp
530 :
531 : real (kind=dbl_kind), dimension (nblyr) :: &
532 0 : vel , & ! advective velocity times dt (m)
533 0 : lapidus_diff , & ! lapidus term and
534 0 : flux_corr , &
535 0 : lapA , &
536 0 : 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 0 : cflag = .false.
548 0 : test_conservation = .false.
549 0 : iDin_p(:) = c0
550 0 : Din_p(:) = c0
551 0 : lapA(:) = c1
552 0 : lapB(:) = c1
553 0 : lapA(nblyr) = c0
554 0 : lapB(1) = c0
555 0 : V_s(:) = c0
556 0 : U_s(:) = c0
557 0 : Q_s(:) = c0
558 0 : C_s(:) = c0
559 : ! Ci_s(:) = c0
560 0 : F_s(:) = c0
561 0 : Ui_s(:) = c0
562 : ! Vi_s(:) = c0
563 0 : ivel(:) = c0
564 0 : vel(:) = c0
565 0 : dh = c0
566 0 : dbgrid = c2
567 0 : fzsaln = c0
568 0 : fzsaln_g = c0
569 :
570 : !-----------------------------------------------------------------
571 : ! Find brine density gradient for gravity drainage parameterization
572 : !-----------------------------------------------------------------
573 :
574 : call calculate_drho(nblyr, igrid, bgrid,&
575 0 : brine_rho, ibrine_rho, drho)
576 0 : 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 0 : do k = 2, nblyr+1
586 0 : iDin_p(k) = k_o*gravit*l_skS/viscos_dynamic*drho(k)/(hbri_old**2)
587 0 : Din_p(k) = (iDin_p(k)*(igrid(k)-bgrid(k)) &
588 0 : + 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 0 : bgrid_temp(:) = bgrid(:)
599 0 : Din_p(nblyr+2) = iDin_p(nblyr+1)
600 0 : if (.NOT. Rayleigh .AND. hbrin < Ra_c) then
601 0 : Din_p(:) = c0
602 0 : iDin_p(:) = c0
603 : else
604 0 : Rayleigh = .true.
605 : endif
606 :
607 : if (hbri_old > thinS .AND. hbrin > thinS .and. &
608 0 : hice_old > thinS .AND. .NOT. first_ice) then
609 :
610 0 : cflag = .true.
611 :
612 0 : bgrid_temp(1) = c0
613 0 : bgrid_temp(nblyr+2) = c1
614 0 : dbgrid = igrid(2)/(bgrid_temp(2)-bgrid_temp(1))
615 :
616 : !-----------------------------------
617 : ! surface boundary terms
618 : !-----------------------------------
619 :
620 0 : lapidus = lapidus_g/real(nblyr,kind=dbl_kind)**2
621 0 : ivel(1) = dht/hbri_old
622 0 : U_s (1) = ivel(1)/dt*dts
623 0 : Ui_s(1) = U_s(1)
624 : ! Ci_s(1) = c0
625 0 : 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 0 : ivel(nblyr+1) = dhb/hbri_old
632 0 : Ui_s(nblyr+1) = ivel(nblyr+1)/dt*dts
633 0 : dSbdx(nblyr) = (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) &
634 0 : - ibrine_sal(nblyr)*ibrine_rho(nblyr)) &
635 0 : / (igrid(nblyr+1)-igrid(nblyr))
636 0 : C_s(nblyr+1) = Dm/brine_sal(nblyr+1)/brine_rho(nblyr+1)*dts/hbri_old**2 &
637 0 : * (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) &
638 0 : - ibrine_sal(nblyr)*ibrine_rho(nblyr)) &
639 0 : / (igrid(nblyr+1)-igrid(nblyr))
640 0 : F_s(nblyr+1) = darcy_V*dts/hbri_old/bphin(nblyr+1)
641 0 : F_s(nblyr+2) = darcy_V*dts/hbri_old/bphin(nblyr+2)
642 0 : vel(nblyr) =(bgrid(nblyr+1)*(dhb) -(bgrid(nblyr+1) - c1)*dht)/hbri_old
643 0 : U_s(nblyr+1) = vel(nblyr)/dt*dts
644 0 : V_s(nblyr+1) = Din_p(nblyr+1)/rhosi &
645 0 : * (rhosi/brine_sal(nblyr+1)/brine_rho(nblyr+1))**exp_h &
646 0 : * dts*dSbdx(nblyr)
647 0 : dSbdx(nblyr+1) = (brine_sal(nblyr+2)*brine_rho(nblyr+2) &
648 0 : - brine_sal(nblyr+1)*brine_rho(nblyr+1)) &
649 0 : / (bgrid(nblyr+2)-bgrid(nblyr+1)+ grid_oS/hbri_old )
650 0 : C_s( nblyr+2) = Dm/brine_sal(nblyr+2)/brine_rho(nblyr+2)*dts/hbri_old**2 &
651 0 : * (brine_sal(nblyr+2)*brine_rho(nblyr+2) &
652 0 : - brine_sal(nblyr+1)*brine_rho(nblyr+1)) &
653 0 : / (bgrid(nblyr+2)-bgrid(nblyr+1) + grid_oS/hbri_old )
654 0 : U_s(nblyr+2) = ivel(nblyr+1)/dt*dts
655 0 : V_s(nblyr+2) = Din_p(nblyr+2)/rhosi &
656 0 : * (bphin(nblyr+1)/bSin(nblyr+2))**exp_h &
657 0 : * dts*dSbdx(nblyr+1)
658 : ! Ci_s(nblyr+1) = C_s(nblyr+2)
659 : ! Vi_s(nblyr+1) = V_s(nblyr+2)
660 0 : dh = (dhb-dht)/dt*dts
661 :
662 0 : do k = 2, nblyr
663 0 : ivel(k) = (igrid(k)*dhb - (igrid(k)-c1)*dht)/hbri_old
664 0 : 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 &
667 : ! * (brine_sal(k+1)*brine_rho(k+1) &
668 : ! - brine_sal(k)*brine_rho(k)) &
669 : ! / (bgrid(k+1)-bgrid(k))
670 0 : dSbdx(k-1) = (ibrine_sal(k)*ibrine_rho(k) &
671 0 : - ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1))
672 0 : F_s(k) = darcy_V*dts/hbri_old/bphin(k)
673 0 : C_s(k) = Dm/brine_sal(k)/brine_rho(k)*dts/hbri_old**2 &
674 0 : * (ibrine_sal(k)*ibrine_rho(k) &
675 0 : - 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) &
678 : ! - brine_sal(k)*brine_rho(k))/(bgrid(k+1)-bgrid(k))
679 0 : vel(k-1) = (bgrid(k)*(dhb) - (bgrid(k) - c1)* dht)/hbri_old
680 0 : U_s(k) = vel(k-1)/dt*dts
681 0 : V_s(k) = Din_p(k)/rhosi &
682 0 : * (rhosi/brine_sal(k)/brine_rho(k))**exp_h*dts*dSbdx(k-1)
683 0 : C_s(2) = c0
684 0 : V_s(2) = c0
685 : enddo !k
686 :
687 : !-----------------------------------------------------------------
688 : ! Solve
689 : !-----------------------------------------------------------------
690 :
691 0 : do m = 1, nint
692 :
693 0 : Sintemp(:) = bSin(:)
694 0 : pre_sin(:) = bSin(:)
695 0 : pre_sinb(:) = bSin(:)
696 0 : Ssum_old = bSin(nblyr+1)*(igrid(nblyr+1)-igrid(nblyr))
697 :
698 : ! forward-difference
699 :
700 0 : do k = 2, nblyr
701 0 : Ssum_old = Ssum_old + bSin(k)*(igrid(k)-igrid(k-1))
702 :
703 0 : pre_sin(k) =bSin(k) + (Ui_s(k)*(bSin(k+1) - bSin(k)) + &
704 0 : V_s(k+1)*bSin(k+1)**3 - V_s(k)*bSin(k)**3 + &
705 0 : (C_s(k+1)+F_s(k+1))*bSin(k+1)-&
706 0 : (C_s(k)+F_s(k))*bSin(k))/(bgrid_temp(k+1)-bgrid_temp(k))
707 : enddo !k
708 :
709 0 : pre_sin(nblyr+1) = bSin(nblyr+1) + (Ui_s(nblyr+1)*(bSin(nblyr+2) - &
710 0 : bSin(nblyr+1)) + V_s(nblyr+2)*bSin(nblyr+2)**3 - &
711 0 : V_s(nblyr+1)*bSin(nblyr+1)**3+ (C_s(nblyr+2)+F_s(nblyr+2))*&
712 0 : bSin(nblyr+2)-(C_s(nblyr+1)+F_s(nblyr+1))*bSin(nblyr+1) )/&
713 0 : (bgrid_temp(nblyr+2)- bgrid_temp(nblyr+1))
714 :
715 : ! backward-difference
716 :
717 0 : pre_sinb(2) = p5*(bSin(2) + pre_sin(2) + (Ui_s(1)&
718 0 : *(pre_sin(2) -pre_sin(1)) + &
719 0 : V_s(2)*pre_sin(2)**3 - &
720 0 : V_s(1)*pre_sin(1)**3 + (C_s(2)+F_s(2))*pre_sin(2)-&
721 0 : (C_s(1)+F_s(1))*pre_sin(1) )/(bgrid_temp(2)-bgrid_temp(1)) )
722 :
723 0 : do k = nblyr+1, 3, -1 !nblyr+1
724 0 : pre_sinb(k) =p5*(bSin(k) + pre_sin(k) + (Ui_s(k-1)&
725 0 : *(pre_sin(k) - pre_sin(k-1)) + &
726 0 : V_s(k)*pre_sin(k)**3 - &
727 0 : V_s(k-1)*pre_sin(k-1)**3 + (C_s(k)+F_s(k))*pre_sin(k)-&
728 0 : (C_s(k-1)+F_s(k-1))*pre_sin(k-1))/(bgrid_temp(k)-bgrid_temp(k-1)) )
729 :
730 0 : Q_s(k) = V_s(k)*pre_sin(k)**2 + U_s(k) + C_s(k) + F_s(k)
731 : enddo !k
732 :
733 0 : Q_s(2) = V_s(2)*pre_sin(2)**2 + U_s(2) + C_s(2) + F_s(2)
734 0 : Q_s(1) = V_s(1)*pre_sin(2)**2 + Ui_s(1) + C_s(1)+ F_s(1)
735 0 : Q_s(nblyr+2) = V_s(nblyr+2)*pre_sin(nblyr+1)**2 + &
736 0 : 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 0 : lapidus_diff(:) = c0
744 0 : flux_corr(:) = c0
745 0 : Ssum_new = c0
746 0 : Ssum_corr = c0
747 0 : fluxcorr = c0
748 0 : fluxg = c0
749 0 : fluxb = c0
750 0 : fluxm = c0
751 :
752 0 : do k = 2, nblyr+1
753 :
754 0 : lapidus_diff(k-1) = lapidus/& ! lapidus/real(nblyr,kind=dbl_kind)**2/&
755 0 : (igrid(k)-igrid(k-1))* &
756 0 : ( lapA(k-1)*ABS(Q_s(k+1)-Q_s(k))*(abs(pre_sinb(k+1))-abs(pre_sinb(k)))/&
757 0 : (bgrid_temp(k+1)-bgrid_temp(k) )**2 - &
758 0 : lapB(k-1)*ABS(Q_s(k)-Q_s(k-1))*(abs(pre_sinb(k))-abs(pre_sinb(k-1)))/&
759 0 : (bgrid_temp(k)-bgrid_temp(k-1))**2)
760 :
761 0 : bSin(k) = pre_sinb(k) + lapidus_diff(k-1)
762 :
763 0 : if (bSin(k) < min_salin) then
764 0 : flux_corr(k-1) = min_salin - bSin(k) ! flux into the ice
765 0 : bSin(k) = min_salin
766 0 : elseif (bSin(k) > -bTin(k)/depressT) then
767 0 : flux_corr(k-1) = bSin(k)+bTin(k)/depressT ! flux into the ice
768 0 : bSin(k) = -bTin(k)/depressT
769 0 : 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 0 : if (k == nblyr+1) bSin(nblyr+2) = S_bot(1)*bSin(nblyr+1) &
775 0 : + S_bot(2)*bSin(nblyr+2)
776 :
777 0 : Ssum_new = Ssum_new + bSin(k)*(igrid(k)-igrid(k-1))
778 0 : fluxcorr = fluxcorr + (flux_corr(k-1) &
779 0 : + lapidus_diff(k-1))*(igrid(k)-igrid(k-1))
780 :
781 : enddo !k
782 :
783 0 : Ssum_tmp = Ssum_old
784 :
785 : call calc_salt_fluxes (nint, nblyr, &
786 : Ui_s, dh,dbgrid,hbri_old,Sintemp, &
787 : pre_sin, fluxb,fluxg,fluxm,V_s, &
788 : C_s, F_s, Ssum_corr,fzsaln_g,fzsaln, &
789 0 : Ssum_tmp,dts, Ssum_new)
790 0 : if (icepack_warnings_aborted(subname)) return
791 :
792 0 : if (test_conservation) then
793 : call check_conserve_salt(nint, m, dt, &
794 : Ssum_tmp, Ssum_new, Ssum_corr,&
795 : fluxcorr, fluxb, fluxg, fluxm, &
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 0 : sum_old = c0
806 0 : sum_new = c0
807 0 : dh_dt = hbrin-hbri_old
808 0 : dS_dt = c0
809 0 : if (dh_dt > c0) then
810 0 : dS_dt = sss*dh_dt*salt_loss
811 0 : do k = 2, nblyr+1
812 0 : bSin(k) = max(min_salin,(bSin(k)*hbri_old + dS_dt)/hbrin)
813 : enddo !k
814 : else
815 0 : do k = 2, nblyr+1
816 0 : sum_old = sum_old + bSin(k)*hbri_old*(igrid(k)-igrid(k-1))
817 0 : bSin(k) = max(min_salin,bSin(k) * (c1-abs(dh_dt)/hbri_old))
818 0 : sum_new = sum_new + bSin(k)*hbrin*(igrid(k)-igrid(k-1))
819 : enddo !k
820 : endif
821 0 : fzsaln = rhosi*(sum_old-sum_new + dS_dt)*p001/dt !kg/m^2/s
822 0 : 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 0 : iDin(:) = c0
833 0 : iphin(:) = c1
834 0 : ikin(:) = c0
835 :
836 0 : do k = 1, nblyr+1
837 0 : if (k < nblyr+1) bphin(k+1) = min(c1,max(puny, &
838 0 : bSin(k+1)*rhosi/(brine_sal(k+1)*brine_rho(k+1))))
839 0 : if (k == 1) then
840 0 : bphin(k) = min(c1,max(puny, bSin(k)*rhosi/(brine_sal(k)*brine_rho(k))))
841 0 : iphin(k) = bphin(2)
842 0 : elseif (k == nblyr+1) then
843 0 : iphin(nblyr+1) = bphin(nblyr+1)
844 : else
845 0 : iphin(k) = min(c1, max(c0,(bphin(k+1) - bphin(k))/(bgrid(k+1) &
846 0 : - bgrid(k))*(igrid(k)-bgrid(k)) + bphin(k)))
847 : endif
848 0 : ikin(k) = k_o*iphin(k)**exp_h
849 : enddo !k
850 :
851 0 : if (cflag) then
852 :
853 0 : do k = 2, nblyr+1
854 0 : iDin(k) = iphin(k)*Dm/hbri_old**2
855 0 : if (Rayleigh .AND. hbrin .GE. Ra_c) &
856 0 : iDin(k) = iDin(k) + l_sk*ikin(k)*gravit/viscos_dynamic &
857 0 : * drho(k)/hbri_old**2
858 : enddo !k
859 : else ! .not. cflag
860 0 : do k = 2, nblyr+1
861 0 : 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 0 : subroutine calc_salt_fluxes (mint, nblyr, &
872 0 : Ui_s,dh,dbgrid,hbri_old,Sintemp,pre_sin,&
873 0 : fluxb,fluxg,fluxm,V_s,&
874 0 : C_s,F_s,Ssum_corr,fzsaln_g,fzsaln,Ssum_old, &
875 : ! fluxcorr,dts, Ssum_new)
876 : dts, Ssum_new)
877 :
878 : integer(kind=int_kind), intent(in) :: &
879 : nblyr, & ! number of bio layers
880 : mint ! current iteration
881 :
882 : real (kind=dbl_kind), intent(in) :: &
883 : dts , & ! halodynamic timesteps (s)
884 : ! hbrin , & ! new brine height after all iterations (m)
885 : dh , & ! (m) change in hbrine over dts
886 : dbgrid , & ! ratio of grid space to spacing across boundary
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
896 : pre_sin , & ! estimate of salinity of layers
897 : C_s , & ! Functions in continuity equation
898 : F_s , &
899 : V_s
900 :
901 : real (kind=dbl_kind), intent(in) :: &
902 : Ssum_old , & ! initial integrated salt content (ppt)/h
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)
907 : fzsaln_g , & ! gravity drainage flux of salt over timestep(kg/m^2/s)
908 : Ssum_corr, & ! boundary flux correction due to numerics
909 : fluxb , & ! total boundary salt flux into the ice (+ into ice)
910 : fluxg , & ! total gravity drainage salt flux into the ice (+ into ice)
911 : fluxm ! total molecular diffusive salt flux into the ice (+ into ice)
912 :
913 : ! local variables
914 :
915 : real (kind=dbl_kind) :: &
916 0 : Ssum_corr_flux , & ! numerical boundary flux correction
917 0 : fluxb_b, fluxb_t, & ! bottom, top and total boundary flux (g/kg/m^2)
918 0 : fluxg_b, fluxg_t, & ! bottom, top and total gravity drainage flux (g/kg/m^2)
919 0 : fluxm_b, fluxm_t ! bottom, top and total molecular diffusion flux (g/kg/m^2)
920 :
921 0 : real (kind=dbl_kind) :: hin_old, hin_next, dhtmp !, dh
922 :
923 : character(len=*),parameter :: subname='(calc_salt_fluxes)'
924 :
925 0 : dhtmp = c1-dh/hbri_old
926 0 : hin_next = hbri_old + real(mint,kind=dbl_kind)*dh
927 0 : 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 0 : fluxb_b = p5* Ui_s(nblyr+1) * (dhtmp*Sintemp(nblyr+2)*dbgrid &
937 0 : + pre_sin(nblyr+1) &
938 0 : + dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)) &
939 0 : + p5*(F_s(nblyr+2) * dhtmp*Sintemp(nblyr+2)*dbgrid &
940 0 : + F_s(nblyr+1) * (pre_sin(nblyr+1) &
941 0 : + dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)))
942 :
943 0 : fluxb_t = -p5*Ui_s(1)*(pre_sin(1)*dbgrid + &
944 0 : dhtmp*Sintemp(2) - &
945 0 : (dbgrid-c1)*pre_sin(2)) - &
946 0 : p5*(dbgrid*F_s(1)*pre_sin(1) + &
947 0 : F_s(2)*(dhtmp*Sintemp(2) &
948 0 : +(c1-dbgrid)*pre_sin(2)))
949 :
950 0 : 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 0 : V_s(nblyr+2)*Sintemp(nblyr+1)**3 + &
960 0 : V_s(nblyr+1)*(pre_sin(nblyr+1)**3 - &
961 : dhtmp*(dbgrid - c1)* &
962 0 : Sintemp(nblyr+1)**3))
963 :
964 0 : fluxg_t = -p5*(dbgrid*V_s(1)*pre_sin(1)**3 + &
965 0 : V_s(2)*(dhtmp*Sintemp(2)**3- &
966 0 : (dbgrid-c1)*pre_sin(2)**3))
967 :
968 0 : 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 0 : fluxm_b = p5*(dhtmp*C_s(nblyr+2)* Sintemp(nblyr+2)*dbgrid &
977 0 : + C_s(nblyr+1)*(pre_sin(nblyr+1) &
978 0 : + dhtmp * Sintemp(nblyr+1)*(c1-dbgrid)))
979 :
980 0 : fluxm_t = -p5 * (C_s(1) * pre_sin(1)*dbgrid &
981 0 : + C_s(2) * (pre_sin(2)*(c1-dbgrid) + dhtmp*Sintemp(2)))
982 :
983 0 : fluxm = fluxm_b + fluxm_t
984 :
985 0 : Ssum_corr = (-dh/hbri_old + p5*(dh/hbri_old)**2)*Ssum_old
986 0 : Ssum_corr_flux = dh*Ssum_old/hin_next + Ssum_corr
987 0 : Ssum_corr = Ssum_corr_flux
988 :
989 : fzsaln_g = fzsaln_g - hin_next * fluxg_b &
990 0 : * rhosi*p001/dts
991 :
992 : !approximate fluxes
993 : !fzsaln = fzsaln - hin_next * (fluxg &
994 : ! + fluxb + fluxm + fluxcorr + Ssum_corr_flux) &
995 : ! * rhosi*p001/dts
996 :
997 : fzsaln = fzsaln + (Ssum_old*hin_old - Ssum_new*hin_next) &
998 0 : * rhosi*p001/dts ! positive into the ocean
999 :
1000 0 : 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, &
1008 : fluxcorr, fluxb, fluxg, fluxm, &
1009 : hbrin, hbri_old)
1010 :
1011 : integer(kind=int_kind), intent(in) :: &
1012 : mint , & ! current iteration
1013 : mmax ! maximum number of iterations
1014 :
1015 : real (kind=dbl_kind), intent(in) :: &
1016 : dt , & ! thermodynamic and halodynamic timesteps (s)
1017 : hbrin , & ! (m) final brine height
1018 : hbri_old , & ! (m) initial brine height
1019 : Ssum_old , & ! initial integrated salt content
1020 : Ssum_new , & ! final integrated salt content
1021 : fluxcorr , & ! flux correction to ensure S >= min_salin
1022 : Ssum_corr , & ! boundary flux correction due to numerics
1023 : fluxb , & ! total boundary salt flux into the ice (+ into ice)
1024 : fluxg , & ! total gravity drainage salt flux into the ice (+ into ice)
1025 : fluxm !
1026 :
1027 : ! local variables
1028 :
1029 : real (kind=dbl_kind):: &
1030 0 : diff2 , & !
1031 0 : dsum_flux , & ! salt change in kg/m^2/s
1032 0 : flux_tot , & ! fluxg + fluxb
1033 0 : order , & !
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)* &
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 0 : subroutine merge_zsal_fluxes(aicenS, &
1083 : zsal_totn, zsal_tot, &
1084 : fzsal, fzsaln, &
1085 : fzsal_g, fzsaln_g)
1086 :
1087 : ! single category fluxes
1088 : real (kind=dbl_kind), intent(in):: &
1089 : aicenS , & ! concentration of ice
1090 : fzsaln , & ! salt flux (kg/m**2/s)
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)
1098 : fzsal , & ! salt flux (kg/m**2/s)
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 0 : zsal_tot = zsal_tot + zsal_totn ! already *aicenS
1108 :
1109 : ! ocean tot and gravity drainage salt fluxes
1110 0 : fzsal = fzsal + fzsaln * aicenS
1111 0 : fzsal_g = fzsal_g + fzsaln_g * aicenS
1112 :
1113 0 : 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 0 : subroutine column_sum_zsal (zsal_totn, nblyr, &
1121 0 : 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)
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 0 : do k = 1, nblyr
1144 : zsal_totn = zsal_totn &
1145 0 : + rhosi * trcrn_S(k) &
1146 : * fbri &
1147 0 : * vicenS/real(nblyr,kind=dbl_kind)
1148 : enddo ! k
1149 :
1150 0 : end subroutine column_sum_zsal
1151 :
1152 : !=======================================================================
1153 :
1154 : end module icepack_zsalinity
1155 :
1156 : !=======================================================================
|
