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