Line data Source code
1 : !=======================================================================
2 : !
3 : ! Thermo calculations after call to coupler, related to ITD:
4 : ! ice thickness redistribution, lateral growth and melting.
5 : !
6 : ! NOTE: The thermodynamic calculation is split in two for load balancing.
7 : ! First icepack_therm_vertical computes vertical growth rates and coupler
8 : ! fluxes. Then icepack_therm_itd does thermodynamic calculations not
9 : ! needed for coupling.
10 : !
11 : ! authors William H. Lipscomb, LANL
12 : ! C. M. Bitz, UW
13 : ! Elizabeth C. Hunke, LANL
14 : !
15 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
16 : ! 2004: Block structure added by William Lipscomb
17 : ! 2006: Streamlined for efficiency by Elizabeth Hunke
18 : ! 2014: Column package created by Elizabeth Hunke
19 : !
20 : module icepack_therm_itd
21 :
22 : use icepack_kinds
23 :
24 : use icepack_fsd, only: floe_rad_c, floe_binwidth
25 :
26 : use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10
27 : use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
28 : use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
29 : use icepack_parameters, only: phi_init, dsin0_frazil, salt_loss
30 : use icepack_parameters, only: Tliquidus_max
31 : use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist
32 : use icepack_parameters, only: kitd, ktherm
33 : use icepack_parameters, only: z_tracers, hfrazilmin, hi_min
34 : use icepack_parameters, only: cpl_frazil, update_ocn_f, saltflux_option
35 : use icepack_parameters, only: icepack_chkoptargflag
36 :
37 : use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, nfsd
38 : use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
39 : use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
40 : use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos, nt_sice
41 : use icepack_tracers, only: nt_alvl, nt_vlvl
42 : use icepack_tracers, only: tr_pond_lvl, tr_pond_topo
43 : use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
44 : use icepack_tracers, only: n_aero, n_iso
45 : use icepack_tracers, only: bio_index
46 :
47 : use icepack_warnings, only: warnstr, icepack_warnings_add
48 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
49 :
50 : use icepack_fsd, only: fsd_weld_thermo, icepack_cleanup_fsd, get_subdt_fsd
51 : use icepack_itd, only: reduce_area, cleanup_itd
52 : use icepack_itd, only: aggregate_area, shift_ice
53 : use icepack_itd, only: column_sum, column_conservation_check
54 : use icepack_isotope, only: isoice_alpha, isotope_frac_method
55 : use icepack_mushy_physics, only: liquidus_temperature_mush, icepack_enthalpy_mush
56 : use icepack_zbgc, only: add_new_ice_bgc
57 : use icepack_zbgc, only: lateral_melt_bgc
58 : use icepack_zbgc_shared, only: bgrid, cgrid, igrid
59 :
60 : implicit none
61 :
62 : private
63 : public :: icepack_step_therm2
64 :
65 : !=======================================================================
66 :
67 : contains
68 :
69 : !=======================================================================
70 : !
71 : ! ITD scheme that shifts ice among categories
72 : !
73 : ! See Lipscomb, W. H. Remapping the thickness distribution in sea
74 : ! ice models. 2001, J. Geophys. Res., Vol 106, 13989--14000.
75 : !
76 : ! Using the thermodynamic "velocities", interpolate to find the
77 : ! velocities in thickness space at the category boundaries, and
78 : ! compute the new locations of the boundaries. Then for each
79 : ! category, compute the thickness distribution function, g(h),
80 : ! between hL and hR, the left and right boundaries of the category.
81 : ! Assume g(h) is a linear polynomial that satisfies two conditions:
82 : !
83 : ! (1) The ice area implied by g(h) equals aicen(n).
84 : ! (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
85 : !
86 : ! Given g(h), at each boundary compute the ice area and volume lying
87 : ! between the original and new boundary locations. Transfer area
88 : ! and volume across each boundary in the appropriate direction, thus
89 : ! restoring the original boundaries.
90 : !
91 : ! authors: William H. Lipscomb, LANL
92 : ! Elizabeth C. Hunke, LANL
93 :
94 0 : subroutine linear_itd (hin_max, trcr_depend, &
95 1790590 : trcr_base, n_trcr_strata,&
96 1790590 : nt_strata, &
97 1790590 : aicen_init, vicen_init, &
98 3581180 : aicen, trcrn, &
99 1790590 : vicen, vsnon, &
100 : aice, aice0, &
101 : fpond, Tf )
102 :
103 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
104 : hin_max ! category boundaries (m)
105 :
106 : integer (kind=int_kind), dimension (:), intent(in) :: &
107 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
108 : n_trcr_strata ! number of underlying tracer layers
109 :
110 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
111 : trcr_base ! = 0 or 1 depending on tracer dependency
112 : ! argument 2: (1) aice, (2) vice, (3) vsno
113 :
114 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
115 : nt_strata ! indices of underlying tracer layers
116 :
117 : real (kind=dbl_kind), intent(in) :: &
118 : Tf ! freezing temperature
119 :
120 : real (kind=dbl_kind), dimension(:), intent(in) :: &
121 : aicen_init, & ! initial ice concentration (before vertical thermo)
122 : vicen_init ! initial ice volume (m)
123 :
124 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
125 : aicen , & ! ice concentration
126 : vicen , & ! volume per unit area of ice (m)
127 : vsnon ! volume per unit area of snow (m)
128 :
129 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
130 : trcrn ! ice tracers
131 :
132 : real (kind=dbl_kind), intent(inout) :: &
133 : aice , & ! concentration of ice
134 : aice0 , & ! concentration of open water
135 : fpond ! fresh water flux to ponds (kg/m^2/s)
136 :
137 : ! local variables
138 :
139 : integer (kind=int_kind) :: &
140 : n, nd , & ! category indices
141 : k ! ice layer index
142 :
143 : real (kind=dbl_kind) :: &
144 : slope , & ! rate of change of dhice with hice
145 : dh0 , & ! change in ice thickness at h = 0
146 : da0 , & ! area melting from category 1
147 : damax , & ! max allowed reduction in category 1 area
148 : etamin, etamax,& ! left and right limits of integration
149 : x1 , & ! etamax - etamin
150 : x2 , & ! (etamax^2 - etamin^2) / 2
151 : x3 , & ! (etamax^3 - etamin^3) / 3
152 : wk1, wk2 ! temporary variables
153 :
154 : real (kind=dbl_kind), dimension(0:ncat) :: &
155 3581180 : hbnew ! new boundary locations
156 :
157 : real (kind=dbl_kind), dimension(ncat) :: &
158 3581180 : g0 , & ! constant coefficient in g(h)
159 3581180 : g1 , & ! linear coefficient in g(h)
160 3581180 : hL , & ! left end of range over which g(h) > 0
161 3581180 : hR ! right end of range over which g(h) > 0
162 :
163 : real (kind=dbl_kind), dimension(ncat) :: &
164 3581180 : hicen , & ! ice thickness for each cat (m)
165 3581180 : hicen_init , & ! initial ice thickness for each cat (pre-thermo)
166 3581180 : dhicen , & ! thickness change for remapping (m)
167 1790590 : daice , & ! ice area transferred across boundary
168 3581180 : dvice ! ice volume transferred across boundary
169 :
170 : real (kind=dbl_kind), dimension (ncat) :: &
171 3581180 : eicen, & ! energy of melting for each ice layer (J/m^2)
172 3581180 : esnon, & ! energy of melting for each snow layer (J/m^2)
173 3581180 : vbrin, & ! ice volume defined by brine height (m)
174 3581180 : sicen ! Bulk salt in h ice (ppt*m)
175 :
176 : real (kind=dbl_kind) :: &
177 : vice_init, vice_final, & ! ice volume summed over categories
178 : vsno_init, vsno_final, & ! snow volume summed over categories
179 : eice_init, eice_final, & ! ice energy summed over categories
180 : esno_init, esno_final, & ! snow energy summed over categories
181 : sice_init, sice_final, & ! ice bulk salinity summed over categories
182 : vbri_init, vbri_final ! briny ice volume summed over categories
183 :
184 : ! NOTE: Third index of donor, daice, dvice should be ncat-1,
185 : ! except that compilers would have trouble when ncat = 1
186 : integer (kind=int_kind), dimension(ncat) :: &
187 1790590 : donor ! donor category index
188 :
189 : logical (kind=log_kind) :: &
190 : remap_flag ! remap ITD if remap_flag is true
191 :
192 : character (len=char_len) :: &
193 : fieldid ! field identifier
194 :
195 : logical (kind=log_kind), parameter :: &
196 : print_diags = .false. ! if true, prints when remap_flag=F
197 :
198 : character(len=*),parameter :: subname='(linear_itd)'
199 :
200 : !-----------------------------------------------------------------
201 : ! Initialize
202 : !-----------------------------------------------------------------
203 :
204 10743540 : do n = 1, ncat
205 8952950 : donor(n) = 0
206 8952950 : daice(n) = c0
207 10743540 : dvice(n) = c0
208 : enddo
209 :
210 : !-----------------------------------------------------------------
211 : ! Compute volume and energy sums that linear remapping should
212 : ! conserve.
213 : !-----------------------------------------------------------------
214 :
215 1790590 : if (conserv_check) then
216 :
217 289584 : do n = 1, ncat
218 :
219 241320 : eicen(n) = c0
220 241320 : esnon(n) = c0
221 241320 : vbrin(n) = c0
222 241320 : sicen(n) = c0
223 :
224 1930560 : do k = 1, nilyr
225 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
226 1930560 : * vicen(n)/real(nilyr,kind=dbl_kind)
227 : enddo
228 1447920 : do k = 1, nslyr
229 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
230 1447920 : * vsnon(n)/real(nslyr,kind=dbl_kind)
231 : enddo
232 :
233 241320 : if (tr_brine) then
234 : vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
235 0 : * vicen(n)
236 : endif
237 :
238 1978824 : do k = 1, nilyr
239 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
240 1930560 : * vicen(n)/real(nilyr,kind=dbl_kind)
241 : enddo
242 :
243 : enddo ! n
244 :
245 48264 : call column_sum (ncat, vicen, vice_init)
246 48264 : if (icepack_warnings_aborted(subname)) return
247 48264 : call column_sum (ncat, vsnon, vsno_init)
248 48264 : if (icepack_warnings_aborted(subname)) return
249 48264 : call column_sum (ncat, eicen, eice_init)
250 48264 : if (icepack_warnings_aborted(subname)) return
251 48264 : call column_sum (ncat, esnon, esno_init)
252 48264 : if (icepack_warnings_aborted(subname)) return
253 48264 : call column_sum (ncat, sicen, sice_init)
254 48264 : if (icepack_warnings_aborted(subname)) return
255 48264 : call column_sum (ncat, vbrin, vbri_init)
256 48264 : if (icepack_warnings_aborted(subname)) return
257 :
258 : endif ! conserv_check
259 :
260 : !-----------------------------------------------------------------
261 : ! Initialize remapping flag.
262 : ! Remapping is done wherever remap_flag = .true.
263 : ! In rare cases the category boundaries may shift too far for the
264 : ! remapping algorithm to work, and remap_flag is set to .false.
265 : ! In these cases the simpler 'rebin' subroutine will shift ice
266 : ! between categories if needed.
267 : !-----------------------------------------------------------------
268 :
269 1790590 : remap_flag = .true.
270 :
271 : !-----------------------------------------------------------------
272 : ! Compute thickness change in each category.
273 : !-----------------------------------------------------------------
274 :
275 10743540 : do n = 1, ncat
276 :
277 8952950 : if (aicen_init(n) > puny) then
278 8371988 : hicen_init(n) = vicen_init(n) / aicen_init(n)
279 : else
280 580962 : hicen_init(n) = c0
281 : endif ! aicen_init > puny
282 :
283 10743540 : if (aicen (n) > puny) then
284 8371984 : hicen (n) = vicen(n) / aicen(n)
285 8371984 : dhicen(n) = hicen(n) - hicen_init(n)
286 : else
287 580966 : hicen (n) = c0
288 580966 : dhicen(n) = c0
289 : endif ! aicen > puny
290 :
291 : enddo ! n
292 :
293 : !-----------------------------------------------------------------
294 : ! Compute new category boundaries, hbnew, based on changes in
295 : ! ice thickness from vertical thermodynamics.
296 : !-----------------------------------------------------------------
297 :
298 1790590 : hbnew(0) = hin_max(0)
299 :
300 8952950 : do n = 1, ncat-1
301 :
302 7162360 : if (hicen_init(n) > puny .and. &
303 : hicen_init(n+1) > puny) then
304 :
305 6580990 : if ((hicen_init(n+1) - hicen_init(n))>0) then
306 :
307 : ! interpolate between adjacent category growth rates
308 : slope = (dhicen(n+1) - dhicen(n)) / &
309 6580990 : (hicen_init(n+1) - hicen_init(n))
310 : hbnew(n) = hin_max(n) + dhicen(n) &
311 6580990 : + slope * (hin_max(n) - hicen_init(n))
312 :
313 : else
314 :
315 0 : write(warnstr,*) subname, &
316 0 : 'ITD Thermodynamics: hicen_init(n+1) <= hicen_init(n)'
317 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
318 0 : call icepack_warnings_add(warnstr)
319 :
320 : endif
321 :
322 581370 : elseif (hicen_init(n) > puny) then ! hicen_init(n+1)=0
323 123902 : hbnew(n) = hin_max(n) + dhicen(n)
324 457468 : elseif (hicen_init(n+1) > puny) then ! hicen_init(n)=0
325 295756 : hbnew(n) = hin_max(n) + dhicen(n+1)
326 : else
327 161712 : hbnew(n) = hin_max(n)
328 : endif
329 :
330 : !-----------------------------------------------------------------
331 : ! Check that each boundary lies between adjacent values of hicen.
332 : ! If not, set remap_flag = .false.
333 : ! Write diagnosis outputs if remap_flag was changed to false
334 : !-----------------------------------------------------------------
335 :
336 7162360 : if (aicen(n) > puny .and. hicen(n) >= hbnew(n)) then
337 2 : remap_flag = .false.
338 :
339 : if (print_diags) then
340 : write(warnstr,*) subname, 'ITD: hicen(n) > hbnew(n)'
341 : call icepack_warnings_add(warnstr)
342 : write(warnstr,*) subname, 'cat ',n
343 : call icepack_warnings_add(warnstr)
344 : write(warnstr,*) subname, 'hicen(n) =', hicen(n)
345 : call icepack_warnings_add(warnstr)
346 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
347 : call icepack_warnings_add(warnstr)
348 : endif
349 :
350 7162358 : elseif (aicen(n+1) > puny .and. hicen(n+1) <= hbnew(n)) then
351 0 : remap_flag = .false.
352 :
353 : if (print_diags) then
354 : write(warnstr,*) subname, 'ITD: hicen(n+1) < hbnew(n)'
355 : call icepack_warnings_add(warnstr)
356 : write(warnstr,*) subname, 'cat ',n
357 : call icepack_warnings_add(warnstr)
358 : write(warnstr,*) subname, 'hicen(n+1) =', hicen(n+1)
359 : call icepack_warnings_add(warnstr)
360 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
361 : call icepack_warnings_add(warnstr)
362 : endif
363 : endif
364 :
365 : !-----------------------------------------------------------------
366 : ! Check that hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
367 : ! If not, set remap_flag = .false.
368 : ! (In principle we could allow this, but it would make the code
369 : ! more complicated.)
370 : ! Write diagnosis outputs if remap_flag was changed to false
371 : !-----------------------------------------------------------------
372 :
373 7162360 : if (hbnew(n) > hin_max(n+1)) then
374 0 : remap_flag = .false.
375 :
376 : if (print_diags) then
377 : write(warnstr,*) subname, 'ITD hbnew(n) > hin_max(n+1)'
378 : call icepack_warnings_add(warnstr)
379 : write(warnstr,*) subname, 'cat ',n
380 : call icepack_warnings_add(warnstr)
381 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
382 : call icepack_warnings_add(warnstr)
383 : write(warnstr,*) subname, 'hin_max(n+1) =', hin_max(n+1)
384 : call icepack_warnings_add(warnstr)
385 : endif
386 : endif
387 :
388 8952950 : if (hbnew(n) < hin_max(n-1)) then
389 0 : remap_flag = .false.
390 :
391 : if (print_diags) then
392 : write(warnstr,*) subname, 'ITD: hbnew(n) < hin_max(n-1)'
393 : call icepack_warnings_add(warnstr)
394 : write(warnstr,*) subname, 'cat ',n
395 : call icepack_warnings_add(warnstr)
396 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
397 : call icepack_warnings_add(warnstr)
398 : write(warnstr,*) subname, 'hin_max(n-1) =', hin_max(n-1)
399 : call icepack_warnings_add(warnstr)
400 : endif
401 : endif
402 :
403 : enddo ! boundaries, 1 to ncat-1
404 :
405 : !-----------------------------------------------------------------
406 : ! Compute hbnew(ncat)
407 : !-----------------------------------------------------------------
408 :
409 1790590 : if (aicen(ncat) > puny) then
410 1667096 : hbnew(ncat) = c3*hicen(ncat) - c2*hbnew(ncat-1)
411 : else
412 123494 : hbnew(ncat) = hin_max(ncat)
413 : endif
414 1790590 : hbnew(ncat) = max(hbnew(ncat),hin_max(ncat-1))
415 :
416 : !-----------------------------------------------------------------
417 : ! Identify cells where the ITD is to be remapped
418 : !-----------------------------------------------------------------
419 :
420 1790590 : if (remap_flag) then
421 :
422 : !-----------------------------------------------------------------
423 : ! Compute g(h) for category 1 at start of time step
424 : ! (hicen = hicen_init)
425 : !-----------------------------------------------------------------
426 :
427 : call fit_line(aicen(1), hicen_init(1), &
428 : hbnew(0), hin_max (1), &
429 : g0 (1), g1 (1), &
430 1790588 : hL (1), hR (1))
431 1790588 : if (icepack_warnings_aborted(subname)) return
432 :
433 : !-----------------------------------------------------------------
434 : ! Find area lost due to melting of thin (category 1) ice
435 : !-----------------------------------------------------------------
436 :
437 1790588 : if (aicen(1) > puny) then
438 :
439 1495238 : dh0 = dhicen(1)
440 1495238 : if (dh0 < c0) then ! remove area from category 1
441 376514 : dh0 = min(-dh0,hin_max(1)) ! dh0 --> |dh0|
442 :
443 : !-----------------------------------------------------------------
444 : ! Integrate g(1) from 0 to dh0 to estimate area melted
445 : !-----------------------------------------------------------------
446 :
447 : ! right integration limit (left limit = 0)
448 376514 : etamax = min(dh0,hR(1)) - hL(1)
449 :
450 376514 : if (etamax > c0) then
451 241701 : x1 = etamax
452 241701 : x2 = p5 * etamax*etamax
453 241701 : da0 = g1(1)*x2 + g0(1)*x1 ! ice area removed
454 :
455 : ! constrain new thickness <= hicen_init
456 241701 : damax = aicen(1) * (c1-hicen(1)/hicen_init(1)) ! damax > 0
457 241701 : da0 = min (da0, damax)
458 :
459 : ! remove area, conserving volume
460 241701 : hicen(1) = hicen(1) * aicen(1) / (aicen(1)-da0)
461 241701 : aicen(1) = aicen(1) - da0
462 :
463 241701 : if (tr_pond_topo) &
464 : fpond = fpond - (da0 * trcrn(nt_apnd,1) &
465 46889 : * trcrn(nt_hpnd,1))
466 :
467 : endif ! etamax > 0
468 :
469 : else ! dh0 >= 0
470 1118724 : hbnew(0) = min(dh0,hin_max(1)) ! shift hbnew(0) to right
471 : endif
472 :
473 : endif ! aicen(1) > puny
474 :
475 : !-----------------------------------------------------------------
476 : ! Compute g(h) for each ice thickness category.
477 : !-----------------------------------------------------------------
478 :
479 10743528 : do n = 1, ncat
480 :
481 : call fit_line(aicen(n), hicen(n), &
482 : hbnew(n-1), hbnew(n), &
483 : g0 (n), g1 (n), &
484 8952940 : hL (n), hR (n))
485 8952940 : if (icepack_warnings_aborted(subname)) return
486 :
487 : !-----------------------------------------------------------------
488 : ! Compute area and volume to be shifted across each boundary.
489 : !-----------------------------------------------------------------
490 :
491 8952940 : donor(n) = 0
492 8952940 : daice(n) = c0
493 10743528 : dvice(n) = c0
494 : enddo
495 :
496 8952940 : do n = 1, ncat-1
497 :
498 7162352 : if (hbnew(n) > hin_max(n)) then ! transfer from n to n+1
499 :
500 : ! left and right integration limits in eta space
501 4193543 : etamin = max(hin_max(n), hL(n)) - hL(n)
502 4193543 : etamax = min(hbnew(n), hR(n)) - hL(n)
503 4193543 : donor(n) = n
504 :
505 : else ! hbnew(n) <= hin_max(n); transfer from n+1 to n
506 :
507 : ! left and right integration limits in eta space
508 2968809 : etamin = c0
509 2968809 : etamax = min(hin_max(n), hR(n+1)) - hL(n+1)
510 2968809 : donor(n) = n+1
511 :
512 : endif ! hbnew(n) > hin_max(n)
513 :
514 7162352 : if (etamax > etamin) then
515 5852358 : x1 = etamax - etamin
516 5852358 : wk1 = etamin*etamin
517 5852358 : wk2 = etamax*etamax
518 5852358 : x2 = p5 * (wk2 - wk1)
519 5852358 : wk1 = wk1*etamin
520 5852358 : wk2 = wk2*etamax
521 5852358 : x3 = p333 * (wk2 - wk1)
522 5852358 : nd = donor(n)
523 5852358 : daice(n) = g1(nd)*x2 + g0(nd)*x1
524 5852358 : dvice(n) = g1(nd)*x3 + g0(nd)*x2 + daice(n)*hL(nd)
525 : endif ! etamax > etamin
526 :
527 : ! If daice or dvice is very small, shift no ice.
528 :
529 7162352 : nd = donor(n)
530 :
531 7162352 : if (daice(n) < aicen(nd)*puny) then
532 924050 : daice(n) = c0
533 924050 : dvice(n) = c0
534 924050 : donor(n) = 0
535 : endif
536 :
537 7162352 : if (dvice(n) < vicen(nd)*puny) then
538 924050 : daice(n) = c0
539 924050 : dvice(n) = c0
540 924050 : donor(n) = 0
541 : endif
542 :
543 : ! If daice is close to aicen or dvice is close to vicen,
544 : ! shift entire category
545 :
546 7162352 : if (daice(n) > aicen(nd)*(c1-puny)) then
547 6073 : daice(n) = aicen(nd)
548 6073 : dvice(n) = vicen(nd)
549 : endif
550 :
551 8952940 : if (dvice(n) > vicen(nd)*(c1-puny)) then
552 6073 : daice(n) = aicen(nd)
553 6073 : dvice(n) = vicen(nd)
554 : endif
555 :
556 : enddo ! boundaries, 1 to ncat-1
557 :
558 : !-----------------------------------------------------------------
559 : ! Shift ice between categories as necessary
560 : !-----------------------------------------------------------------
561 :
562 : ! maintain qsno negative definiteness
563 10743528 : do n = 1, ncat
564 25528668 : do k = nt_qsno, nt_qsno+nslyr-1
565 23738080 : trcrn(k,n) = trcrn(k,n) + rhos*Lfresh
566 : enddo
567 : enddo
568 : ! maintain rhos_cmp positive definiteness
569 1790588 : if (snwredist(1:3) == 'ITD') then
570 1302444 : do n = 1, ncat
571 6729294 : do k = nt_rhos, nt_rhos+nslyr-1
572 6512220 : trcrn(k,n) = max(trcrn(k,n)-rhosmin, c0)
573 : ! trcrn(k,n) = trcrn(k,n) - rhosmin
574 : enddo
575 : enddo
576 : endif
577 :
578 : call shift_ice (trcr_depend, &
579 : trcr_base, &
580 : n_trcr_strata, &
581 : nt_strata, &
582 : aicen, trcrn, &
583 : vicen, vsnon, &
584 : hicen, donor, &
585 1790588 : daice, dvice, Tf )
586 1790588 : if (icepack_warnings_aborted(subname)) return
587 :
588 : ! maintain qsno negative definiteness
589 10743528 : do n = 1, ncat
590 25528668 : do k = nt_qsno, nt_qsno+nslyr-1
591 23738080 : trcrn(k,n) = trcrn(k,n) - rhos*Lfresh
592 : enddo
593 : enddo
594 : ! maintain rhos_cmp positive definiteness
595 1790588 : if (snwredist(1:3) == 'ITD') then
596 1302444 : do n = 1, ncat
597 6729294 : do k = nt_rhos, nt_rhos+nslyr-1
598 6512220 : trcrn(k,n) = trcrn(k,n) + rhosmin
599 : enddo
600 : enddo
601 : endif
602 :
603 : !-----------------------------------------------------------------
604 : ! Make sure hice(1) >= minimum ice thickness hi_min.
605 : !-----------------------------------------------------------------
606 :
607 1790588 : if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then
608 :
609 1721 : da0 = aicen(1) * (c1 - hicen(1)/hi_min)
610 1721 : aicen(1) = aicen(1) - da0
611 1721 : hicen(1) = hi_min
612 :
613 1721 : if (tr_pond_topo) &
614 : fpond = fpond - (da0 * trcrn(nt_apnd,1) &
615 0 : * trcrn(nt_hpnd,1))
616 : endif
617 :
618 : endif ! remap_flag
619 :
620 : !-----------------------------------------------------------------
621 : ! Update fractional ice area in each grid cell.
622 : !-----------------------------------------------------------------
623 :
624 1790590 : call aggregate_area (aicen, aice, aice0)
625 1790590 : if (icepack_warnings_aborted(subname)) return
626 :
627 : !-----------------------------------------------------------------
628 : ! Check volume and energy conservation.
629 : !-----------------------------------------------------------------
630 :
631 1790590 : if (conserv_check) then
632 :
633 289584 : do n = 1, ncat
634 :
635 241320 : eicen(n) = c0
636 241320 : esnon(n) = c0
637 241320 : vbrin(n) = c0
638 241320 : sicen(n) = c0
639 :
640 1930560 : do k = 1, nilyr
641 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
642 1930560 : * vicen(n)/real(nilyr,kind=dbl_kind)
643 : enddo
644 1447920 : do k = 1, nslyr
645 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
646 1447920 : * vsnon(n)/real(nslyr,kind=dbl_kind)
647 : enddo
648 :
649 241320 : if (tr_brine) then
650 : vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
651 0 : * vicen(n)
652 : endif
653 :
654 1978824 : do k = 1, nilyr
655 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
656 1930560 : * vicen(n)/real(nilyr,kind=dbl_kind)
657 : enddo
658 :
659 : enddo ! n
660 :
661 48264 : call column_sum (ncat, vicen, vice_final)
662 48264 : if (icepack_warnings_aborted(subname)) return
663 48264 : call column_sum (ncat, vsnon, vsno_final)
664 48264 : if (icepack_warnings_aborted(subname)) return
665 48264 : call column_sum (ncat, eicen, eice_final)
666 48264 : if (icepack_warnings_aborted(subname)) return
667 48264 : call column_sum (ncat, esnon, esno_final)
668 48264 : if (icepack_warnings_aborted(subname)) return
669 48264 : call column_sum (ncat, sicen, sice_final)
670 48264 : if (icepack_warnings_aborted(subname)) return
671 48264 : call column_sum (ncat, vbrin, vbri_final)
672 48264 : if (icepack_warnings_aborted(subname)) return
673 :
674 48264 : fieldid = subname//':vice'
675 : call column_conservation_check (fieldid, &
676 : vice_init, vice_final, &
677 48264 : puny)
678 48264 : if (icepack_warnings_aborted(subname)) return
679 48264 : fieldid = subname//':vsno'
680 : call column_conservation_check (fieldid, &
681 : vsno_init, vsno_final, &
682 48264 : puny)
683 48264 : if (icepack_warnings_aborted(subname)) return
684 48264 : fieldid = subname//':eice'
685 : call column_conservation_check (fieldid, &
686 : eice_init, eice_final, &
687 48264 : puny*Lfresh*rhoi)
688 48264 : if (icepack_warnings_aborted(subname)) return
689 48264 : fieldid = subname//':esno'
690 : call column_conservation_check (fieldid, &
691 : esno_init, esno_final, &
692 48264 : puny*Lfresh*rhos)
693 48264 : if (icepack_warnings_aborted(subname)) return
694 48264 : fieldid = subname//':sicen'
695 : call column_conservation_check (fieldid, &
696 : sice_init, sice_final, &
697 48264 : puny)
698 48264 : if (icepack_warnings_aborted(subname)) return
699 48264 : fieldid = subname//':vbrin'
700 : call column_conservation_check (fieldid, &
701 : vbri_init, vbri_final, &
702 48264 : puny*c10)
703 48264 : if (icepack_warnings_aborted(subname)) return
704 :
705 : endif ! conservation check
706 :
707 : end subroutine linear_itd
708 :
709 : !=======================================================================
710 : !
711 : ! Fit g(h) with a line, satisfying area and volume constraints.
712 : ! To reduce roundoff errors caused by large values of g0 and g1,
713 : ! we actually compute g(eta), where eta = h - hL, and hL is the
714 : ! left boundary.
715 : !
716 : ! authors: William H. Lipscomb, LANL
717 : ! Elizabeth C. Hunke, LANL
718 :
719 10743528 : subroutine fit_line (aicen, hice, &
720 : hbL, hbR, &
721 : g0, g1, &
722 : hL, hR)
723 :
724 : real (kind=dbl_kind), intent(in) :: &
725 : aicen ! concentration of ice
726 :
727 : real (kind=dbl_kind), intent(in) :: &
728 : hbL, hbR , & ! left and right category boundaries
729 : hice ! ice thickness
730 :
731 : real (kind=dbl_kind), intent(out):: &
732 : g0, g1 , & ! coefficients in linear equation for g(eta)
733 : hL , & ! min value of range over which g(h) > 0
734 : hR ! max value of range over which g(h) > 0
735 :
736 : ! local variables
737 :
738 : real (kind=dbl_kind) :: &
739 : h13 , & ! hbL + 1/3 * (hbR - hbL)
740 : h23 , & ! hbL + 2/3 * (hbR - hbL)
741 : dhr , & ! 1 / (hR - hL)
742 : wk1, wk2 ! temporary variables
743 :
744 : character(len=*),parameter :: subname='(fit_line)'
745 :
746 : !-----------------------------------------------------------------
747 : ! Compute g0, g1, hL, and hR for each category to be remapped.
748 : !-----------------------------------------------------------------
749 :
750 10743528 : if (aicen > puny .and. hbR - hbL > puny) then
751 :
752 : ! Initialize hL and hR
753 :
754 9867214 : hL = hbL
755 9867214 : hR = hbR
756 :
757 : ! Change hL or hR if hicen(n) falls outside central third of range
758 :
759 9867214 : h13 = p333 * (c2*hL + hR)
760 9867214 : h23 = p333 * (hL + c2*hR)
761 9867214 : if (hice < h13) then
762 1437319 : hR = c3*hice - c2*hL
763 8429895 : elseif (hice > h23) then
764 1660779 : hL = c3*hice - c2*hR
765 : endif
766 :
767 : ! Compute coefficients of g(eta) = g0 + g1*eta
768 :
769 9867214 : dhr = c1 / (hR - hL)
770 9867214 : wk1 = c6 * aicen * dhr
771 9867214 : wk2 = (hice - hL) * dhr
772 9867214 : g0 = wk1 * (p666 - wk2)
773 9867214 : g1 = c2*dhr * wk1 * (wk2 - p5)
774 :
775 : else
776 :
777 876314 : g0 = c0
778 876314 : g1 = c0
779 876314 : hL = c0
780 876314 : hR = c0
781 :
782 : endif ! aicen > puny
783 :
784 10743528 : end subroutine fit_line
785 :
786 : !=======================================================================
787 : !
788 : ! Given some added new ice to the base of the existing ice, recalculate
789 : ! vertical tracer so that new grid cells are all the same size.
790 : !
791 : ! author: A. K. Turner, LANL
792 : !
793 328456 : subroutine update_vertical_tracers(trc, h1, h2, trc0)
794 :
795 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
796 : trc ! vertical tracer
797 :
798 : real (kind=dbl_kind), intent(in) :: &
799 : h1, & ! old thickness
800 : h2, & ! new thickness
801 : trc0 ! tracer value of added ice on ice bottom
802 :
803 : ! local variables
804 :
805 656912 : real(kind=dbl_kind), dimension(nilyr) :: trc2 ! updated tracer temporary
806 :
807 : ! vertical indices for old and new grid
808 : integer :: k1, k2
809 :
810 : real (kind=dbl_kind) :: &
811 : z1a, z1b, & ! upper, lower boundary of old cell/added new ice at bottom
812 : z2a, z2b, & ! upper, lower boundary of new cell
813 : overlap , & ! overlap between old and new cell
814 : rnilyr
815 :
816 : character(len=*),parameter :: subname='(update_vertical_tracers)'
817 :
818 328456 : rnilyr = real(nilyr,dbl_kind)
819 :
820 : ! loop over new grid cells
821 2627648 : do k2 = 1, nilyr
822 :
823 : ! initialize new tracer
824 2299192 : trc2(k2) = c0
825 :
826 : ! calculate upper and lower boundary of new cell
827 2299192 : z2a = ((k2 - 1) * h2) / rnilyr
828 2299192 : z2b = (k2 * h2) / rnilyr
829 :
830 : ! loop over old grid cells
831 18393536 : do k1 = 1, nilyr
832 :
833 : ! calculate upper and lower boundary of old cell
834 16094344 : z1a = ((k1 - 1) * h1) / rnilyr
835 16094344 : z1b = (k1 * h1) / rnilyr
836 :
837 : ! calculate overlap between old and new cell
838 16094344 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
839 :
840 : ! aggregate old grid cell contribution to new cell
841 18393536 : trc2(k2) = trc2(k2) + overlap * trc(k1)
842 :
843 : enddo ! k1
844 :
845 : ! calculate upper and lower boundary of added new ice at bottom
846 2299192 : z1a = h1
847 2299192 : z1b = h2
848 :
849 : ! calculate overlap between added ice and new cell
850 2299192 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
851 : ! aggregate added ice contribution to new cell
852 2299192 : trc2(k2) = trc2(k2) + overlap * trc0
853 : ! renormalize new grid cell
854 2627648 : trc2(k2) = (rnilyr * trc2(k2)) / h2
855 :
856 : enddo ! k2
857 :
858 : ! update vertical tracer array with the adjusted tracer
859 2627648 : trc = trc2
860 :
861 328456 : end subroutine update_vertical_tracers
862 :
863 : !=======================================================================
864 : !
865 : ! Given the fraction of ice melting laterally in each grid cell
866 : ! (computed in subroutine frzmlt_bottom_lateral), melt ice.
867 : !
868 : ! author: C. M. Bitz, UW
869 : ! 2003: Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL
870 : ! 2016 Lettie Roach, NIWA/VUW, added floe size dependence
871 : !
872 1972116 : subroutine lateral_melt (dt, fpond, &
873 : fresh, fsalt, &
874 0 : fhocn, faero_ocn, &
875 1972116 : fiso_ocn, &
876 1972116 : rsiden, meltl, &
877 : wlat, &
878 3944232 : aicen, vicen, &
879 1972116 : vsnon, trcrn, &
880 1972116 : flux_bio, d_afsd_latm)
881 :
882 : real (kind=dbl_kind), intent(in) :: &
883 : dt ! time step (s)
884 :
885 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
886 : aicen , & ! concentration of ice
887 : vicen , & ! volume per unit area of ice (m)
888 : vsnon ! volume per unit area of snow (m)
889 :
890 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
891 : trcrn ! tracer array
892 :
893 : real (kind=dbl_kind), dimension (:), intent(in) :: &
894 : rsiden ! fraction of ice that melts laterally
895 :
896 : real (kind=dbl_kind), intent(in), optional :: &
897 : wlat ! lateral melt rate (m/s)
898 :
899 : real (kind=dbl_kind), intent(inout) :: &
900 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
901 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
902 : fsalt , & ! salt flux to ocean (kg/m^2/s)
903 : fhocn , & ! net heat flux to ocean (W/m^2)
904 : meltl ! lateral ice melt (m/step-->cm/day)
905 :
906 : real (kind=dbl_kind), dimension(nbtrcr), intent(inout) :: &
907 : flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s)
908 :
909 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
910 : faero_ocn ! aerosol flux to ocean (kg/m^2/s)
911 :
912 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
913 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
914 :
915 : real (kind=dbl_kind), dimension (:), intent(out), optional :: &
916 : d_afsd_latm ! change in fsd due to lateral melt (m)
917 :
918 : ! local variables
919 :
920 : integer (kind=int_kind) :: &
921 : n , & ! thickness category index
922 : k , & ! layer index
923 : nsubt ! sub timesteps for FSD tendency
924 :
925 : real (kind=dbl_kind) :: &
926 : dfhocn , & ! change in fhocn
927 : dfpond , & ! change in fpond
928 : dfresh , & ! change in fresh
929 : dfsalt , & ! change in fsalt
930 : dvssl , & ! snow surface layer volume
931 : dvint , & ! snow interior layer
932 : tmp
933 :
934 : real (kind=dbl_kind), dimension (ncat) :: &
935 3944232 : aicen_init, & ! initial area fraction
936 3944232 : vicen_init, & ! volume per unit area of ice (m)
937 3944232 : vsnon_init, & ! volume per unit area of snow (m)
938 3944232 : G_radialn ! rate of lateral melt (m/s)
939 :
940 : real (kind=dbl_kind), dimension (:,:), allocatable :: &
941 1972116 : afsdn , & ! floe size distribution tracer
942 1972116 : afsdn_init ! initial value
943 :
944 : real (kind=dbl_kind), dimension (:), allocatable :: &
945 1972116 : df_flx , & ! finite difference for FSD
946 1972116 : afsd_tmp , & !
947 1972116 : d_afsd_tmp, & !
948 1972116 : f_flx !
949 :
950 : real (kind=dbl_kind) :: &
951 : sicen, &
952 : etot, & ! column energy per itd cat, for FSD code
953 : elapsed_t, & ! FSD subcycling
954 : subdt ! FSD timestep (s)
955 :
956 : character(len=*), parameter :: subname='(lateral_melt)'
957 :
958 2867148 : if (tr_fsd) d_afsd_latm = c0
959 :
960 5385160 : if (.not. any(rsiden(:) > c0)) return ! no lateral melt, get out now
961 :
962 1232072 : dfhocn = c0
963 1232072 : dfpond = c0
964 1232072 : dfresh = c0
965 1232072 : dfsalt = c0
966 1232072 : dvssl = c0
967 1232072 : dvint = c0
968 7390024 : vicen_init(:) = vicen(:)
969 7390024 : vsnon_init(:) = vsnon(:)
970 :
971 1232072 : if (tr_fsd) then
972 60510 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:))
973 60510 : if (icepack_warnings_aborted(subname)) return
974 :
975 60510 : allocate(afsdn(nfsd,ncat))
976 60510 : allocate(afsdn_init(nfsd,ncat))
977 60510 : allocate(df_flx(nfsd))
978 60510 : allocate(afsd_tmp(nfsd))
979 60510 : allocate(d_afsd_tmp(nfsd))
980 60510 : allocate(f_flx(nfsd+1))
981 :
982 363060 : aicen_init = aicen
983 3062960 : afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:)
984 3062960 : afsdn_init = afsdn ! for diagnostics
985 588388 : df_flx = c0
986 648898 : f_flx = c0
987 : end if
988 :
989 7390024 : do n = 1, ncat
990 :
991 : !-----------------------------------------------------------------
992 : ! Melt the ice and increment fluxes.
993 : !-----------------------------------------------------------------
994 :
995 : ! fluxes to coupler
996 : ! dfresh > 0, dfsalt > 0, dfpond > 0
997 :
998 6157952 : dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt
999 6157952 : if (saltflux_option == 'prognostic') then
1000 134250 : sicen = c0
1001 1074000 : do k=1,nilyr
1002 1074000 : sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
1003 : enddo
1004 134250 : dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt
1005 : else
1006 6023702 : dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
1007 : endif
1008 6157952 : fresh = fresh + dfresh
1009 6157952 : fsalt = fsalt + dfsalt
1010 :
1011 6157952 : if (tr_pond_topo) then
1012 378640 : dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
1013 378640 : fpond = fpond - dfpond
1014 : endif
1015 :
1016 : ! history diagnostics
1017 6157952 : meltl = meltl + vicen_init(n)*rsiden(n)
1018 :
1019 : ! state variables
1020 6157952 : aicen(n) = aicen(n) * (c1 - rsiden(n))
1021 6157952 : vicen(n) = vicen(n) * (c1 - rsiden(n))
1022 6157952 : vsnon(n) = vsnon(n) * (c1 - rsiden(n))
1023 :
1024 : ! floe size distribution
1025 6157952 : if (tr_fsd) then
1026 302550 : if (rsiden(n) > puny) then
1027 297633 : if (aicen(n) > puny) then
1028 :
1029 : ! adaptive subtimestep
1030 282592 : elapsed_t = c0
1031 2762379 : afsd_tmp(:) = afsdn_init(:,n)
1032 2762379 : d_afsd_tmp(:) = c0
1033 282592 : nsubt = 0
1034 :
1035 565184 : DO WHILE (elapsed_t.lt.dt)
1036 :
1037 282592 : nsubt = nsubt + 1
1038 282592 : if (nsubt.gt.100) then
1039 0 : write(warnstr,*) subname, 'latm not converging'
1040 0 : call icepack_warnings_add(warnstr)
1041 : endif
1042 :
1043 : ! finite differences
1044 2762379 : df_flx(:) = c0
1045 3044971 : f_flx (:) = c0
1046 282592 : G_radialn(n) = -wlat
1047 2479787 : do k = 2, nfsd
1048 2396940 : f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
1049 : end do
1050 :
1051 2762379 : do k = 1, nfsd
1052 2762379 : df_flx(k) = f_flx(k+1) - f_flx(k)
1053 : end do
1054 :
1055 2762379 : if (abs(sum(df_flx(:))) > puny) then
1056 0 : write(warnstr,*) subname, 'sum(df_flx) /= 0'
1057 0 : call icepack_warnings_add(warnstr)
1058 : endif
1059 :
1060 : ! this term ensures area conservation
1061 2762379 : tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
1062 :
1063 : ! fsd tendency
1064 2762379 : do k = 1, nfsd
1065 : d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
1066 2762379 : * (c1/floe_rad_c(k) - tmp)
1067 : end do
1068 :
1069 : ! timestep required for this
1070 282592 : subdt = get_subdt_fsd(afsd_tmp(:), d_afsd_tmp(:))
1071 282592 : subdt = MIN(subdt, dt)
1072 :
1073 : ! update fsd and elapsed time
1074 2762379 : afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
1075 282592 : elapsed_t = elapsed_t + subdt
1076 :
1077 : END DO
1078 :
1079 2762379 : afsdn(:,n) = afsd_tmp(:)
1080 :
1081 : end if ! aicen
1082 : end if ! rside > 0, otherwise do nothing
1083 :
1084 : end if ! tr_fsd
1085 :
1086 : ! fluxes
1087 47842864 : do k = 1, nilyr
1088 : ! enthalpy tracers do not change (e/v constant)
1089 : ! heat flux to coupler for ice melt (dfhocn < 0)
1090 : dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
1091 41684912 : * vicen_init(n)/real(nilyr,kind=dbl_kind)
1092 47842864 : fhocn = fhocn + dfhocn
1093 : enddo ! nilyr
1094 :
1095 16990104 : do k = 1, nslyr
1096 : ! heat flux to coupler for snow melt (dfhocn < 0)
1097 : dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
1098 10832152 : * vsnon_init(n)/real(nslyr,kind=dbl_kind)
1099 16990104 : fhocn = fhocn + dfhocn
1100 : enddo ! nslyr
1101 :
1102 6157952 : if (tr_aero) then
1103 443840 : do k = 1, n_aero
1104 : faero_ocn(k) = faero_ocn(k) &
1105 : + (vsnon_init(n) * (trcrn(nt_aero +4*(k-1),n) &
1106 : + trcrn(nt_aero+1+4*(k-1),n)) &
1107 : + vicen_init(n) * (trcrn(nt_aero+2+4*(k-1),n) &
1108 : + trcrn(nt_aero+3+4*(k-1),n))) &
1109 443840 : * rsiden(n) / dt
1110 : enddo ! k
1111 : endif ! tr_aero
1112 :
1113 6157952 : if (tr_iso) then
1114 537000 : do k = 1, n_iso
1115 : fiso_ocn(k) = fiso_ocn(k) &
1116 : + (vsnon_init(n)*trcrn(nt_isosno+k-1,n) &
1117 : + vicen_init(n)*trcrn(nt_isoice+k-1,n)) &
1118 537000 : * rsiden(n) / dt
1119 : enddo ! k
1120 : endif ! tr_iso
1121 :
1122 : !-----------------------------------------------------------------
1123 : ! Biogeochemistry
1124 : !-----------------------------------------------------------------
1125 :
1126 7390024 : if (z_tracers) then ! snow tracers
1127 709290 : dvssl = p5*vsnon_init(n)/real(nslyr,kind=dbl_kind) ! snow surface layer
1128 709290 : dvint = vsnon_init(n) - dvssl ! snow interior
1129 12055470 : do k = 1, nbtrcr
1130 : flux_bio(k) = flux_bio(k) &
1131 : + (trcrn(bio_index(k)+nblyr+1,n)*dvssl &
1132 : + trcrn(bio_index(k)+nblyr+2,n)*dvint) &
1133 12055470 : * rsiden(n) / dt
1134 : enddo
1135 : endif
1136 :
1137 : enddo ! n
1138 :
1139 1232072 : if (z_tracers) &
1140 : call lateral_melt_bgc(dt, &
1141 : rsiden, vicen_init, &
1142 141858 : trcrn, flux_bio)
1143 1232072 : if (icepack_warnings_aborted(subname)) return
1144 :
1145 :
1146 1232072 : if (tr_fsd) then
1147 :
1148 3002450 : trcrn(nt_fsd:nt_fsd+nfsd-1,:) = afsdn
1149 :
1150 60510 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
1151 60510 : if (icepack_warnings_aborted(subname)) return
1152 :
1153 : ! diagnostics
1154 588388 : do k = 1, nfsd
1155 527878 : d_afsd_latm(k) = c0
1156 3227778 : do n = 1, ncat
1157 : d_afsd_latm(k) = d_afsd_latm(k) &
1158 3167268 : + afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n)
1159 : end do
1160 : end do
1161 :
1162 60510 : deallocate(afsdn)
1163 60510 : deallocate(afsdn_init)
1164 60510 : deallocate(df_flx)
1165 60510 : deallocate(afsd_tmp)
1166 60510 : deallocate(d_afsd_tmp)
1167 60510 : deallocate(f_flx)
1168 :
1169 : end if
1170 :
1171 1972116 : end subroutine lateral_melt
1172 :
1173 : !=======================================================================
1174 : !
1175 : ! Given the volume of new ice grown in open water, compute its area
1176 : ! and thickness and add it to the appropriate category or categories.
1177 : !
1178 : ! NOTE: Usually all the new ice is added to category 1. An exception is
1179 : ! made if there is no open water or if the new ice is too thick
1180 : ! for category 1, in which case ice is distributed evenly over the
1181 : ! entire cell. Subroutine rebin should be called in case the ice
1182 : ! thickness lies outside category bounds after new ice formation.
1183 : !
1184 : ! When ice must be added to categories above category 1, the mushy
1185 : ! formulation (ktherm=2) adds it only to the bottom of the ice. When
1186 : ! added to only category 1, all formulations combine the new ice and
1187 : ! existing ice tracers as bulk quantities.
1188 : !
1189 : ! authors: William H. Lipscomb, LANL
1190 : ! Elizabeth C. Hunke, LANL
1191 : ! Adrian Turner, LANL
1192 : ! Lettie Roach, NIWA/VUW
1193 : !
1194 1972116 : subroutine add_new_ice (dt, &
1195 1972116 : hin_max, ktherm, &
1196 1972116 : aicen, trcrn, &
1197 1972116 : vicen, &
1198 : aice0, aice, &
1199 : frzmlt, frazil, &
1200 : frz_onset, yday, &
1201 : fresh, fsalt, &
1202 : Tf, sss, &
1203 1972116 : salinz, phi_init, &
1204 : dSin0_frazil, &
1205 0 : flux_bio, &
1206 1972116 : ocean_bio, &
1207 : frazil_diag, &
1208 1972116 : fiso_ocn, &
1209 : HDO_ocn, H2_16O_ocn, &
1210 : H2_18O_ocn, &
1211 : wave_sig_ht, &
1212 1972116 : wave_spectrum, &
1213 1972116 : wavefreq, &
1214 1972116 : dwavefreq, &
1215 1972116 : d_afsd_latg, &
1216 1972116 : d_afsd_newi)
1217 :
1218 : use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice
1219 :
1220 : integer (kind=int_kind), intent(in) :: &
1221 : ktherm ! type of thermodynamics (-1 none, 1 BL99, 2 mushy)
1222 :
1223 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
1224 : hin_max ! category boundaries (m)
1225 :
1226 : real (kind=dbl_kind), intent(in) :: &
1227 : dt , & ! time step (s)
1228 : aice , & ! total concentration of ice
1229 : frzmlt, & ! freezing/melting potential (W/m^2)
1230 : Tf , & ! freezing temperature (C)
1231 : sss ! sea surface salinity (ppt)
1232 :
1233 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1234 : aicen , & ! concentration of ice
1235 : vicen ! volume per unit area of ice (m)
1236 :
1237 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
1238 : trcrn ! ice tracers
1239 : ! 1: surface temperature
1240 :
1241 : real (kind=dbl_kind), intent(inout) :: &
1242 : aice0 , & ! concentration of open water
1243 : frazil , & ! frazil ice growth (m/step-->cm/day)
1244 : frazil_diag,& ! frazil ice growth diagnostic (m/step-->cm/day)
1245 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
1246 : fsalt ! salt flux to ocean (kg/m^2/s)
1247 :
1248 : real (kind=dbl_kind), intent(inout), optional :: &
1249 : frz_onset ! day of year that freezing begins (congel or frazil)
1250 :
1251 : real (kind=dbl_kind), intent(in), optional :: &
1252 : yday ! day of year
1253 :
1254 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1255 : salinz ! initial salinity profile
1256 :
1257 : real (kind=dbl_kind), intent(in) :: &
1258 : phi_init , & ! initial frazil liquid fraction
1259 : dSin0_frazil ! initial frazil bulk salinity reduction from sss
1260 :
1261 : ! BGC
1262 :
1263 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1264 : flux_bio ! tracer flux to ocean from biology (mmol/m^2/s)
1265 :
1266 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1267 : ocean_bio ! ocean concentration of biological tracer
1268 :
1269 : ! water isotopes
1270 :
1271 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
1272 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
1273 :
1274 : real (kind=dbl_kind), intent(in), optional :: &
1275 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
1276 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
1277 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
1278 :
1279 : ! floe size distribution
1280 : real (kind=dbl_kind), intent(in), optional :: &
1281 : wave_sig_ht ! significant height of waves globally (m)
1282 :
1283 : real (kind=dbl_kind), dimension(:), intent(in), optional :: &
1284 : wave_spectrum ! ocean surface wave spectrum, E(f) (m^2 s)
1285 :
1286 : real(kind=dbl_kind), dimension(:), intent(in), optional :: &
1287 : wavefreq, & ! wave frequencies (s^-1)
1288 : dwavefreq ! wave frequency bin widths (s^-1)
1289 :
1290 : real (kind=dbl_kind), dimension(:), intent(out), optional :: &
1291 : ! change in thickness distribution (area)
1292 : d_afsd_latg , & ! due to fsd lateral growth
1293 : d_afsd_newi ! new ice formation
1294 :
1295 : ! local variables
1296 :
1297 : integer (kind=int_kind) :: &
1298 : ncats , & ! max categories to which new ice is added, initially
1299 : n , & ! ice category index
1300 : k , & ! ice layer index
1301 : it ! aerosol tracer index
1302 :
1303 : real (kind=dbl_kind) :: &
1304 : ai0new , & ! area of new ice added to cat 1
1305 : vi0new , & ! volume of new ice added to cat 1
1306 : hsurp , & ! thickness of new ice added to each cat
1307 : fnew , & ! heat flx to open water for new ice (W/m^2)
1308 : hi0new , & ! thickness of new ice
1309 : hi0max , & ! max allowed thickness of new ice
1310 : vsurp , & ! volume of new ice added to each cat
1311 : vtmp , & ! total volume of new and old ice
1312 : area1 , & ! starting fractional area of existing ice
1313 : alvl , & ! starting level ice area
1314 : dfresh , & ! change in fresh
1315 : dfsalt , & ! change in fsalt
1316 : vi0tmp , & ! frzmlt part of frazil
1317 : Ti , & ! frazil temperature
1318 : qi0new , & ! frazil ice enthalpy
1319 : Si0new , & ! frazil ice bulk salinity
1320 : vi0_init , & ! volume of new ice
1321 : vice1 , & ! starting volume of existing ice
1322 : vice_init, vice_final, & ! ice volume summed over categories
1323 : eice_init, eice_final ! ice energy summed over categories
1324 :
1325 : real (kind=dbl_kind) :: frazil_conc
1326 :
1327 : real (kind=dbl_kind), dimension (nilyr) :: &
1328 3944232 : Sprofile ! salinity profile used for new ice additions
1329 :
1330 : character (len=char_len) :: &
1331 : fieldid ! field identifier
1332 :
1333 : real (kind=dbl_kind), dimension (ncat) :: &
1334 3944232 : eicen, & ! energy of melting for each ice layer (J/m^2)
1335 3944232 : aicen_init, & ! fractional area of ice
1336 3944232 : vicen_init ! volume per unit area of ice (m)
1337 :
1338 : ! floe size distribution
1339 : real (kind=dbl_kind), dimension (:,:), allocatable :: &
1340 1972116 : afsdn ! floe size distribution tracer (originally areal_mfstd_init)
1341 :
1342 : ! real (kind=dbl_kind), dimension (nfsd) :: &
1343 : ! afsd , & ! fsd tracer for each thickness category
1344 :
1345 : real (kind=dbl_kind), dimension(ncat) :: & ! for now
1346 : ! change in thickness distribution (area)
1347 3944232 : d_an_latg , & ! due to fsd lateral growth
1348 3944232 : d_an_newi ! new ice formation
1349 :
1350 : real (kind=dbl_kind), dimension (ncat) :: &
1351 3944232 : d_an_tot, & ! change in the ITD due to lateral growth and new ice
1352 3944232 : area2 ! area after lateral growth and before new ice formation
1353 :
1354 : real (kind=dbl_kind), dimension (ncat) :: &
1355 3944232 : vin0new ! volume of new ice added to any thickness cat
1356 :
1357 : real (kind=dbl_kind) :: &
1358 : latsurf_area, & ! fractional area of ice on sides of floes
1359 : lead_area , & ! fractional area of ice in lead region
1360 : G_radial , & ! lateral melt rate (m/s)
1361 : tot_latg , & ! total fsd lateral growth in open water
1362 : ai0mod ! ai0new - tot_latg
1363 :
1364 : character(len=*),parameter :: subname='(add_new_ice)'
1365 :
1366 : !-----------------------------------------------------------------
1367 : ! initialize
1368 : !-----------------------------------------------------------------
1369 :
1370 1972116 : hsurp = c0
1371 1972116 : hi0new = c0
1372 1972116 : ai0new = c0
1373 11543112 : d_an_latg(:) = c0
1374 11543112 : d_an_tot(:) = c0
1375 11543112 : d_an_newi(:) = c0
1376 11543112 : vin0new(:) = c0
1377 :
1378 1972116 : if (tr_fsd) then
1379 993708 : d_afsd_latg(:) = c0 ! diagnostics
1380 993708 : d_afsd_newi(:) = c0
1381 : end if
1382 :
1383 11543112 : area2(:) = aicen(:)
1384 1972116 : lead_area = c0
1385 1972116 : latsurf_area = c0
1386 1972116 : G_radial = c0
1387 1972116 : tot_latg = c0
1388 1972116 : if (ncat > 1) then
1389 1899720 : hi0max = hin_max(1)*0.9_dbl_kind ! not too close to boundary
1390 : else
1391 72396 : hi0max = bignum ! big number
1392 : endif
1393 :
1394 1972116 : if (tr_fsd) then
1395 98676 : allocate(afsdn(nfsd,ncat))
1396 5067216 : afsdn(:,:) = c0
1397 98676 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:))
1398 98676 : if (icepack_warnings_aborted(subname)) return
1399 : endif
1400 :
1401 11543112 : do n = 1, ncat
1402 9570996 : aicen_init(n) = aicen(n)
1403 9570996 : vicen_init(n) = vicen(n)
1404 9570996 : eicen(n) = c0
1405 11543112 : if (tr_fsd) then
1406 4968540 : do k = 1, nfsd
1407 4968540 : afsdn(k,n) = trcrn(nt_fsd+k-1,n)
1408 : enddo
1409 : endif
1410 : enddo
1411 :
1412 1972116 : if (conserv_check) then
1413 :
1414 434376 : do n = 1, ncat
1415 2968236 : do k = 1, nilyr
1416 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
1417 2895840 : * vicen(n)/real(nilyr,kind=dbl_kind)
1418 : enddo
1419 : enddo
1420 72396 : call column_sum (ncat, vicen, vice_init)
1421 72396 : if (icepack_warnings_aborted(subname)) return
1422 72396 : call column_sum (ncat, eicen, eice_init)
1423 72396 : if (icepack_warnings_aborted(subname)) return
1424 :
1425 : endif ! conserv_check
1426 :
1427 : !-----------------------------------------------------------------
1428 : ! Compute average enthalpy of new ice.
1429 : ! Sprofile is the salinity profile used when adding new ice to
1430 : ! all categories, for ktherm/=2, and to category 1 for all ktherm.
1431 : !
1432 : ! NOTE: POP assumes new ice is fresh!
1433 : !-----------------------------------------------------------------
1434 :
1435 1972116 : if (ktherm == 2) then ! mushy
1436 1537740 : if (sss > c2 * dSin0_frazil) then
1437 1537740 : Si0new = sss - dSin0_frazil
1438 : else
1439 0 : Si0new = sss**2 / (c4*dSin0_frazil)
1440 : endif
1441 12301920 : do k = 1, nilyr
1442 12301920 : Sprofile(k) = Si0new
1443 : enddo
1444 1537740 : Ti = min(liquidus_temperature_mush(Si0new/phi_init), Tliquidus_max)
1445 1537740 : qi0new = icepack_enthalpy_mush(Ti, Si0new)
1446 : else
1447 2606256 : do k = 1, nilyr
1448 2606256 : Sprofile(k) = salinz(k)
1449 : enddo
1450 434376 : qi0new = -rhoi*Lfresh
1451 : endif ! ktherm
1452 :
1453 : !-----------------------------------------------------------------
1454 : ! Compute the volume, area, and thickness of new ice.
1455 : !-----------------------------------------------------------------
1456 :
1457 1972116 : fnew = max (frzmlt, c0) ! fnew > 0 iff frzmlt > 0
1458 1972116 : vi0new = -fnew*dt / qi0new ! note sign convention, qi < 0
1459 1972116 : vi0_init = vi0new ! for bgc
1460 :
1461 : ! increment ice volume and energy
1462 1972116 : if (conserv_check) then
1463 72396 : vice_init = vice_init + vi0new
1464 72396 : eice_init = eice_init + vi0new*qi0new
1465 : endif
1466 :
1467 : ! history diagnostics
1468 1972116 : frazil = vi0new
1469 :
1470 1972116 : if (present(frz_onset) .and. present(yday)) then
1471 1972116 : if (frazil > puny .and. frz_onset < puny) frz_onset = yday
1472 : endif
1473 :
1474 : !-----------------------------------------------------------------
1475 : ! Update fresh water and salt fluxes.
1476 : !
1477 : ! NOTE: POP assumes fresh water and salt flux due to frzmlt > 0
1478 : ! is NOT included in fluxes fresh and fsalt.
1479 : !-----------------------------------------------------------------
1480 :
1481 1972116 : dfresh = c0
1482 1972116 : dfsalt = c0
1483 1972116 : if (cpl_frazil == 'external') then
1484 : ! do nothing here, calculations are in the coupler or elsewhere
1485 : else
1486 1972116 : if (update_ocn_f) then
1487 72396 : dfresh = -rhoi*vi0new/dt
1488 1899720 : elseif (cpl_frazil == 'fresh_ice_correction' .and. ktherm == 2) then
1489 : ! correct frazil fluxes for mushy
1490 1537740 : vi0tmp = fnew*dt / (rhoi*Lfresh) ! ocn/cpl assumes frazil volume is pure, fresh ice
1491 1537740 : dfresh = -rhoi*(vi0new - vi0tmp)/dt
1492 1537740 : frazil_diag = frazil - vi0tmp
1493 : ! else
1494 : ! do nothing - other correction options could be implemented in the future
1495 : endif
1496 :
1497 1972116 : if (saltflux_option == 'prognostic') then
1498 46116 : dfsalt = Si0new*p001*dfresh
1499 : else
1500 1926000 : dfsalt = ice_ref_salinity*p001*dfresh
1501 : endif
1502 1972116 : fresh = fresh + dfresh
1503 1972116 : fsalt = fsalt + dfsalt
1504 : endif
1505 :
1506 : !-----------------------------------------------------------------
1507 : ! Decide how to distribute the new ice.
1508 : !-----------------------------------------------------------------
1509 :
1510 1972116 : if (vi0new > c0) then
1511 :
1512 703244 : if (tr_fsd) then ! lateral growth of existing ice
1513 : ! calculate change in conc due to lateral growth
1514 : ! update vi0new, without change to afsdn or aicen
1515 : call fsd_lateral_growth(dt, aice, &
1516 : aicen, vicen, &
1517 : vi0new, frazil, &
1518 : afsdn, &
1519 : lead_area, latsurf_area, &
1520 : G_radial, d_an_latg, &
1521 38145 : tot_latg)
1522 38145 : if (icepack_warnings_aborted(subname)) return
1523 :
1524 : ! volume added to each category from lateral growth
1525 228870 : do n = 1, ncat
1526 228870 : if (aicen(n) > puny) then
1527 190692 : vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
1528 : endif
1529 : end do
1530 :
1531 : ! check lateral growth doesn't exceed total growth
1532 : ! if it does, adjust it
1533 228870 : if (SUM(vin0new)>vi0new) then
1534 0 : write(warnstr,*) subname,'lateral growth warning ',vi0new,SUM(vin0new)
1535 0 : call icepack_warnings_add(warnstr)
1536 0 : vin0new(:) = vin0new(:)*vi0new/SUM(vin0new)
1537 : end if
1538 :
1539 228870 : vi0new = vi0new - SUM(vin0new)
1540 228870 : frazil = frazil - SUM(vin0new)
1541 :
1542 : endif
1543 :
1544 703244 : ai0mod = aice0
1545 : ! separate frazil ice growth from lateral ice growth
1546 703244 : if (tr_fsd) ai0mod = aice0-tot_latg
1547 :
1548 : ! new ice area and thickness
1549 : ! hin_max(0) < new ice thickness < hin_max(1)
1550 703244 : if (ai0mod > puny) then
1551 703214 : hi0new = max(vi0new/ai0mod, hfrazilmin)
1552 703214 : if (hi0new > hi0max .and. ai0mod+puny < c1) then
1553 : ! distribute excess volume over all categories (below)
1554 32838 : hi0new = hi0max
1555 32838 : ai0new = ai0mod
1556 32838 : vsurp = vi0new - ai0new*hi0new
1557 32838 : hsurp = vsurp / aice
1558 32838 : vi0new = ai0new*hi0new
1559 : else
1560 : ! put ice in a single category, with hsurp = 0
1561 670376 : ai0new = vi0new/hi0new
1562 : endif
1563 : else ! aice0 < puny
1564 30 : hsurp = vi0new/aice ! new thickness in each cat
1565 30 : vi0new = c0
1566 : endif ! aice0 > puny
1567 :
1568 : ! combine areal change from new ice growth and lateral growth
1569 703244 : d_an_newi(1) = ai0new
1570 3229076 : d_an_tot(2:ncat) = d_an_latg(2:ncat)
1571 703244 : d_an_tot(1) = d_an_latg(1) + d_an_newi(1)
1572 703244 : if (tr_fsd) then
1573 38145 : vin0new(1) = vin0new(1) + ai0new*hi0new ! not BFB
1574 : else
1575 665099 : vin0new(1) = vi0new
1576 : endif
1577 :
1578 : endif ! vi0new > puny
1579 :
1580 : !-----------------------------------------------------------------
1581 : ! Distribute excess ice volume among ice categories by increasing
1582 : ! ice thickness, leaving ice area unchanged.
1583 : !
1584 : ! NOTE: If new ice contains globally conserved tracers
1585 : ! (e.g., isotopes from seawater), code must be added here.
1586 : !
1587 : ! The mushy formulation (ktherm=2) puts the new ice only at the
1588 : ! bottom of existing ice and adjusts the layers accordingly.
1589 : ! The other formulations distribute the new ice throughout the
1590 : ! existing ice column.
1591 : !-----------------------------------------------------------------
1592 :
1593 1972116 : if (hsurp > c0) then ! add ice to all categories
1594 :
1595 197208 : do n = 1, ncat
1596 :
1597 164340 : vsurp = hsurp * aicen(n)
1598 :
1599 : ! update ice age due to freezing (new ice age = dt)
1600 164340 : vtmp = vicen(n) + vsurp
1601 164340 : if (tr_iage .and. vtmp > puny) &
1602 : trcrn(nt_iage,n) = &
1603 3925 : (trcrn(nt_iage,n)*vicen(n) + dt*vsurp) / vtmp
1604 :
1605 164340 : if (tr_lvl .and. vicen(n) > puny) then
1606 : trcrn(nt_vlvl,n) = &
1607 : (trcrn(nt_vlvl,n)*vicen(n) + &
1608 160303 : trcrn(nt_alvl,n)*vsurp) / vtmp
1609 : endif
1610 :
1611 164340 : if (tr_aero .and. vtmp > puny) then
1612 18990 : do it = 1, n_aero
1613 : trcrn(nt_aero+2+4*(it-1),n) = &
1614 9495 : trcrn(nt_aero+2+4*(it-1),n)*vicen(n) / vtmp
1615 : trcrn(nt_aero+3+4*(it-1),n) = &
1616 18990 : trcrn(nt_aero+3+4*(it-1),n)*vicen(n) / vtmp
1617 : enddo
1618 : endif
1619 :
1620 164340 : if (tr_iso .and. vtmp > puny) then
1621 15700 : do it=1,n_iso
1622 11775 : frazil_conc = c0
1623 11775 : if (it==1) &
1624 3925 : frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
1625 11775 : if (it==2) &
1626 3925 : frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
1627 11775 : if (it==3) &
1628 3925 : frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
1629 :
1630 : ! dilution and uptake in the ice
1631 : trcrn(nt_isoice+it-1,n) &
1632 : = (trcrn(nt_isoice+it-1,n)*vicen(n) &
1633 : + frazil_conc*rhoi*vsurp) &
1634 11775 : / vtmp
1635 :
1636 : fiso_ocn(it) = fiso_ocn(it) &
1637 15700 : - frazil_conc*rhoi*vsurp/dt
1638 : enddo
1639 : endif
1640 :
1641 : ! update category volumes
1642 164340 : vicen(n) = vtmp
1643 :
1644 197208 : if (ktherm == 2) then
1645 164340 : vsurp = hsurp * aicen(n) ! note - save this above?
1646 164340 : vtmp = vicen(n) - vsurp ! vicen is the new volume
1647 164340 : if (vicen(n) > c0) then
1648 : call update_vertical_tracers( &
1649 : trcrn(nt_qice:nt_qice+nilyr-1,n), &
1650 164228 : vtmp, vicen(n), qi0new)
1651 164228 : if (icepack_warnings_aborted(subname)) return
1652 : call update_vertical_tracers( &
1653 : trcrn(nt_sice:nt_sice+nilyr-1,n), &
1654 164228 : vtmp, vicen(n), Si0new)
1655 164228 : if (icepack_warnings_aborted(subname)) return
1656 : endif
1657 : else
1658 0 : do k = 1, nilyr
1659 : ! factor of nilyr cancels out
1660 0 : vsurp = hsurp * aicen(n) ! note - save this above?
1661 0 : vtmp = vicen(n) - vsurp ! vicen is the new volume
1662 0 : if (vicen(n) > c0) then
1663 : ! enthalpy
1664 : trcrn(nt_qice+k-1,n) = &
1665 0 : (trcrn(nt_qice+k-1,n)*vtmp + qi0new*vsurp) / vicen(n)
1666 : ! salinity
1667 : trcrn(nt_sice+k-1,n) = &
1668 0 : (trcrn(nt_sice+k-1,n)*vtmp + Sprofile(k)*vsurp) / vicen(n)
1669 : endif
1670 : enddo ! k
1671 : endif ! ktherm
1672 :
1673 : enddo ! n
1674 :
1675 : endif ! hsurp > 0
1676 :
1677 : !-----------------------------------------------------------------
1678 : ! Combine new ice grown in open water with ice categories.
1679 : ! Using the floe size distribution, ice is added laterally to all
1680 : ! categories; otherwise it is added to category 1.
1681 : ! Assume that vsnon and esnon are unchanged.
1682 : ! The mushy formulation assumes salt from frazil is added uniformly
1683 : ! to category 1, while the others use a salinity profile.
1684 : !-----------------------------------------------------------------
1685 :
1686 1972116 : ncats = 1 ! add new ice to category 1 by default
1687 1972116 : if (tr_fsd) ncats = ncat ! add new ice laterally to all categories
1688 :
1689 4338936 : do n = 1, ncats
1690 :
1691 4338936 : if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then ! add ice to category n
1692 :
1693 854047 : area1 = aicen(n) ! save
1694 854047 : vice1 = vicen(n) ! save
1695 854047 : area2(n) = aicen_init(n) + d_an_latg(n) ! save area after latg, before newi
1696 854047 : aicen(n) = aicen(n) + d_an_tot(n) ! after lateral growth and new ice growth
1697 :
1698 854047 : aice0 = aice0 - d_an_tot(n)
1699 854047 : vicen(n) = vicen(n) + vin0new(n)
1700 :
1701 854047 : trcrn(nt_Tsfc,n) = (trcrn(nt_Tsfc,n)*area1 + Tf*d_an_tot(n))/aicen(n)
1702 854047 : trcrn(nt_Tsfc,n) = min (trcrn(nt_Tsfc,n), c0)
1703 :
1704 854047 : if (tr_FY) then
1705 19997 : trcrn(nt_FY,n) = (trcrn(nt_FY,n)*area1 + d_an_tot(n))/aicen(n)
1706 19997 : trcrn(nt_FY,n) = min(trcrn(nt_FY,n), c1)
1707 : endif
1708 :
1709 854047 : if (tr_fsd) then ! evolve the floe size distribution
1710 : ! both new frazil ice and lateral growth
1711 : call fsd_add_new_ice (n, &
1712 : dt, ai0new, &
1713 : d_an_latg, d_an_newi, &
1714 : G_radial, area2, &
1715 : wave_sig_ht, &
1716 : wave_spectrum, &
1717 : wavefreq, &
1718 : dwavefreq, &
1719 : d_afsd_latg, &
1720 : d_afsd_newi, &
1721 : afsdn, aicen_init, &
1722 188978 : aicen, trcrn)
1723 188978 : if (icepack_warnings_aborted(subname)) return
1724 : endif
1725 :
1726 854047 : if (vicen(n) > puny) then
1727 854047 : if (tr_iage) &
1728 19259 : trcrn(nt_iage,n) = (trcrn(nt_iage,n)*vice1 + dt*vin0new(n))/vicen(n)
1729 :
1730 854047 : if (tr_aero) then
1731 55996 : do it = 1, n_aero
1732 : trcrn(nt_aero+2+4*(it-1),n) = &
1733 27998 : trcrn(nt_aero+2+4*(it-1),n)*vice1/vicen(n)
1734 : trcrn(nt_aero+3+4*(it-1),n) = &
1735 55996 : trcrn(nt_aero+3+4*(it-1),n)*vice1/vicen(n)
1736 : enddo
1737 : endif
1738 :
1739 854047 : if (tr_iso) then
1740 77036 : do it=1,n_iso
1741 57777 : frazil_conc = c0
1742 57777 : if (it==1) &
1743 19259 : frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
1744 57777 : if (it==2) &
1745 19259 : frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
1746 57777 : if (it==3) &
1747 19259 : frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
1748 :
1749 : trcrn(nt_isoice+it-1,1) &
1750 : = (trcrn(nt_isoice+it-1,1)*vice1) &
1751 57777 : + frazil_conc*rhoi*vi0new/vicen(1)
1752 :
1753 : fiso_ocn(it) = fiso_ocn(it) &
1754 77036 : - frazil_conc*rhoi*vi0new/dt
1755 : enddo
1756 : endif ! if iso
1757 :
1758 854047 : if (tr_lvl) then
1759 835406 : alvl = trcrn(nt_alvl,n)
1760 : trcrn(nt_alvl,n) = &
1761 835406 : (trcrn(nt_alvl,n)*area1 + d_an_tot(n))/aicen(n)
1762 : trcrn(nt_vlvl,n) = &
1763 835406 : (trcrn(nt_vlvl,n)*vice1 + vin0new(n))/vicen(n)
1764 : endif
1765 :
1766 854047 : if (tr_pond_topo) then
1767 : trcrn(nt_apnd,n) = &
1768 18641 : trcrn(nt_apnd,n)*area1/aicen(n)
1769 835406 : elseif (tr_pond_lvl) then
1770 712394 : if (trcrn(nt_alvl,n) > puny) then
1771 : trcrn(nt_apnd,n) = &
1772 712394 : trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n))
1773 : endif
1774 : endif
1775 : endif ! vicen > 0
1776 :
1777 6250748 : do k = 1, nilyr
1778 6250748 : if (vicen(n) > c0) then
1779 : ! factor of nilyr cancels out
1780 : ! enthalpy
1781 : trcrn(nt_qice+k-1,n) = &
1782 5396701 : (trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n)
1783 : ! salinity
1784 : trcrn(nt_sice+k-1,n) = &
1785 5396701 : (trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n)
1786 : endif
1787 : enddo ! nilyr
1788 : endif ! vin0new > c0
1789 : enddo ! n
1790 :
1791 1972116 : if (conserv_check) then
1792 :
1793 434376 : do n = 1, ncat
1794 361980 : eicen(n) = c0
1795 2968236 : do k = 1, nilyr
1796 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
1797 2895840 : * vicen(n)/real(nilyr,kind=dbl_kind)
1798 : enddo
1799 : enddo
1800 72396 : call column_sum (ncat, vicen, vice_final)
1801 72396 : if (icepack_warnings_aborted(subname)) return
1802 72396 : call column_sum (ncat, eicen, eice_final)
1803 72396 : if (icepack_warnings_aborted(subname)) return
1804 :
1805 72396 : fieldid = subname//':vice'
1806 : call column_conservation_check (fieldid, &
1807 : vice_init, vice_final, &
1808 72396 : puny)
1809 72396 : if (icepack_warnings_aborted(subname)) return
1810 72396 : fieldid = subname//':eice'
1811 : call column_conservation_check (fieldid, &
1812 : eice_init, eice_final, &
1813 72396 : puny*Lfresh*rhoi)
1814 72396 : if (icepack_warnings_aborted(subname)) return
1815 :
1816 : endif ! conserv_check
1817 :
1818 : !-----------------------------------------------------------------
1819 : ! Biogeochemistry
1820 : !-----------------------------------------------------------------
1821 1972116 : if (tr_brine .or. nbtrcr > 0) then
1822 : call add_new_ice_bgc(dt, ncats, &
1823 : aicen_init, vicen_init, vi0_init, &
1824 : aicen, vicen, vin0new, &
1825 : trcrn, &
1826 : ocean_bio, flux_bio, hsurp, &
1827 182052 : d_an_tot)
1828 182052 : if (icepack_warnings_aborted(subname)) return
1829 : endif
1830 :
1831 1972116 : if (tr_fsd) then
1832 98676 : deallocate(afsdn)
1833 : endif
1834 :
1835 1972116 : end subroutine add_new_ice
1836 :
1837 : !=======================================================================
1838 : !autodocument_start icepack_step_therm2
1839 : ! Driver for thermodynamic changes not needed for coupling:
1840 : ! transport in thickness space, lateral growth and melting.
1841 : !
1842 : ! authors: William H. Lipscomb, LANL
1843 : ! Elizabeth C. Hunke, LANL
1844 :
1845 3944232 : subroutine icepack_step_therm2(dt, hin_max, &
1846 0 : aicen, &
1847 1972116 : vicen, vsnon, &
1848 1972116 : aicen_init, vicen_init, &
1849 1972116 : trcrn, &
1850 : aice0, aice, &
1851 1972116 : trcr_depend, &
1852 1972116 : trcr_base, n_trcr_strata, &
1853 1972116 : nt_strata, &
1854 : Tf, sss, &
1855 1972116 : salinz, &
1856 0 : rsiden, meltl, &
1857 : wlat, &
1858 : frzmlt, frazil, &
1859 : frain, fpond, &
1860 : fresh, fsalt, &
1861 : fhocn, update_ocn_f, &
1862 1972116 : faero_ocn, &
1863 1972116 : first_ice, fzsal, &
1864 3944232 : flux_bio, ocean_bio, &
1865 : frazil_diag, &
1866 : frz_onset, yday, &
1867 1972116 : fiso_ocn, HDO_ocn, &
1868 : H2_16O_ocn, H2_18O_ocn, &
1869 : wave_sig_ht, &
1870 1972116 : wave_spectrum, &
1871 1972116 : wavefreq, &
1872 1972116 : dwavefreq, &
1873 1972116 : d_afsd_latg, d_afsd_newi, &
1874 1972116 : d_afsd_latm, d_afsd_weld)
1875 :
1876 : use icepack_parameters, only: icepack_init_parameters
1877 :
1878 : logical (kind=log_kind), intent(in), optional :: &
1879 : update_ocn_f ! if true, update fresh water and salt fluxes
1880 :
1881 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
1882 : hin_max ! category boundaries (m)
1883 :
1884 : real (kind=dbl_kind), intent(in) :: &
1885 : dt , & ! time step
1886 : Tf , & ! freezing temperature (C)
1887 : sss , & ! sea surface salinity (ppt)
1888 : frzmlt ! freezing/melting potential (W/m^2)
1889 :
1890 : integer (kind=int_kind), dimension (:), intent(in) :: &
1891 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
1892 : n_trcr_strata ! number of underlying tracer layers
1893 :
1894 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
1895 : trcr_base ! = 0 or 1 depending on tracer dependency
1896 : ! argument 2: (1) aice, (2) vice, (3) vsno
1897 :
1898 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
1899 : nt_strata ! indices of underlying tracer layers
1900 :
1901 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1902 : rsiden , & ! fraction of ice that melts laterally
1903 : salinz , & ! initial salinity profile
1904 : ocean_bio ! ocean concentration of biological tracer
1905 :
1906 : real (kind=dbl_kind), intent(inout) :: &
1907 : aice , & ! sea ice concentration
1908 : aice0 , & ! concentration of open water
1909 : frain , & ! rainfall rate (kg/m^2 s)
1910 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
1911 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
1912 : fsalt , & ! salt flux to ocean (kg/m^2/s)
1913 : fhocn , & ! net heat flux to ocean (W/m^2)
1914 : meltl , & ! lateral ice melt (m/step-->cm/day)
1915 : frazil , & ! frazil ice growth (m/step-->cm/day)
1916 : frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day)
1917 :
1918 : real (kind=dbl_kind), intent(inout), optional :: &
1919 : fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) (deprecated)
1920 :
1921 : real (kind=dbl_kind), intent(in), optional :: &
1922 : wlat ! lateral melt rate (m/s)
1923 :
1924 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1925 : aicen_init,& ! initial concentration of ice
1926 : vicen_init,& ! initial volume per unit area of ice (m)
1927 : aicen , & ! concentration of ice
1928 : vicen , & ! volume per unit area of ice (m)
1929 : vsnon , & ! volume per unit area of snow (m)
1930 : faero_ocn, & ! aerosol flux to ocean (kg/m^2/s)
1931 : flux_bio ! all bio fluxes to ocean
1932 :
1933 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
1934 : trcrn ! tracers
1935 :
1936 : logical (kind=log_kind), dimension(:), intent(inout) :: &
1937 : first_ice ! true until ice forms
1938 :
1939 : real (kind=dbl_kind), intent(inout), optional :: &
1940 : frz_onset ! day of year that freezing begins (congel or frazil)
1941 :
1942 : real (kind=dbl_kind), intent(in), optional :: &
1943 : yday ! day of year
1944 :
1945 : ! water isotopes
1946 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
1947 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
1948 :
1949 : real (kind=dbl_kind), intent(in), optional :: &
1950 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
1951 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
1952 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
1953 :
1954 : real (kind=dbl_kind), intent(in), optional :: &
1955 : wave_sig_ht ! significant height of waves in ice (m)
1956 :
1957 : real (kind=dbl_kind), dimension(:), intent(in), optional :: &
1958 : wave_spectrum ! ocean surface wave spectrum E(f) (m^2 s)
1959 :
1960 : real(kind=dbl_kind), dimension(:), intent(in), optional :: &
1961 : wavefreq, & ! wave frequencies (s^-1)
1962 : dwavefreq ! wave frequency bin widths (s^-1)
1963 :
1964 : real (kind=dbl_kind), dimension(:), intent(out), optional :: &
1965 : ! change in floe size distribution (area)
1966 : d_afsd_latg, & ! due to fsd lateral growth
1967 : d_afsd_newi, & ! new ice formation
1968 : d_afsd_latm, & ! lateral melt
1969 : d_afsd_weld ! welding
1970 :
1971 : !autodocument_end
1972 :
1973 : ! local variables
1974 :
1975 : logical (kind=log_kind), save :: &
1976 : first_call = .true. ! first call flag
1977 :
1978 : character(len=*),parameter :: subname='(icepack_step_therm2)'
1979 :
1980 : !-----------------------------------------------------------------
1981 : ! Check optional arguments and set local values
1982 : !-----------------------------------------------------------------
1983 :
1984 1972116 : if (present(update_ocn_f)) then
1985 0 : call icepack_init_parameters(update_ocn_f_in=update_ocn_f)
1986 : endif
1987 1972116 : if (icepack_chkoptargflag(first_call)) then
1988 83 : if (tr_iso) then
1989 2 : if (.not.(present(fiso_ocn) .and. &
1990 : present(HDO_ocn) .and. &
1991 : present(H2_16O_ocn) .and. &
1992 : present(H2_18O_ocn))) then
1993 0 : call icepack_warnings_add(subname//' error in iso arguments, tr_iso=T')
1994 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1995 0 : return
1996 : endif
1997 : endif
1998 83 : if (tr_fsd) then
1999 4 : if (.not.(present(wlat) .and. &
2000 : present(wave_sig_ht) .and. &
2001 : present(wave_spectrum) .and. &
2002 : present(wavefreq) .and. &
2003 : present(dwavefreq) .and. &
2004 : present(d_afsd_latg) .and. &
2005 : present(d_afsd_newi) .and. &
2006 : present(d_afsd_latm) .and. &
2007 : present(d_afsd_weld))) then
2008 0 : call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T')
2009 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
2010 0 : return
2011 : endif
2012 : endif
2013 : endif
2014 :
2015 : !-----------------------------------------------------------------
2016 : ! Let rain drain through to the ocean.
2017 : !-----------------------------------------------------------------
2018 :
2019 1972116 : fresh = fresh + frain * aice
2020 :
2021 : !-----------------------------------------------------------------
2022 : ! Given thermodynamic growth rates, transport ice between
2023 : ! thickness categories.
2024 : !-----------------------------------------------------------------
2025 :
2026 : ! call ice_timer_start(timer_catconv) ! category conversions
2027 :
2028 : !-----------------------------------------------------------------
2029 : ! Compute fractional ice area in each grid cell.
2030 : !-----------------------------------------------------------------
2031 :
2032 4563756 : flux_bio(:) = c0
2033 :
2034 1972116 : call aggregate_area (aicen, aice, aice0)
2035 1972116 : if (icepack_warnings_aborted(subname)) return
2036 :
2037 1972116 : if (kitd == 1) then
2038 :
2039 : !-----------------------------------------------------------------
2040 : ! Identify grid cells with ice.
2041 : !-----------------------------------------------------------------
2042 :
2043 1827324 : if (aice > puny) then
2044 :
2045 : call linear_itd (hin_max, &
2046 : trcr_depend, &
2047 : trcr_base, &
2048 : n_trcr_strata, &
2049 : nt_strata, &
2050 : aicen_init, &
2051 : vicen_init, &
2052 : aicen, &
2053 : trcrn, &
2054 : vicen, &
2055 : vsnon, &
2056 : aice , &
2057 : aice0 , &
2058 1790590 : fpond, Tf )
2059 1790590 : if (icepack_warnings_aborted(subname)) return
2060 :
2061 : endif ! aice > puny
2062 :
2063 : endif ! kitd = 1
2064 :
2065 : ! call ice_timer_stop(timer_catconv) ! category conversions
2066 :
2067 : !-----------------------------------------------------------------
2068 : ! Add frazil ice growing in leads.
2069 : !-----------------------------------------------------------------
2070 :
2071 : ! identify ice-ocean cells
2072 :
2073 : call add_new_ice (dt, &
2074 : hin_max, ktherm, &
2075 : aicen, trcrn, &
2076 : vicen, &
2077 : aice0, aice, &
2078 : frzmlt, frazil, &
2079 : frz_onset, yday, &
2080 : fresh, fsalt, &
2081 : Tf, sss, &
2082 : salinz, phi_init, &
2083 : dSin0_frazil, &
2084 : flux_bio, &
2085 : ocean_bio, &
2086 : frazil_diag, fiso_ocn, &
2087 : HDO_ocn, H2_16O_ocn, &
2088 : H2_18O_ocn, &
2089 : wave_sig_ht, &
2090 : wave_spectrum, &
2091 : wavefreq, dwavefreq, &
2092 1972116 : d_afsd_latg, d_afsd_newi)
2093 :
2094 1972116 : if (icepack_warnings_aborted(subname)) return
2095 :
2096 : !-----------------------------------------------------------------
2097 : ! Melt ice laterally.
2098 : !-----------------------------------------------------------------
2099 :
2100 : call lateral_melt (dt, fpond, &
2101 : fresh, fsalt, &
2102 : fhocn, faero_ocn, &
2103 : fiso_ocn, &
2104 : rsiden, meltl, &
2105 : wlat, &
2106 : aicen, vicen, &
2107 : vsnon, trcrn, &
2108 : flux_bio, &
2109 1972116 : d_afsd_latm)
2110 1972116 : if (icepack_warnings_aborted(subname)) return
2111 :
2112 : ! Floe welding during freezing conditions
2113 1972116 : if (tr_fsd) then
2114 : call fsd_weld_thermo (dt, frzmlt, &
2115 : aicen, trcrn, &
2116 98676 : d_afsd_weld)
2117 98676 : if (icepack_warnings_aborted(subname)) return
2118 : endif
2119 :
2120 : !-----------------------------------------------------------------
2121 : ! For the special case of a single category, adjust the area and
2122 : ! volume (assuming that half the volume change decreases the
2123 : ! thickness, and the other half decreases the area).
2124 : !-----------------------------------------------------------------
2125 :
2126 1972116 : if (ncat==1) &
2127 : call reduce_area (hin_max (0), &
2128 : aicen (1), vicen (1), &
2129 72396 : aicen_init(1), vicen_init(1))
2130 1972116 : if (icepack_warnings_aborted(subname)) return
2131 :
2132 : !-----------------------------------------------------------------
2133 : ! ITD cleanup: Rebin thickness categories if necessary, and remove
2134 : ! categories with very small areas.
2135 : !-----------------------------------------------------------------
2136 :
2137 : call cleanup_itd (dt, hin_max, &
2138 : aicen, trcrn(1:ntrcr,:), &
2139 : vicen, vsnon, &
2140 : aice0, aice, &
2141 : tr_aero, &
2142 : tr_pond_topo, &
2143 : first_ice, &
2144 : trcr_depend, trcr_base, &
2145 : n_trcr_strata, nt_strata, &
2146 : fpond, fresh, &
2147 : fsalt, fhocn, &
2148 : faero_ocn, fiso_ocn, &
2149 1972116 : flux_bio, Tf )
2150 1972116 : if (icepack_warnings_aborted(subname)) return
2151 :
2152 1972116 : first_call = .false.
2153 :
2154 : end subroutine icepack_step_therm2
2155 :
2156 : !=======================================================================
2157 :
2158 : end module icepack_therm_itd
2159 :
2160 : !=======================================================================
|
