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 198 : 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 12 : 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 520 : floe_rad
122 :
123 : real (kind=dbl_kind), dimension(:), allocatable :: &
124 198 : 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 198 : l_write_diags = .true.
134 198 : if (present(write_diags)) then
135 198 : l_write_diags = write_diags
136 : endif
137 :
138 198 : 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, & ! LCOV_EXCL_LINE
144 : 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, & ! LCOV_EXCL_LINE
145 : 9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, & ! LCOV_EXCL_LINE
146 : 3.35434988e+03, 4.55051413e+03, 6.17323164e+03, 8.37461170e+03, & ! LCOV_EXCL_LINE
147 : 1.13610059e+04, 1.54123510e+04, 2.09084095e+04, 2.83643675e+04, & ! LCOV_EXCL_LINE
148 0 : 3.84791270e+04 /)
149 :
150 198 : 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, & ! LCOV_EXCL_LINE
156 : 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, & ! LCOV_EXCL_LINE
157 : 9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, & ! LCOV_EXCL_LINE
158 0 : 3.35434988e+03 /)
159 :
160 198 : else if (nfsd.eq.12) then
161 :
162 194 : 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, & ! LCOV_EXCL_LINE
166 : 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, & ! LCOV_EXCL_LINE
167 2716 : 9.45812834e+02 /)
168 :
169 4 : else if (nfsd.eq.1) then ! default case
170 :
171 4 : allocate(lims(1+1))
172 :
173 12 : 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) ! LCOV_EXCL_LINE
186 : floe_area_l (nfsd), & ! fsd area at lower bound (m^2) ! LCOV_EXCL_LINE
187 : floe_area_h (nfsd), & ! fsd area at higher bound (m^2) ! LCOV_EXCL_LINE
188 : floe_area_c (nfsd), & ! fsd area at bin centre (m^2) ! LCOV_EXCL_LINE
189 : floe_area_binwidth (nfsd), & ! floe area bin width (m^2) ! LCOV_EXCL_LINE
190 : iweld (nfsd, nfsd), & ! fsd categories that can weld ! LCOV_EXCL_LINE
191 198 : stat=ierr)
192 198 : 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 2530 : floe_rad_l = lims(1:nfsd )
199 2530 : floe_rad_h = lims(2:nfsd+1)
200 2530 : floe_rad_c = (floe_rad_h+floe_rad_l)/c2
201 :
202 2530 : floe_area_l = c4*floeshape*floe_rad_l**2
203 2530 : floe_area_c = c4*floeshape*floe_rad_c**2
204 2530 : floe_area_h = c4*floeshape*floe_rad_h**2
205 :
206 2530 : floe_binwidth = floe_rad_h - floe_rad_l
207 :
208 2530 : floe_area_binwidth = floe_area_h - floe_area_l
209 :
210 : ! floe size categories that can combine during welding
211 30470 : iweld(:,:) = -999
212 2530 : do n = 1, nfsd
213 30470 : do m = 1, nfsd
214 27940 : test = floe_area_c(n) + floe_area_c(m)
215 335236 : do k = 1, nfsd-1
216 307296 : if ((test >= floe_area_l(k)) .and. (test < floe_area_h(k))) &
217 50444 : iweld(n,m) = k
218 : end do
219 30272 : if (test >= floe_area_l(nfsd)) iweld(n,m) = nfsd
220 : end do
221 : end do
222 :
223 198 : if (allocated(lims)) deallocate(lims)
224 :
225 : ! write fsd bounds
226 198 : floe_rad(0) = floe_rad_l(1)
227 2530 : do n = 1, nfsd
228 2332 : floe_rad(n) = floe_rad_h(n)
229 : ! Save character string to write to history file
230 2332 : write (c_nf, '(i2)') n
231 2332 : write (c_fsd1, '(f6.3)') floe_rad(n-1)
232 2332 : write (c_fsd2, '(f6.3)') floe_rad(n)
233 2530 : c_fsd_range(n)=c_fsd1//'m < fsd Cat '//c_nf//' < '//c_fsd2//'m'
234 : enddo
235 :
236 198 : if (l_write_diags) then
237 18 : write(warnstr,*) ' '
238 18 : call icepack_warnings_add(warnstr)
239 18 : write(warnstr,*) subname
240 18 : call icepack_warnings_add(warnstr)
241 18 : write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
242 18 : call icepack_warnings_add(warnstr)
243 223 : do n = 1, nfsd
244 205 : write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
245 223 : call icepack_warnings_add(warnstr)
246 : enddo
247 18 : write(warnstr,*) ' '
248 18 : call icepack_warnings_add(warnstr)
249 : endif
250 :
251 198 : 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 108 : subroutine icepack_init_fsd(nfsd, ice_ic, &
272 : floe_rad_c, & ! fsd size bin centre in m (radius) ! LCOV_EXCL_LINE
273 : floe_binwidth, & ! fsd size bin width in m (radius) ! LCOV_EXCL_LINE
274 116 : 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) ! LCOV_EXCL_LINE
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 16 : real (kind=dbl_kind) :: alpha, totfrac
294 :
295 : integer (kind=int_kind) :: k
296 :
297 : real (kind=dbl_kind), dimension (nfsd) :: &
298 276 : num_fsd ! number distribution of floes
299 :
300 116 : 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 116 : alpha = 2.1_dbl_kind
309 116 : totfrac = c0 ! total fraction of floes
310 1464 : do k = 1, nfsd
311 1348 : num_fsd(k) = (2*floe_rad_c(k))**(-alpha-c1) ! number distribution of floes
312 1348 : afsd (k) = num_fsd(k)*floe_area_c(k)*floe_binwidth(k) ! fraction distribution of floes
313 1464 : totfrac = totfrac + afsd(k)
314 : enddo
315 1464 : afsd = afsd/totfrac ! normalize
316 :
317 : endif ! ice_ic
318 :
319 116 : 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 67686183 : subroutine icepack_cleanup_fsd (ncat, nfsd, afsdn)
329 :
330 : integer (kind=int_kind), intent(in) :: &
331 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
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 67686183 : if (tr_fsd) then
347 :
348 406117098 : do n = 1, ncat
349 338430915 : call icepack_cleanup_fsdn(nfsd, afsdn(:,n))
350 406117098 : 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 464807120 : 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 78624321 : tot
379 :
380 5829273455 : do k = 1, nfsd
381 5829273455 : if (afsd(k) < puny) afsd(k) = c0
382 : enddo
383 :
384 5829273455 : tot = sum(afsd(:))
385 464807120 : if (tot > puny) then
386 1432614744 : do k = 1, nfsd
387 1432614744 : afsd(k) = afsd(k) / tot ! normalize
388 : enddo
389 : else ! represents ice-free cell, so set to zero
390 4396658711 : afsd(:) = c0
391 : endif
392 :
393 464807120 : 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 2213453 : subroutine partition_area (ncat, nfsd, &
404 : floe_rad_c, aice, & ! LCOV_EXCL_LINE
405 : aicen, vicen, & ! LCOV_EXCL_LINE
406 : afsdn, lead_area, & ! LCOV_EXCL_LINE
407 : latsurf_area)
408 :
409 : integer (kind=int_kind), intent(in) :: &
410 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
434 : k ! floe size index
435 :
436 : real (kind=dbl_kind) :: &
437 : width_leadreg, & ! width of lead region (m) ! LCOV_EXCL_LINE
438 337297 : thickness ! actual thickness of ice in thickness cat (m)
439 :
440 : !-----------------------------------------------------------------
441 : ! Initialize
442 : !-----------------------------------------------------------------
443 2213453 : lead_area = c0
444 2213453 : 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 2213453 : width_leadreg = floe_rad_c(1)
449 :
450 : ! Only calculate these areas where there is some ice
451 2213453 : if (aice > puny) then
452 :
453 : ! lead area = sum of areas of annuli around all floes
454 13279386 : do n = 1, ncat
455 141468866 : do k = 1, nfsd
456 15633320 : lead_area = lead_area + aicen(n) * afsdn(k,n) &
457 : * (c2*width_leadreg /floe_rad_c(k) & ! LCOV_EXCL_LINE
458 139255635 : + 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 2213231 : lead_area=MIN(lead_area, c1-aice)
464 2213231 : if (lead_area < c0) then
465 3656 : if (lead_area < -puny) then
466 : ! stop 'lead_area lt0 in partition_area'
467 : else
468 3656 : 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 13279386 : do n = 1, ncat
474 11066155 : thickness = c0
475 11066155 : if (aicen(n) > c0) thickness = vicen(n)/aicen(n)
476 :
477 141468866 : do k = 1, nfsd
478 : latsurf_area = latsurf_area &
479 139255635 : + 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 2213453 : 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 2213453 : subroutine fsd_lateral_growth (ncat, nfsd, &
504 : dt, aice, & ! LCOV_EXCL_LINE
505 : aicen, vicen, & ! LCOV_EXCL_LINE
506 : vi0new, & ! LCOV_EXCL_LINE
507 : frazil, floe_rad_c, & ! LCOV_EXCL_LINE
508 : afsdn, & ! LCOV_EXCL_LINE
509 : lead_area, latsurf_area, & ! LCOV_EXCL_LINE
510 : G_radial, d_an_latg, & ! LCOV_EXCL_LINE
511 : tot_latg)
512 :
513 : integer (kind=int_kind), intent(in) :: &
514 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
515 : nfsd ! number of floe size categories
516 :
517 : real (kind=dbl_kind), intent(in) :: &
518 : dt , & ! time step (s) ! LCOV_EXCL_LINE
519 : aice ! total concentration of ice
520 :
521 : real (kind=dbl_kind), dimension (:), intent(in) :: &
522 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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 337297 : 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 ! LCOV_EXCL_LINE
553 : latsurf_area ! the area covered by lateral surface of floes
554 :
555 : character(len=*),parameter :: subname='(fsd_lateral_growth)'
556 :
557 2213453 : lead_area = c0
558 2213453 : latsurf_area = c0
559 2213453 : G_radial = c0
560 2213453 : tot_latg = c0
561 13280718 : d_an_latg = c0
562 :
563 : ! partition volume into lateral growth and frazil
564 : call partition_area (ncat, nfsd, &
565 : floe_rad_c, aice, & ! LCOV_EXCL_LINE
566 : aicen, vicen, & ! LCOV_EXCL_LINE
567 : afsdn, lead_area, & ! LCOV_EXCL_LINE
568 2213453 : latsurf_area)
569 2213453 : if (icepack_warnings_aborted(subname)) return
570 :
571 2213453 : vi0new_lat = c0
572 2213453 : if (latsurf_area > puny) then
573 2213231 : vi0new_lat = vi0new * lead_area / (c1 + aice/latsurf_area)
574 : end if
575 :
576 : ! for history/diagnostics
577 2213453 : frazil = vi0new - vi0new_lat
578 :
579 : ! lateral growth increment
580 2213453 : if (vi0new_lat > puny) then
581 2053906 : G_radial = vi0new_lat/dt
582 12323436 : do n = 1, ncat
583 :
584 131730511 : do k = 1, nfsd
585 14493235 : d_an_latg(n) = d_an_latg(n) &
586 129676605 : + 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 12323436 : 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 2213453 : vi0new = vi0new - vi0new_lat
601 13280718 : 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 10321314 : subroutine fsd_add_new_ice (ncat, n, nfsd, &
623 : dt, ai0new, & ! LCOV_EXCL_LINE
624 : d_an_latg, d_an_newi, & ! LCOV_EXCL_LINE
625 : floe_rad_c, floe_binwidth, & ! LCOV_EXCL_LINE
626 : G_radial, area2, & ! LCOV_EXCL_LINE
627 : wave_sig_ht, & ! LCOV_EXCL_LINE
628 : wave_spectrum, & ! LCOV_EXCL_LINE
629 : wavefreq, & ! LCOV_EXCL_LINE
630 : dwavefreq, & ! LCOV_EXCL_LINE
631 : d_afsd_latg, & ! LCOV_EXCL_LINE
632 : d_afsd_newi, & ! LCOV_EXCL_LINE
633 : afsdn, aicen_init, & ! LCOV_EXCL_LINE
634 20642628 : aicen, trcrn)
635 :
636 : integer (kind=int_kind), intent(in) :: &
637 : n , & ! thickness category number ! LCOV_EXCL_LINE
638 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
639 : nfsd ! number of floe size categories
640 :
641 : real (kind=dbl_kind), intent(in) :: &
642 : dt , & ! time step (s) ! LCOV_EXCL_LINE
643 : ai0new , & ! area of new ice added to cat 1 ! LCOV_EXCL_LINE
644 : G_radial , & ! lateral melt rate (m/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
661 : aicen , & ! after update ! LCOV_EXCL_LINE
662 : floe_rad_c , & ! fsd size bin centre in m (radius) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
681 : new_size, & ! index for floe size of new ice ! LCOV_EXCL_LINE
682 : nsubt ! number of subtimesteps
683 :
684 : real (kind=dbl_kind) :: &
685 3081140 : elapsed_t, subdt ! elapsed time, subtimestep (s)
686 :
687 : real (kind=dbl_kind), dimension (nfsd,ncat) :: &
688 99595473 : afsdn_latg ! fsd after lateral growth
689 :
690 : real (kind=dbl_kind), dimension (nfsd) :: &
691 : dafsd_tmp, & ! tmp FSD ! LCOV_EXCL_LINE
692 : df_flx , & ! finite differences for fsd ! LCOV_EXCL_LINE
693 36741311 : afsd_ni ! fsd after new ice added
694 :
695 : real (kind=dbl_kind), dimension(nfsd+1) :: &
696 27960567 : f_flx ! finite differences in floe size
697 :
698 : character(len=*),parameter :: subname='(fsd_add_new_ice)'
699 :
700 130248355 : afsdn_latg(:,n) = afsdn(:,n) ! default
701 :
702 10321314 : if (d_an_latg(n) > puny) then ! lateral growth
703 :
704 : ! adaptive timestep
705 7248820 : elapsed_t = c0
706 7248820 : nsubt = 0
707 :
708 14497640 : DO WHILE (elapsed_t.lt.dt)
709 :
710 7248820 : nsubt = nsubt + 1
711 7248820 : if (nsubt.gt.100) print *, 'latg not converging'
712 :
713 : ! finite differences
714 93977469 : df_flx(:) = c0 ! NB could stay zero if all in largest FS cat
715 101226289 : f_flx (:) = c0
716 86728649 : do k = 2, nfsd
717 86728649 : f_flx(k) = G_radial * afsdn_latg(k-1,n) / floe_binwidth(k-1)
718 : end do
719 93977469 : do k = 1, nfsd
720 93977469 : 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 93977469 : dafsd_tmp(:) = c0
726 93977469 : do k = 1, nfsd
727 19722490 : dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) &
728 1144325311 : * (c1/floe_rad_c(k) - SUM(afsdn_latg(:,n)/floe_rad_c(:))) )
729 :
730 : end do
731 :
732 : ! timestep required for this
733 7248820 : subdt = get_subdt_fsd(nfsd, afsdn_latg(:,n), dafsd_tmp(:))
734 7248820 : subdt = MIN(subdt, dt)
735 :
736 : ! update fsd and elapsed time
737 93977469 : afsdn_latg(:,n) = afsdn_latg(:,n) + subdt*dafsd_tmp(:)
738 14497640 : elapsed_t = elapsed_t + subdt
739 :
740 : END DO
741 :
742 7248820 : call icepack_cleanup_fsdn (nfsd, afsdn_latg(:,n))
743 7248820 : if (icepack_warnings_aborted(subname)) return
744 93977469 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n)
745 :
746 : end if ! lat growth
747 :
748 10321314 : new_size = nfsd
749 10321314 : if (n == 1) then
750 : ! add new frazil ice to smallest thickness
751 2112131 : if (d_an_newi(n) > puny) then
752 :
753 26590749 : afsd_ni(:) = c0
754 :
755 26590749 : if (SUM(afsdn_latg(:,n)) > puny) then ! fsd exists
756 :
757 2111862 : if (wave_spec) then
758 2033050 : if (wave_sig_ht > puny) then
759 0 : call wave_dep_growth (nfsd, wave_spectrum, &
760 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
761 14935 : new_size)
762 14935 : if (icepack_warnings_aborted(subname)) return
763 : end if
764 :
765 : ! grow in new_size category
766 721173 : afsd_ni(new_size) = (afsdn_latg(new_size,n)*area2(n) + ai0new) &
767 2513832 : / (area2(n) + ai0new)
768 24305677 : do k = 1, new_size-1 ! diminish other floe cats accordingly
769 24305677 : afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new)
770 : end do
771 2123973 : do k = new_size+1, nfsd ! diminish other floe cats accordingly
772 2123973 : 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 236436 : afsd_ni(1) = (afsdn_latg(1,n)*area2(n) + ai0new) &
777 236436 : / (area2(n) + ai0new)
778 78812 : do k = 2, nfsd ! diminish other floe cats accordingly
779 78812 : 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 269 : if (wave_spec) then
786 267 : if (wave_sig_ht > puny) then
787 0 : call wave_dep_growth (nfsd, wave_spectrum, &
788 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
789 0 : new_size)
790 0 : if (icepack_warnings_aborted(subname)) return
791 : end if
792 :
793 267 : 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 26590749 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_ni(:)
801 2112131 : call icepack_cleanup_fsdn (nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,n))
802 2112131 : if (icepack_warnings_aborted(subname)) return
803 : endif ! d_an_newi > puny
804 : endif ! n = 1
805 :
806 : ! history/diagnostics
807 130248355 : do k = 1, nfsd
808 : ! sum over n
809 14558113 : d_afsd_latg(k) = d_afsd_latg(k) &
810 : + area2(n)*afsdn_latg(k,n) & ! after latg ! LCOV_EXCL_LINE
811 149043267 : - aicen_init(n)*afsdn(k,n) ! at start
812 14558113 : d_afsd_newi(k) = d_afsd_newi(k) &
813 : + aicen(n)*trcrn(nt_fsd+k-1,n) & ! after latg and newi ! LCOV_EXCL_LINE
814 173922694 : - 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 14935 : subroutine wave_dep_growth (nfsd, local_wave_spec, &
830 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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 : w_amp, & ! wave amplitude (m) ! LCOV_EXCL_LINE
855 : f_peak, & ! peak frequency (s^-1) ! LCOV_EXCL_LINE
856 0 : r_max ! floe radius (m)
857 :
858 : integer (kind=int_kind) :: k
859 :
860 388310 : w_amp = c2* SQRT(SUM(local_wave_spec*dwavefreq)) ! sig wave amplitude
861 388310 : f_peak = wavefreq(MAXLOC(local_wave_spec, DIM=1)) ! peak frequency
862 :
863 : ! tensile failure
864 14935 : if (w_amp > puny .and. f_peak > puny) then
865 14935 : 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 14935 : new_size = nfsd
871 179220 : do k = nfsd-1, 1, -1
872 179220 : if (r_max <= floe_rad_h(k)) new_size = k
873 : end do
874 :
875 14935 : 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 22480848 : subroutine fsd_weld_thermo (ncat, nfsd, &
891 : dt, frzmlt, & ! LCOV_EXCL_LINE
892 : aicen, trcrn, & ! LCOV_EXCL_LINE
893 22480848 : d_afsd_weld)
894 :
895 : integer (kind=int_kind), intent(in) :: &
896 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
927 : n , & ! thickness category index ! LCOV_EXCL_LINE
928 : k, kx, ky, i, j ! floe size category indices
929 :
930 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
931 238066416 : afsdn ! floe size distribution tracer
932 :
933 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
934 238066416 : d_afsdn_weld ! change in afsdn due to welding
935 :
936 : real (kind=dbl_kind), dimension(nfsd) :: &
937 : stability , & ! check for stability ! LCOV_EXCL_LINE
938 : nfsd_tmp , & ! number fsd ! LCOV_EXCL_LINE
939 : afsd_init , & ! initial values ! LCOV_EXCL_LINE
940 : afsd_tmp , & ! work array ! LCOV_EXCL_LINE
941 101259888 : gain, loss ! welding tendencies
942 :
943 : real(kind=dbl_kind) :: &
944 : prefac , & ! multiplies kernel ! LCOV_EXCL_LINE
945 : kern , & ! kernel ! LCOV_EXCL_LINE
946 : subdt , & ! subcycling time step for stability (s) ! LCOV_EXCL_LINE
947 3842880 : elapsed_t ! elapsed subcycling time
948 :
949 : character(len=*), parameter :: subname='(fsd_weld_thermo)'
950 :
951 :
952 1430896368 : afsdn (:,:) = c0
953 281683104 : afsd_init(:) = c0
954 281683104 : stability = c0
955 22480848 : prefac = p5
956 :
957 134885088 : do n = 1, ncat
958 :
959 1408415520 : d_afsd_weld (:) = c0
960 1408415520 : d_afsdn_weld(:,n) = c0
961 1408415520 : afsdn(:,n) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
962 112404240 : call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
963 112404240 : 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 : (aicen(n) > aminweld) .and. & ! low concentrations unlikely to weld ! LCOV_EXCL_LINE
969 41695248 : (SUM(afsdn(1:nfsd-1,n)) > puny)) then ! some ice in nfsd-1 categories
970 :
971 59943182 : afsd_init(:) = afsdn(:,n) ! save initial values
972 59943182 : afsd_tmp (:) = afsd_init(:) ! work array
973 :
974 : ! in case of minor numerical errors
975 59943182 : WHERE(afsd_tmp < puny) afsd_tmp = c0
976 115275350 : afsd_tmp = afsd_tmp/SUM(afsd_tmp)
977 :
978 : ! adaptive sub-timestep
979 4611014 : elapsed_t = c0
980 137302863 : DO WHILE (elapsed_t < dt)
981 :
982 : ! calculate sub timestep
983 1724994037 : nfsd_tmp = afsd_tmp/floe_area_c
984 0 : WHERE (afsd_tmp > puny) &
985 1724994037 : stability = nfsd_tmp/(c_weld*afsd_tmp*aicen(n))
986 1724994037 : WHERE (stability < puny) stability = bignum
987 1724994037 : subdt = MINVAL(stability)
988 132691849 : subdt = MIN(subdt,dt)
989 :
990 1724994037 : loss(:) = c0
991 1724994037 : gain(:) = c0
992 :
993 1724994037 : do i = 1, nfsd ! consider loss from this category
994 20832620293 : do j = 1, nfsd ! consider all interaction partners
995 19107626256 : k = iweld(i,j) ! product of i+j
996 20699928444 : if (k > i) then
997 10615347920 : kern = c_weld * floe_area_c(i) * aicen(n)
998 10615347920 : loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j)
999 10615347920 : if (i.eq.j) prefac = c1 ! otherwise 0.5
1000 10615347920 : 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 1724994037 : afsd_tmp(:) = afsd_tmp(:) + subdt*(gain(:) - loss(:))
1011 :
1012 : ! in case of minor numerical errors
1013 1724994037 : WHERE(afsd_tmp < puny) afsd_tmp = c0
1014 3317296225 : afsd_tmp = afsd_tmp/SUM(afsd_tmp)
1015 :
1016 : ! update time
1017 132691849 : elapsed_t = elapsed_t + subdt
1018 :
1019 : ! stop if all in largest floe size cat
1020 137302863 : if (afsd_tmp(nfsd) > (c1-puny)) exit
1021 :
1022 : END DO ! time
1023 :
1024 4611014 : call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
1025 4611014 : if (icepack_warnings_aborted(subname)) return
1026 :
1027 59943182 : do k = 1, nfsd
1028 55332168 : afsdn(k,n) = afsd_tmp(k)
1029 55332168 : trcrn(nt_fsd+k-1,n) = afsdn(k,n)
1030 : ! history/diagnostics
1031 59943182 : 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 281683104 : do k = 1, nfsd
1038 259202256 : d_afsd_weld(k) = c0
1039 1577694384 : do n = 1, ncat
1040 1555213536 : 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 20712489 : function get_subdt_fsd(nfsd, afsd_init, d_afsd) &
1055 3300189 : 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 70689774 : check_dt ! to compute subcycle timestep (s)
1070 :
1071 : integer (kind=int_kind) :: k
1072 :
1073 262225074 : check_dt(:) = bignum
1074 262225074 : do k = 1, nfsd
1075 241512585 : if (d_afsd(k) > puny) check_dt(k) = (1-afsd_init(k))/d_afsd(k)
1076 262225074 : if (d_afsd(k) < -puny) check_dt(k) = afsd_init(k)/ABS(d_afsd(k))
1077 : end do
1078 :
1079 262225074 : subdt = MINVAL(check_dt)
1080 :
1081 20712489 : end function get_subdt_fsd
1082 :
1083 :
1084 : !=======================================================================
1085 :
1086 : end module icepack_fsd
1087 :
1088 : !=======================================================================
1089 :
|
