Line data Source code
1 : !=========================================================================
2 : !
3 : ! This module contains the subroutines required to define
4 : ! a floe size distribution tracer for sea ice
5 : !
6 : ! Theory based on:
7 : !
8 : ! Horvat, C., & Tziperman, E. (2015). A prognostic model of the sea-ice
9 : ! floe size and thickness distribution. The Cryosphere, 9(6), 2119-2134.
10 : ! doi:10.5194/tc-9-2119-2015
11 : !
12 : ! and implementation described in:
13 : !
14 : ! Roach, L. A., Horvat, C., Dean, S. M., & Bitz, C. M. (2018). An emergent
15 : ! sea ice floe size distribution in a global coupled ocean--sea ice model.
16 : ! Journal of Geophysical Research: Oceans, 123(6), 4322-4337.
17 : ! doi:10.1029/2017JC013692
18 : !
19 : ! with some modifications.
20 : !
21 : ! For floe welding parameter and tensile mode parameter values, see
22 : !
23 : ! Roach, L. A., Smith, M. M., & Dean, S. M. (2018). Quantifying
24 : ! growth of pancake sea ice floes using images from drifting buoys.
25 : ! Journal of Geophysical Research: Oceans, 123(4), 2851-2866.
26 : ! doi: 10.1002/2017JC013693
27 : !
28 : ! Variable naming convention:
29 : ! for k = 1, nfsd and n = 1, ncat
30 : ! afsdn(k,n) = trcrn(:,:,nt_nfsd+k-1,n,:)
31 : ! afsd (k) is per thickness category or averaged over n
32 : ! afs is the associated scalar value for (k,n)
33 : !
34 : ! authors: Lettie Roach, VUW/NIWA
35 : ! C. M. Bitz, UW
36 : !
37 : ! 2016: CMB started
38 : ! 2016-8: LR worked on most of it
39 : ! 2019: ECH ported to Icepack
40 :
41 : !-----------------------------------------------------------------
42 :
43 : module icepack_fsd
44 :
45 : use icepack_kinds
46 : use icepack_parameters, only: c0, c1, c2, c4, p01, p1, p5, puny
47 : use icepack_parameters, only: pi, floeshape, wave_spec, bignum, gravit, rhoi
48 : use icepack_tracers, only: nt_fsd, tr_fsd, nfsd, ncat
49 : use icepack_warnings, only: warnstr, icepack_warnings_add
50 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
51 :
52 : implicit none
53 :
54 : private
55 : public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, &
56 : fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo, get_subdt_fsd
57 :
58 : real(kind=dbl_kind), dimension(:), allocatable, public :: &
59 : floe_rad_h, & ! fsd size higher bound in m (radius)
60 : floe_rad_c, & ! fsd size center in m (radius)
61 : floe_rad_l, & ! fsd size lower bound in m (radius)
62 : floe_binwidth, & ! fsd size binwidth in m (radius)
63 : floe_area_l, & ! fsd area at lower bound (m^2)
64 : floe_area_h, & ! fsd area at higher bound (m^2)
65 : floe_area_c, & ! fsd area at bin centre (m^2)
66 : floe_area_binwidth ! floe area bin width (m^2)
67 :
68 : character (len=35), dimension(:), allocatable, public :: &
69 : c_fsd_range ! string for history output
70 :
71 : integer(kind=int_kind), dimension(:,:), allocatable, public :: &
72 : iweld ! floe size categories that can combine
73 : ! during welding (dimensionless)
74 :
75 : !=======================================================================
76 :
77 : contains
78 :
79 : !=======================================================================
80 : ! Note that radius widths cannot be larger than twice previous
81 : ! category width or floe welding will not have an effect
82 : !
83 : ! Note also that the bound of the lowest floe size category is used
84 : ! to define the lead region width and the domain spacing for wave fracture
85 : !
86 : !autodocument_start icepack_init_fsd_bounds
87 : ! Initialize ice fsd bounds (call whether or not restarting)
88 : ! Define the bounds, midpoints and widths of floe size
89 : ! categories in area and radius
90 : !
91 : ! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
92 :
93 4 : subroutine icepack_init_fsd_bounds( &
94 4 : floe_rad_l_out, & ! fsd size lower bound in m (radius)
95 4 : floe_rad_c_out, & ! fsd size bin centre in m (radius)
96 4 : floe_binwidth_out, & ! fsd size bin width in m (radius)
97 4 : c_fsd_range_out, & ! string for history output
98 : write_diags) ! flag for writing diagnostics
99 :
100 : real(kind=dbl_kind), dimension(:), intent(out), optional :: &
101 : floe_rad_l_out, & ! fsd size lower bound in m (radius)
102 : floe_rad_c_out, & ! fsd size bin centre in m (radius)
103 : floe_binwidth_out ! fsd size bin width in m (radius)
104 :
105 : character (len=35), dimension(:), intent(out), optional :: &
106 : c_fsd_range_out ! string for history output
107 :
108 : logical (kind=log_kind), intent(in), optional :: &
109 : write_diags ! write diags flag
110 :
111 : !autodocument_end
112 :
113 : ! local variables
114 :
115 : integer (kind=int_kind) :: n, m, k
116 : integer (kind=int_kind) :: ierr
117 :
118 : real (kind=dbl_kind) :: test
119 :
120 : real (kind=dbl_kind), dimension(:), allocatable :: &
121 4 : lims, floe_rad
122 :
123 : character(len=8) :: c_fsd1,c_fsd2
124 : character(len=2) :: c_nf
125 : character(len=*), parameter :: subname='(icepack_init_fsd_bounds)'
126 :
127 4 : if (nfsd.eq.24) then
128 :
129 0 : allocate(lims(24+1))
130 :
131 : lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, &
132 : 5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, &
133 : 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, &
134 : 9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, &
135 : 3.35434988e+03, 4.55051413e+03, 6.17323164e+03, 8.37461170e+03, &
136 : 1.13610059e+04, 1.54123510e+04, 2.09084095e+04, 2.83643675e+04, &
137 0 : 3.84791270e+04 /)
138 :
139 4 : elseif (nfsd.eq.16) then
140 :
141 0 : allocate(lims(16+1))
142 :
143 : lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, &
144 : 5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, &
145 : 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, &
146 : 9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, &
147 0 : 3.35434988e+03 /)
148 :
149 4 : else if (nfsd.eq.12) then
150 :
151 3 : allocate(lims(12+1))
152 :
153 : lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, &
154 : 5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, &
155 : 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, &
156 45 : 9.45812834e+02 /)
157 :
158 1 : else if (nfsd.eq.1) then ! default case
159 :
160 1 : allocate(lims(1+1))
161 :
162 4 : lims = (/ 6.65000000e-02, 3.0e+02 /)
163 :
164 : else
165 :
166 : call icepack_warnings_add(subname//&
167 0 : ' floe size categories not defined for nfsd')
168 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
169 0 : return
170 :
171 : end if
172 :
173 : allocate( &
174 : floe_rad_h (nfsd), & ! fsd size higher bound in m (radius)
175 : floe_rad_l (nfsd), & ! fsd size lower bound in m (radius)
176 : floe_rad_c (nfsd), & ! fsd size center in m (radius)
177 : floe_rad (0:nfsd), & ! fsd bounds in m (radius)
178 : floe_area_l (nfsd), & ! fsd area at lower bound (m^2)
179 : floe_area_h (nfsd), & ! fsd area at higher bound (m^2)
180 : floe_area_c (nfsd), & ! fsd area at bin centre (m^2)
181 : floe_area_binwidth (nfsd), & ! floe area bin width (m^2)
182 : floe_binwidth (nfsd), & ! floe bin width (m)
183 : c_fsd_range (nfsd), & !
184 : iweld (nfsd, nfsd), & ! fsd categories that can weld
185 4 : stat=ierr)
186 4 : if (ierr/=0) then
187 0 : call icepack_warnings_add(subname//' Out of Memory fsd')
188 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
189 0 : return
190 : endif
191 :
192 45 : floe_rad_l = lims(1:nfsd )
193 45 : floe_rad_h = lims(2:nfsd+1)
194 45 : floe_rad_c = (floe_rad_h+floe_rad_l)/c2
195 :
196 45 : floe_area_l = c4*floeshape*floe_rad_l**2
197 45 : floe_area_h = c4*floeshape*floe_rad_h**2
198 : ! floe_area_c = c4*floeshape*floe_rad_c**2
199 : ! This is exactly in the middle of floe_area_h and floe_area_l
200 : ! Whereas the above calculation is closer to floe_area_l.
201 45 : floe_area_c = (floe_area_h+floe_area_l)/c2
202 :
203 45 : floe_binwidth = floe_rad_h - floe_rad_l
204 :
205 45 : floe_area_binwidth = floe_area_h - floe_area_l
206 :
207 : ! floe size categories that can combine during welding
208 474 : iweld(:,:) = -999
209 41 : do n = 1, nfsd
210 474 : do m = 1, nfsd
211 433 : test = floe_area_c(n) + floe_area_c(m)
212 5185 : do k = 1, nfsd-1
213 4752 : if ((test >= floe_area_l(k)) .and. (test < floe_area_h(k))) &
214 780 : iweld(n,m) = k
215 : end do
216 470 : if (test >= floe_area_l(nfsd)) iweld(n,m) = nfsd
217 : end do
218 : end do
219 :
220 4 : if (allocated(lims)) deallocate(lims)
221 :
222 : ! write fsd bounds
223 4 : floe_rad(0) = floe_rad_l(1)
224 41 : do n = 1, nfsd
225 37 : floe_rad(n) = floe_rad_h(n)
226 : ! Save character string to write to history file
227 37 : write (c_nf, '(i2)') n
228 37 : write (c_fsd1, '(f7.3)') floe_rad(n-1)
229 37 : write (c_fsd2, '(f7.3)') floe_rad(n)
230 41 : c_fsd_range(n)=c_fsd1//'m < fsd Cat '//c_nf//' < '//c_fsd2//'m'
231 : enddo
232 :
233 4 : if (present(write_diags)) then
234 4 : if (write_diags) then
235 4 : write(warnstr,*) ' '
236 4 : call icepack_warnings_add(warnstr)
237 4 : write(warnstr,*) subname
238 4 : call icepack_warnings_add(warnstr)
239 4 : write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
240 4 : call icepack_warnings_add(warnstr)
241 41 : do n = 1, nfsd
242 37 : write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
243 41 : call icepack_warnings_add(warnstr)
244 : enddo
245 4 : write(warnstr,*) ' '
246 4 : call icepack_warnings_add(warnstr)
247 : endif
248 : endif
249 :
250 4 : if (present(floe_rad_l_out)) then
251 0 : if (size(floe_rad_l_out) /= size(floe_rad_l)) then
252 0 : call icepack_warnings_add(subname//' floe_rad_l_out incorrect size')
253 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
254 0 : return
255 : endif
256 0 : floe_rad_l_out(:) = floe_rad_l(:)
257 : endif
258 :
259 4 : if (present(floe_rad_c_out)) then
260 4 : if (size(floe_rad_c_out) /= size(floe_rad_c)) then
261 0 : call icepack_warnings_add(subname//' floe_rad_c_out incorrect size')
262 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
263 0 : return
264 : endif
265 41 : floe_rad_c_out(:) = floe_rad_c(:)
266 : endif
267 :
268 4 : if (present(floe_binwidth_out)) then
269 0 : if (size(floe_binwidth_out) /= size(floe_binwidth)) then
270 0 : call icepack_warnings_add(subname//' floe_binwidth_out incorrect size')
271 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
272 0 : return
273 : endif
274 0 : floe_binwidth_out(:) = floe_binwidth(:)
275 : endif
276 :
277 4 : if (present(c_fsd_range_out)) then
278 0 : if (size(c_fsd_range_out) /= size(c_fsd_range)) then
279 0 : call icepack_warnings_add(subname//' c_fsd_range_out incorrect size')
280 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
281 0 : return
282 : endif
283 0 : c_fsd_range_out(:) = c_fsd_range(:)
284 : endif
285 :
286 4 : end subroutine icepack_init_fsd_bounds
287 :
288 : !=======================================================================
289 : ! When growing from no-ice conditions, initialize to zero.
290 : ! This allows the FSD to emerge, as described in Roach, Horvat et al. (2018)
291 : !
292 : ! Otherwise initalize with a power law, following Perovich
293 : ! & Jones (2014). The basin-wide applicability of such a
294 : ! prescribed power law has not yet been tested.
295 : !
296 : ! Perovich, D. K., & Jones, K. F. (2014). The seasonal evolution of
297 : ! sea ice floe size distribution. Journal of Geophysical Research: Oceans,
298 : ! 119(12), 8767–8777. doi:10.1002/2014JC010136
299 : !
300 : !autodocument_start icepack_init_fsd
301 : !
302 : ! Initialize the FSD
303 : !
304 : ! authors: Lettie Roach, NIWA/VUW
305 :
306 40 : subroutine icepack_init_fsd(ice_ic, afsd) ! floe size distribution tracer
307 :
308 : character(len=char_len_long), intent(in) :: &
309 : ice_ic ! method of ice cover initialization
310 :
311 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
312 : afsd ! floe size tracer: fraction distribution of floes
313 :
314 : !autodocument_end
315 :
316 : ! local variables
317 :
318 : real (kind=dbl_kind) :: alpha, totfrac
319 :
320 : integer (kind=int_kind) :: k
321 :
322 : real (kind=dbl_kind), dimension (nfsd) :: &
323 80 : num_fsd ! number distribution of floes
324 :
325 40 : if (trim(ice_ic) == 'none') then
326 :
327 0 : afsd(:) = c0
328 :
329 : else ! Perovich (2014)
330 :
331 : ! fraction of ice in each floe size and thickness category
332 : ! same for ALL cells (even where no ice) initially
333 40 : alpha = 2.1_dbl_kind
334 40 : totfrac = c0 ! total fraction of floes
335 410 : do k = 1, nfsd
336 370 : num_fsd(k) = (2*floe_rad_c(k))**(-alpha-c1) ! number distribution of floes
337 370 : afsd (k) = num_fsd(k)*floe_area_c(k)*floe_binwidth(k) ! fraction distribution of floes
338 410 : totfrac = totfrac + afsd(k)
339 : enddo
340 410 : afsd = afsd/totfrac ! normalize
341 :
342 : endif ! ice_ic
343 :
344 40 : end subroutine icepack_init_fsd
345 :
346 : !=======================================================================
347 : !autodocument_start icepack_cleanup_fsd
348 : !
349 : ! Clean up small values and renormalize
350 : !
351 : ! authors: Elizabeth Hunke, LANL
352 : !
353 219696 : subroutine icepack_cleanup_fsd (afsdn)
354 :
355 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
356 : afsdn ! floe size distribution tracer
357 :
358 : !autodocument_end
359 : ! local variables
360 :
361 : integer (kind=int_kind) :: &
362 : n ! thickness category index
363 :
364 : character(len=*), parameter :: subname='(icepack_cleanup_fsd)'
365 :
366 219696 : if (tr_fsd) then
367 :
368 1318176 : do n = 1, ncat
369 1098480 : call icepack_cleanup_fsdn(afsdn(:,n))
370 1318176 : if (icepack_warnings_aborted(subname)) return
371 : enddo
372 :
373 : endif ! tr_fsd
374 :
375 : end subroutine icepack_cleanup_fsd
376 :
377 : !=======================================================================
378 : !
379 : ! Clean up small values and renormalize -- per category
380 : !
381 : ! authors: Elizabeth Hunke, LANL
382 : !
383 :
384 1855021 : subroutine icepack_cleanup_fsdn (afsd)
385 :
386 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
387 : afsd ! floe size distribution tracer
388 :
389 : ! local variables
390 :
391 : integer (kind=int_kind) :: &
392 : k ! floe size category index
393 :
394 : real (kind=dbl_kind) :: &
395 : tot
396 :
397 19083774 : do k = 1, nfsd
398 19083774 : if (afsd(k) < puny) afsd(k) = c0
399 : enddo
400 :
401 19083774 : tot = sum(afsd(:))
402 1855021 : if (tot > puny) then
403 18523606 : do k = 1, nfsd
404 18523606 : afsd(k) = afsd(k) / tot ! normalize
405 : enddo
406 : else ! represents ice-free cell, so set to zero
407 560168 : afsd(:) = c0
408 : endif
409 :
410 1855021 : end subroutine icepack_cleanup_fsdn
411 :
412 : !=======================================================================
413 : !
414 : ! Given the joint ice thickness and floe size distribution, calculate
415 : ! the lead region and the total lateral surface area following Horvat
416 : ! and Tziperman (2015).
417 : !
418 : ! authors: Lettie Roach, NIWA/VUW
419 :
420 38145 : subroutine partition_area (aice, &
421 76290 : aicen, vicen, &
422 38145 : afsdn, lead_area, &
423 : latsurf_area)
424 :
425 : real (kind=dbl_kind), intent(in) :: &
426 : aice ! ice concentration
427 :
428 : real (kind=dbl_kind), dimension(:), intent(in) :: &
429 : aicen , & ! fractional area of ice
430 : vicen ! volume per unit area of ice (m)
431 :
432 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
433 : afsdn ! floe size distribution tracer
434 :
435 : real (kind=dbl_kind), intent(out) :: &
436 : lead_area , & ! the fractional area of the lead region
437 : latsurf_area ! the area covered by lateral surface of floes
438 :
439 : ! local variables
440 :
441 : integer (kind=int_kind) :: &
442 : n , & ! thickness category index
443 : k ! floe size index
444 :
445 : real (kind=dbl_kind) :: &
446 : width_leadreg, & ! width of lead region (m)
447 : thickness ! actual thickness of ice in thickness cat (m)
448 :
449 : !-----------------------------------------------------------------
450 : ! Initialize
451 : !-----------------------------------------------------------------
452 38145 : lead_area = c0
453 38145 : latsurf_area = c0
454 :
455 : ! Set the width of the lead region to be the smallest
456 : ! floe size category, as per Horvat & Tziperman (2015)
457 38145 : width_leadreg = floe_rad_c(1)
458 :
459 : ! Only calculate these areas where there is some ice
460 38145 : if (aice > puny) then
461 :
462 : ! lead area = sum of areas of annuli around all floes
463 228852 : do n = 1, ncat
464 2063622 : do k = 1, nfsd
465 : lead_area = lead_area + aicen(n) * afsdn(k,n) &
466 : * (c2*width_leadreg /floe_rad_c(k) &
467 2025480 : + width_leadreg**2/floe_rad_c(k)**2)
468 : enddo ! k
469 : enddo ! n
470 :
471 : ! cannot be greater than the open water fraction
472 38142 : lead_area=MIN(lead_area, c1-aice)
473 38142 : if (lead_area < c0) then
474 0 : if (lead_area < -puny) then
475 : ! stop 'lead_area lt0 in partition_area'
476 : else
477 0 : lead_area=MAX(lead_area,c0)
478 : end if
479 : end if
480 :
481 : ! Total fractional lateral surface area in each grid (per unit ocean area)
482 228852 : do n = 1, ncat
483 190710 : thickness = c0
484 190710 : if (aicen(n) > c0) thickness = vicen(n)/aicen(n)
485 :
486 2063622 : do k = 1, nfsd
487 : latsurf_area = latsurf_area &
488 2025480 : + afsdn(k,n) * aicen(n) * c2 * thickness/floe_rad_c(k)
489 : end do
490 : end do
491 :
492 : ! check
493 : ! if (latsurf_area < c0) stop 'negative latsurf_area'
494 : ! if (latsurf_area /= latsurf_area) stop 'latsurf_area NaN'
495 :
496 : end if ! aice
497 :
498 38145 : end subroutine partition_area
499 :
500 : !=======================================================================
501 : !
502 : ! Lateral growth at the edges of floes
503 : !
504 : ! Compute the portion of new ice growth that occurs at the edges of
505 : ! floes. The remainder will grow as new ice frazil ice in open water
506 : ! (away from floes).
507 : !
508 : ! See Horvat & Tziperman (2015) and Roach, Horvat et al. (2018).
509 : !
510 : ! authors: Lettie Roach, NIWA/VUW
511 : !
512 38145 : subroutine fsd_lateral_growth (dt, aice, &
513 38145 : aicen, vicen, &
514 : vi0new, &
515 : frazil, &
516 38145 : afsdn, &
517 : lead_area, latsurf_area, &
518 38145 : G_radial, d_an_latg, &
519 : tot_latg)
520 :
521 : real (kind=dbl_kind), intent(in) :: &
522 : dt , & ! time step (s)
523 : aice ! total concentration of ice
524 :
525 : real (kind=dbl_kind), dimension (:), intent(in) :: &
526 : aicen , & ! concentration of ice
527 : vicen ! volume per unit area of ice (m)
528 :
529 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
530 : afsdn ! floe size distribution tracer
531 :
532 : real (kind=dbl_kind), intent(inout) :: &
533 : vi0new , & ! volume of new ice added to cat 1 (m)
534 : frazil ! frazil ice growth (m/step-->cm/day)
535 :
536 : real (kind=dbl_kind), dimension(ncat), intent(out) :: &
537 : d_an_latg ! change in aicen occuring due to lateral growth
538 :
539 : real (kind=dbl_kind), intent(out) :: &
540 : G_radial , & ! lateral melt rate (m/s)
541 : tot_latg ! total change in aice due to
542 : ! lateral growth at the edges of floes
543 :
544 : ! local variables
545 : integer (kind=int_kind) :: &
546 : n, k ! ice category indices
547 :
548 : real (kind=dbl_kind) :: &
549 : vi0new_lat ! volume of new ice added laterally to FSD (m)
550 :
551 : real (kind=dbl_kind), intent(out) :: &
552 : lead_area , & ! the fractional area of the lead region
553 : latsurf_area ! the area covered by lateral surface of floes
554 :
555 : character(len=*),parameter :: subname='(fsd_lateral_growth)'
556 :
557 38145 : lead_area = c0
558 38145 : latsurf_area = c0
559 38145 : G_radial = c0
560 38145 : tot_latg = c0
561 228870 : d_an_latg = c0
562 :
563 : ! partition volume into lateral growth and frazil
564 : call partition_area (aice, &
565 : aicen, vicen, &
566 : afsdn, lead_area, &
567 38145 : latsurf_area)
568 38145 : if (icepack_warnings_aborted(subname)) return
569 :
570 38145 : vi0new_lat = c0
571 38145 : if (latsurf_area > puny) then
572 38142 : vi0new_lat = vi0new * lead_area / (c1 + aice/latsurf_area)
573 : end if
574 :
575 : ! lateral growth increment
576 38145 : if (vi0new_lat > puny) then
577 37712 : G_radial = vi0new_lat/dt
578 226272 : do n = 1, ncat
579 :
580 2052072 : do k = 1, nfsd
581 : d_an_latg(n) = d_an_latg(n) &
582 2014360 : + c2*aicen(n)*afsdn(k,n)*G_radial*dt/floe_rad_c(k)
583 : end do
584 : end do ! n
585 :
586 : ! cannot expand ice laterally beyond lead region
587 226272 : if (SUM(d_an_latg(:)).ge.lead_area) then
588 0 : d_an_latg(:) = d_an_latg(:)/SUM(d_an_latg(:))
589 0 : d_an_latg(:) = d_an_latg(:)*lead_area
590 : end if
591 :
592 : endif ! vi0new_lat > 0
593 :
594 : ! Use remaining ice volume as in standard model,
595 228870 : tot_latg = SUM(d_an_latg(:))
596 :
597 : end subroutine fsd_lateral_growth
598 :
599 : !=======================================================================
600 : !
601 : ! Evolve the FSD subject to lateral growth and the growth of new ice
602 : ! in thickness category 1.
603 : !
604 : ! If ocean surface wave forcing is provided, the floe size of new ice
605 : ! (grown away from floe edges) can computed from the wave field
606 : ! based on the tensile (stretching) stress limitation following
607 : ! Shen et al. (2001). Otherwise, new floes all grow in the smallest
608 : ! floe size category, representing pancake ice formation.
609 : !
610 : ! Shen, H., Ackley, S., & Hopkins, M. (2001). A conceptual model
611 : ! for pancake-ice formation in a wave field.
612 : ! Annals of Glaciology, 33, 361-367. doi:10.3189/172756401781818239
613 : !
614 : ! authors: Lettie Roach, NIWA/VUW
615 : !
616 188978 : subroutine fsd_add_new_ice (n, &
617 : dt, ai0new, &
618 0 : d_an_latg, d_an_newi, &
619 0 : G_radial, area2, &
620 : wave_sig_ht, &
621 188978 : wave_spectrum, &
622 188978 : wavefreq, &
623 188978 : dwavefreq, &
624 188978 : d_afsd_latg, &
625 0 : d_afsd_newi, &
626 377956 : afsdn, aicen_init, &
627 377956 : aicen, trcrn)
628 :
629 : integer (kind=int_kind), intent(in) :: &
630 : n ! thickness category number
631 :
632 : real (kind=dbl_kind), intent(in) :: &
633 : dt , & ! time step (s)
634 : ai0new , & ! area of new ice added to cat 1
635 : G_radial , & ! lateral melt rate (m/s)
636 : wave_sig_ht ! wave significant height (everywhere) (m)
637 :
638 : real (kind=dbl_kind), dimension(:), intent(in) :: &
639 : wave_spectrum ! ocean surface wave spectrum as a function of frequency
640 : ! power spectral density of surface elevation, E(f) (units m^2 s)
641 :
642 : real(kind=dbl_kind), dimension(:), intent(in) :: &
643 : wavefreq , & ! wave frequencies (s^-1)
644 : dwavefreq ! wave frequency bin widths (s^-1)
645 :
646 : real (kind=dbl_kind), dimension(:), intent(in) :: &
647 : d_an_latg , & ! change in aicen due to lateral growth
648 : d_an_newi ! change in aicen due to frazil ice formation
649 :
650 : real (kind=dbl_kind), dimension (:), intent(in) :: &
651 : aicen_init , & ! fractional area of ice
652 : aicen ! after update
653 :
654 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
655 : afsdn ! floe size distribution tracer
656 :
657 : real (kind=dbl_kind), dimension (:), intent(in) :: &
658 : area2 ! area after lateral growth, before new ice formation
659 :
660 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
661 : trcrn ! ice tracers
662 :
663 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
664 : ! change in floe size distribution (area)
665 : d_afsd_latg, & ! due to fsd lateral growth
666 : d_afsd_newi ! new ice formation
667 :
668 : integer (kind=int_kind) :: &
669 : k, & ! floe size category index
670 : new_size, & ! index for floe size of new ice
671 : nsubt ! number of subtimesteps
672 :
673 : real (kind=dbl_kind) :: &
674 : elapsed_t, subdt ! elapsed time, subtimestep (s)
675 :
676 : real (kind=dbl_kind), dimension (nfsd,ncat) :: &
677 377956 : afsdn_latg ! fsd after lateral growth
678 :
679 : real (kind=dbl_kind), dimension (nfsd) :: &
680 377956 : dafsd_tmp, & ! tmp FSD
681 377956 : df_flx , & ! finite differences for fsd
682 188978 : afsd_ni ! fsd after new ice added
683 :
684 : real (kind=dbl_kind), dimension(nfsd+1) :: &
685 188978 : f_flx ! finite differences in floe size
686 :
687 : character(len=*),parameter :: subname='(fsd_add_new_ice)'
688 :
689 2016472 : afsdn_latg(:,n) = afsdn(:,n) ! default
690 :
691 188978 : if (d_an_latg(n) > puny) then ! lateral growth
692 :
693 : ! adaptive timestep
694 135313 : elapsed_t = c0
695 135313 : nsubt = 0
696 :
697 270626 : DO WHILE (elapsed_t.lt.dt)
698 :
699 135313 : nsubt = nsubt + 1
700 135313 : if (nsubt.gt.100) then
701 0 : write(warnstr,*) subname,'latg not converging'
702 0 : call icepack_warnings_add(warnstr)
703 : endif
704 :
705 : ! finite differences
706 1691551 : df_flx(:) = c0 ! NB could stay zero if all in largest FS cat
707 1826864 : f_flx (:) = c0
708 1556238 : do k = 2, nfsd
709 1550100 : f_flx(k) = G_radial * afsdn_latg(k-1,n) / floe_binwidth(k-1)
710 : end do
711 1691551 : do k = 1, nfsd
712 1691551 : df_flx(k) = f_flx(k+1) - f_flx(k)
713 : end do
714 :
715 1691551 : dafsd_tmp(:) = c0
716 1691551 : do k = 1, nfsd
717 : dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) &
718 20298889 : * (c1/floe_rad_c(k) - SUM(afsdn_latg(:,n)/floe_rad_c(:))) )
719 :
720 : end do
721 :
722 : ! timestep required for this
723 135313 : subdt = get_subdt_fsd(afsdn_latg(:,n), dafsd_tmp(:))
724 135313 : subdt = MIN(subdt, dt)
725 :
726 : ! update fsd and elapsed time
727 1691551 : afsdn_latg(:,n) = afsdn_latg(:,n) + subdt*dafsd_tmp(:)
728 135313 : elapsed_t = elapsed_t + subdt
729 :
730 : END DO
731 :
732 135313 : call icepack_cleanup_fsdn (afsdn_latg(:,n))
733 135313 : if (icepack_warnings_aborted(subname)) return
734 1691551 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n)
735 :
736 : end if ! lat growth
737 :
738 188978 : new_size = nfsd
739 188978 : if (n == 1) then
740 : ! add new frazil ice to smallest thickness
741 38145 : if (d_an_newi(n) > puny) then
742 :
743 405124 : afsd_ni(:) = c0
744 :
745 405124 : if (SUM(afsdn_latg(:,n)) > puny) then ! fsd exists
746 :
747 38139 : if (wave_spec) then
748 29890 : if (wave_sig_ht > puny) then
749 : call wave_dep_growth (wave_spectrum, &
750 : wavefreq, dwavefreq, &
751 29890 : new_size)
752 29890 : if (icepack_warnings_aborted(subname)) return
753 : end if
754 :
755 : ! grow in new_size category
756 : afsd_ni(new_size) = (afsdn_latg(new_size,n)*area2(n) + ai0new) &
757 29890 : / (area2(n) + ai0new)
758 29890 : do k = 1, new_size-1 ! diminish other floe cats accordingly
759 0 : afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new)
760 : end do
761 358680 : do k = new_size+1, nfsd ! diminish other floe cats accordingly
762 358680 : afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new)
763 : end do
764 :
765 : else ! grow in smallest floe size category
766 : afsd_ni(1) = (afsdn_latg(1,n)*area2(n) + ai0new) &
767 8249 : / (area2(n) + ai0new)
768 8249 : do k = 2, nfsd ! diminish other floe cats accordingly
769 0 : afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n)+ai0new)
770 : enddo
771 : end if ! wave spec
772 :
773 : else ! no fsd, so entirely new ice
774 :
775 6 : if (wave_spec) then
776 4 : if (wave_sig_ht > puny) then
777 : call wave_dep_growth (wave_spectrum, &
778 : wavefreq, dwavefreq, &
779 4 : new_size)
780 4 : if (icepack_warnings_aborted(subname)) return
781 : end if
782 :
783 4 : afsd_ni(new_size) = c1
784 : else
785 2 : afsd_ni(1) = c1
786 : endif ! wave forcing
787 :
788 : endif ! entirely new ice or not
789 :
790 405124 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_ni(:)
791 38145 : call icepack_cleanup_fsdn (trcrn(nt_fsd:nt_fsd+nfsd-1,n))
792 38145 : if (icepack_warnings_aborted(subname)) return
793 : endif ! d_an_newi > puny
794 : endif ! n = 1
795 :
796 : ! history/diagnostics
797 2016472 : do k = 1, nfsd
798 : ! sum over n
799 : d_afsd_latg(k) = d_afsd_latg(k) &
800 : + area2(n)*afsdn_latg(k,n) & ! after latg
801 1827494 : - aicen_init(n)*afsdn(k,n) ! at start
802 : d_afsd_newi(k) = d_afsd_newi(k) &
803 : + aicen(n)*trcrn(nt_fsd+k-1,n) & ! after latg and newi
804 2016472 : - area2(n)*afsdn_latg(k,n) ! after latg
805 : enddo ! k
806 :
807 : end subroutine fsd_add_new_ice
808 :
809 : !=======================================================================
810 : !
811 : ! Given a wave spectrum, calculate size of new floes based on
812 : ! tensile failure, following Shen et al. (2001)
813 : !
814 : ! The tensile mode parameter is based on in-situ measurements
815 : ! by Roach, Smith & Dean (2018).
816 : !
817 : ! authors: Lettie Roach, NIWA/VUW
818 : !
819 29894 : subroutine wave_dep_growth (local_wave_spec, &
820 29894 : wavefreq, dwavefreq, &
821 : new_size)
822 :
823 : real (kind=dbl_kind), dimension(:), intent(in) :: &
824 : local_wave_spec ! ocean surface wave spectrum as a function of frequency
825 : ! power spectral density of surface elevation, E(f) (units m^2 s)
826 : ! dimension set in ice_forcing
827 :
828 : real(kind=dbl_kind), dimension(:), intent(in) :: &
829 : wavefreq, & ! wave frequencies (s^-1)
830 : dwavefreq ! wave frequency bin widths (s^-1)
831 :
832 : integer (kind=int_kind), intent(out) :: &
833 : new_size ! index of floe size category in which new floes will grow
834 :
835 : ! local variables
836 : real (kind=dbl_kind), parameter :: &
837 : tensile_param = 0.167_dbl_kind ! tensile mode parameter (kg m^-1 s^-2)
838 : ! value from Roach, Smith & Dean (2018)
839 :
840 : real (kind=dbl_kind) :: &
841 : w_amp, & ! wave amplitude (m)
842 : f_peak, & ! peak frequency (s^-1)
843 : r_max ! floe radius (m)
844 :
845 : integer (kind=int_kind) :: k
846 :
847 777244 : w_amp = c2* SQRT(SUM(local_wave_spec*dwavefreq)) ! sig wave amplitude
848 777244 : f_peak = wavefreq(MAXLOC(local_wave_spec, DIM=1)) ! peak frequency
849 :
850 : ! tensile failure
851 29894 : if (w_amp > puny .and. f_peak > puny) then
852 29894 : r_max = p5*SQRT(tensile_param*gravit/(pi**5*rhoi*w_amp*2))/f_peak**2
853 : else
854 0 : r_max = bignum
855 : end if
856 :
857 29894 : new_size = nfsd
858 358728 : do k = nfsd-1, 1, -1
859 358728 : if (r_max <= floe_rad_h(k)) new_size = k
860 : end do
861 :
862 29894 : end subroutine wave_dep_growth
863 :
864 : !=======================================================================
865 : !
866 : ! Floes are perimitted to weld together in freezing conditions, according
867 : ! to their geometric probability of overlap if placed randomly on the
868 : ! domain. The rate per unit area c_weld is the total number
869 : ! of floes that weld with another, per square meter, per unit time, in the
870 : ! case of a fully covered ice surface (aice=1), equal to twice the reduction
871 : ! in total floe number. See Roach, Smith & Dean (2018).
872 : !
873 : !
874 : ! authors: Lettie Roach, NIWA/VUW
875 : !
876 :
877 98676 : subroutine fsd_weld_thermo (dt, frzmlt, &
878 98676 : aicen, trcrn, &
879 98676 : d_afsd_weld)
880 :
881 : real (kind=dbl_kind), intent(in) :: &
882 : dt ! time step (s)
883 :
884 : real (kind=dbl_kind), dimension (:), intent(in) :: &
885 : aicen ! ice concentration
886 :
887 : real (kind=dbl_kind), intent(in) :: &
888 : frzmlt ! freezing/melting potential (W/m^2)
889 :
890 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
891 : trcrn ! ice tracers
892 :
893 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
894 : d_afsd_weld ! change in fsd due to welding
895 :
896 : ! local variables
897 : real (kind=dbl_kind), parameter :: &
898 : aminweld = p1 ! minimum ice concentration likely to weld
899 :
900 : real (kind=dbl_kind), parameter :: &
901 : c_weld = 1.0e-8_dbl_kind
902 : ! constant of proportionality for welding
903 : ! total number of floes that weld with another, per square meter,
904 : ! per unit time, in the case of a fully covered ice surface
905 : ! units m^-2 s^-1, see documentation for details
906 :
907 : integer (kind=int_kind) :: &
908 : n , & ! thickness category index
909 : k, i, j ! floe size category indices
910 :
911 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
912 197352 : afsdn ! floe size distribution tracer
913 :
914 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
915 197352 : d_afsdn_weld ! change in afsdn due to welding
916 :
917 : real (kind=dbl_kind), dimension(nfsd) :: &
918 197352 : stability , & ! check for stability
919 197352 : nfsd_tmp , & ! number fsd
920 98676 : afsd_init , & ! initial values
921 197352 : afsd_tmp , & ! work array
922 98676 : gain, loss ! welding tendencies
923 :
924 : real(kind=dbl_kind) :: &
925 : kern , & ! kernel
926 : subdt , & ! subcycling time step for stability (s)
927 : elapsed_t ! elapsed subcycling time
928 :
929 : character(len=*), parameter :: subname='(fsd_weld_thermo)'
930 :
931 :
932 5067216 : afsdn (:,:) = c0
933 993708 : afsd_init(:) = c0
934 993708 : stability = c0
935 :
936 592056 : do n = 1, ncat
937 :
938 4968540 : d_afsd_weld (:) = c0
939 4968540 : d_afsdn_weld(:,n) = c0
940 4968540 : afsdn(:,n) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
941 493380 : call icepack_cleanup_fsdn (afsdn(:,n))
942 493380 : if (icepack_warnings_aborted(subname)) return
943 :
944 : ! If there is some ice in the lower (nfsd-1) categories
945 : ! and there is freezing potential
946 : if ((frzmlt > puny) .and. & ! freezing potential
947 4475160 : (aicen(n) > aminweld) .and. & ! low concentrations unlikely to weld
948 98676 : (SUM(afsdn(1:nfsd-1,n)) > puny)) then ! some ice in nfsd-1 categories
949 :
950 1166139 : afsd_init(:) = afsdn(:,n) ! save initial values
951 1166139 : afsd_tmp (:) = afsd_init(:) ! work array
952 :
953 : ! in case of minor numerical errors
954 1166139 : WHERE(afsd_tmp < puny) afsd_tmp = c0
955 2242575 : afsd_tmp = afsd_tmp/SUM(afsd_tmp)
956 :
957 : ! adaptive sub-timestep
958 89703 : elapsed_t = c0
959 1977052 : DO WHILE (elapsed_t < dt)
960 :
961 : ! calculate sub timestep
962 24535537 : nfsd_tmp = afsd_tmp/floe_area_c
963 : WHERE (afsd_tmp > puny) &
964 24535537 : stability = nfsd_tmp/(c_weld*afsd_tmp*aicen(n))
965 24535537 : WHERE (stability < puny) stability = bignum
966 24535537 : subdt = MINVAL(stability)
967 1887349 : subdt = MIN(subdt,dt)
968 :
969 24535537 : loss(:) = c0
970 24535537 : gain(:) = c0
971 :
972 24535537 : do i = 1, nfsd ! consider loss from this category
973 296313793 : do j = 1, nfsd ! consider all interaction partners
974 271778256 : k = iweld(i,j) ! product of i+j
975 294426444 : if (k > i) then
976 156649967 : kern = c_weld * floe_area_c(i) * aicen(n)
977 156649967 : loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j)
978 156649967 : gain(k) = gain(k) + kern*afsd_tmp(i)*afsd_tmp(j)
979 : end if
980 : end do
981 : end do
982 :
983 : ! does largest category lose?
984 : ! if (loss(nfsd) > puny) stop 'weld, largest cat losing'
985 : ! if (gain(1) > puny) stop 'weld, smallest cat gaining'
986 :
987 : ! update afsd
988 24535537 : afsd_tmp(:) = afsd_tmp(:) + subdt*(gain(:) - loss(:))
989 :
990 : ! in case of minor numerical errors
991 24535537 : WHERE(afsd_tmp < puny) afsd_tmp = c0
992 47183725 : afsd_tmp = afsd_tmp/SUM(afsd_tmp)
993 :
994 : ! update time
995 1887349 : elapsed_t = elapsed_t + subdt
996 :
997 : ! stop if all in largest floe size cat
998 1887349 : if (afsd_tmp(nfsd) > (c1-puny)) exit
999 :
1000 : END DO ! time
1001 :
1002 1166139 : afsdn(:,n) = afsd_tmp(:)
1003 89703 : call icepack_cleanup_fsdn (afsdn(:,n))
1004 89703 : if (icepack_warnings_aborted(subname)) return
1005 :
1006 1166139 : do k = 1, nfsd
1007 1076436 : trcrn(nt_fsd+k-1,n) = afsdn(k,n)
1008 : ! history/diagnostics
1009 1166139 : d_afsdn_weld(k,n) = afsdn(k,n) - afsd_init(k)
1010 : enddo
1011 : endif ! try to weld
1012 : enddo ! ncat
1013 :
1014 : ! history/diagnostics
1015 993708 : do k = 1, nfsd
1016 895032 : d_afsd_weld(k) = c0
1017 5468868 : do n = 1, ncat
1018 5370192 : d_afsd_weld(k) = d_afsd_weld(k) + aicen(n)*d_afsdn_weld(k,n)
1019 : end do ! n
1020 : end do ! k
1021 :
1022 : end subroutine fsd_weld_thermo
1023 :
1024 : !=======================================================================
1025 : !
1026 : ! Adaptive timestepping (process agnostic)
1027 : ! See reference: Horvat & Tziperman (2017) JGR, Appendix A
1028 : !
1029 : ! authors: 2018 Lettie Roach, NIWA/VUW
1030 : !
1031 : !
1032 417905 : function get_subdt_fsd(afsd_init, d_afsd) &
1033 : result(subdt)
1034 :
1035 : real (kind=dbl_kind), dimension (nfsd), intent(in) :: &
1036 : afsd_init, d_afsd ! floe size distribution tracer
1037 :
1038 : ! output
1039 : real (kind=dbl_kind) :: &
1040 : subdt ! subcycle timestep (s)
1041 :
1042 : ! local variables
1043 : real (kind=dbl_kind), dimension (nfsd) :: &
1044 835810 : check_dt ! to compute subcycle timestep (s)
1045 :
1046 : integer (kind=int_kind) :: k
1047 :
1048 4453930 : check_dt(:) = bignum
1049 4453930 : do k = 1, nfsd
1050 4036025 : if (d_afsd(k) > puny) check_dt(k) = (1-afsd_init(k))/d_afsd(k)
1051 4453930 : if (d_afsd(k) < -puny) check_dt(k) = afsd_init(k)/ABS(d_afsd(k))
1052 : end do
1053 :
1054 4453930 : subdt = MINVAL(check_dt)
1055 :
1056 417905 : end function get_subdt_fsd
1057 :
1058 :
1059 : !=======================================================================
1060 :
1061 : end module icepack_fsd
1062 :
1063 : !=======================================================================
1064 :
|
