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