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