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