Line data Source code
1 : ! This module contains the subroutines required to fracture sea ice
2 : ! by ocean surface waves
3 : !
4 : ! Theory based on:
5 : !
6 : ! Horvat, C., & Tziperman, E. (2015). A prognostic model of the sea-ice
7 : ! floe size and thickness distribution. The Cryosphere, 9(6), 2119–2134.
8 : ! doi:10.5194/tc-9-2119-2015
9 : !
10 : ! and implementation described in:
11 : !
12 : ! Roach, L. A., Horvat, C., Dean, S. M., & Bitz, C. M. (2018). An emergent
13 : ! sea ice floe size distribution in a global coupled ocean--sea ice model.
14 : ! Journal of Geophysical Research: Oceans, 123(6), 4322–4337.
15 : ! doi:10.1029/2017JC013692
16 : !
17 : ! now with some modifications to allow direct input of ocean surface wave spectrum.
18 : !
19 : ! We calculate the fractures that would occur if waves enter a fully ice-covered
20 : ! region defined in one dimension in the direction of propagation, and then apply
21 : ! the outcome proportionally to the ice-covered fraction in each grid cell. Assuming
22 : ! that sea ice flexes with the sea surface height field, strains are computed on this
23 : ! sub-grid-scale 1D domain. If the strain between successive extrema exceeds a critical
24 : ! value new floes are formed with diameters equal to the distance between the extrema.
25 : !
26 : ! authors: 2016-8 Lettie Roach, NIWA/VUW
27 : !
28 : !
29 : module icepack_wavefracspec
30 :
31 : use icepack_kinds
32 : use icepack_parameters, only: p01, p5, c0, c1, c2, c3, c4, c10
33 : use icepack_parameters, only: bignum, puny, gravit, pi
34 : use icepack_tracers, only: nt_fsd
35 : use icepack_warnings, only: warnstr, icepack_warnings_add, icepack_warnings_aborted
36 : use icepack_fsd
37 :
38 : implicit none
39 : private
40 : public :: icepack_init_wave, icepack_step_wavefracture
41 :
42 : real (kind=dbl_kind), parameter :: &
43 : swh_minval = 0.01_dbl_kind, & ! minimum value of wave height (m) ! LCOV_EXCL_LINE
44 : straincrit = 3.e-5_dbl_kind, & ! critical strain ! LCOV_EXCL_LINE
45 : D = 1.e4_dbl_kind, & ! domain size ! LCOV_EXCL_LINE
46 : dx = c1, & ! domain spacing ! LCOV_EXCL_LINE
47 : threshold = c10 ! peak-finding threshold -
48 : ! points are defined to be extrema if they
49 : ! are a local max or min over a distance
50 : ! of 10m on both sides, based on the
51 : ! observations of Toyota et al. (2011) who
52 : ! find this to be the order of the smallest
53 : ! floe size affected by wave fracture
54 :
55 : integer (kind=int_kind), parameter :: &
56 : nx = 10000 ! number of points in domain
57 :
58 : integer (kind=int_kind), parameter :: &
59 : max_no_iter = 100 ! max no of iterations to compute wave fracture
60 :
61 :
62 : !=======================================================================
63 :
64 : contains
65 :
66 : !=======================================================================
67 : !autodocument_start icepack_init_wave
68 : ! Initialize the wave spectrum and frequencies for the FSD
69 : !
70 : ! authors: 2018 Lettie Roach, NIWA/VUW
71 :
72 0 : subroutine icepack_init_wave(nfreq, &
73 : wave_spectrum_profile, & ! LCOV_EXCL_LINE
74 0 : wavefreq, dwavefreq)
75 :
76 : integer(kind=int_kind), intent(in) :: &
77 : nfreq ! number of wave frequencies
78 :
79 : real(kind=dbl_kind), dimension(:), intent(out) :: &
80 : wave_spectrum_profile, & ! ocean surface wave spectrum as a function of frequency ! LCOV_EXCL_LINE
81 : ! power spectral density of surface elevation, E(f) (units m^2 s)
82 : wavefreq, & ! wave frequencies (s^-1)
83 : dwavefreq ! wave frequency bin widths (s^-1)
84 :
85 : !autodocument_end
86 : ! local variables
87 : integer (kind=int_kind) :: k
88 :
89 : real(kind=dbl_kind), dimension(100) :: &
90 0 : wave_spectrum_data ! default values for nfreq profile
91 :
92 : ! set for 25 frequencies
93 :
94 0 : wave_spectrum_data = c0
95 :
96 : ! FOR TESTING ONLY - do not use for actual runs!!
97 0 : wave_spectrum_data(1) = 0.00015429197810590267
98 0 : wave_spectrum_data(2) = 0.002913531381636858
99 0 : wave_spectrum_data(3) = 0.02312942035496235
100 0 : wave_spectrum_data(4) = 0.07201970368623734
101 0 : wave_spectrum_data(5) = 0.06766948103904724
102 0 : wave_spectrum_data(6) = 0.005527883302420378
103 0 : wave_spectrum_data(7) = 3.326293881400488e-05
104 0 : wave_spectrum_data(8) = 6.815936703929992e-10
105 0 : wave_spectrum_data(9) = 2.419401186610744e-20
106 :
107 0 : do k = 1, nfreq
108 0 : wave_spectrum_profile(k) = wave_spectrum_data(k)
109 : enddo
110 :
111 : ! hardwired for wave coupling with NIWA version of Wavewatch
112 : ! From Wavewatch, f(n+1) = C*f(n) where C is a constant set by the user
113 : ! These freq are for C = 1.1
114 0 : wavefreq = (/0.04118, 0.045298, 0.0498278, 0.05481058, 0.06029164, &
115 : 0.06632081, 0.07295289, 0.08024818, 0.08827299, 0.09710029, & ! LCOV_EXCL_LINE
116 : 0.10681032, 0.11749136, 0.1292405, 0.14216454, 0.15638101, & ! LCOV_EXCL_LINE
117 : 0.17201911, 0.18922101, 0.20814312, 0.22895744, 0.25185317, & ! LCOV_EXCL_LINE
118 0 : 0.27703848, 0.30474234, 0.33521661, 0.36873826, 0.40561208/)
119 :
120 : ! boundaries of bin n are at f(n)*sqrt(1/C) and f(n)*sqrt(C)
121 0 : dwavefreq(:) = wavefreq(:)*(SQRT(1.1_dbl_kind) - SQRT(c1/1.1_dbl_kind))
122 :
123 0 : end subroutine icepack_init_wave
124 :
125 : !=======================================================================
126 : !
127 : ! Calculate the change in the FSD arising from wave fracture
128 : !
129 : ! authors: 2017 Lettie Roach, NIWA/VUW
130 : !
131 0 : function get_dafsd_wave(nfsd, afsd_init, fracture_hist, frac) &
132 0 : result(d_afsd)
133 :
134 : integer (kind=int_kind), intent(in) :: &
135 : nfsd ! number of floe size categories
136 :
137 : real (kind=dbl_kind), dimension (:), intent(in) :: &
138 : afsd_init, fracture_hist
139 :
140 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
141 : frac
142 :
143 : ! output
144 : real (kind=dbl_kind), dimension (nfsd) :: &
145 : d_afsd
146 :
147 : ! local variables
148 : real (kind=dbl_kind), dimension (nfsd) :: &
149 0 : loss, gain, omega
150 :
151 : integer (kind=int_kind) :: k
152 :
153 : character(len=*),parameter :: subname='(get_dafsd_wave)'
154 :
155 0 : do k = 1, nfsd
156 : ! fracture_hist is already normalized
157 0 : omega(k) = afsd_init(k)*SUM(fracture_hist(1:k-1))
158 : end do
159 :
160 0 : loss = omega
161 :
162 0 : do k =1,nfsd
163 0 : gain(k) = SUM(omega*frac(:,k))
164 : end do
165 :
166 0 : d_afsd(:) = gain(:) - loss(:)
167 :
168 0 : if (SUM(d_afsd(:)) > puny) then
169 0 : write(warnstr,*) subname, 'area not conserved, waves'
170 0 : call icepack_warnings_add(warnstr)
171 : endif
172 :
173 0 : WHERE (ABS(d_afsd).lt.puny) d_afsd = c0
174 :
175 0 : end function get_dafsd_wave
176 :
177 : !=======================================================================
178 : !autodocument_start icepack_step_wavefracture
179 : !
180 : ! Given fracture histogram computed from local wave spectrum, evolve
181 : ! the floe size distribution
182 : !
183 : ! authors: 2018 Lettie Roach, NIWA/VUW
184 : !
185 0 : subroutine icepack_step_wavefracture(wave_spec_type, &
186 : dt, ncat, nfsd, & ! LCOV_EXCL_LINE
187 : nfreq, & ! LCOV_EXCL_LINE
188 : aice, vice, aicen, & ! LCOV_EXCL_LINE
189 : floe_rad_l, floe_rad_c, & ! LCOV_EXCL_LINE
190 : wave_spectrum, wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
191 0 : trcrn, d_afsd_wave)
192 :
193 :
194 : character (len=char_len), intent(in) :: &
195 : wave_spec_type ! type of wave spectrum forcing
196 :
197 : integer (kind=int_kind), intent(in) :: &
198 : nfreq, & ! number of wave frequency categories ! LCOV_EXCL_LINE
199 : ncat, & ! number of thickness categories ! LCOV_EXCL_LINE
200 : nfsd ! number of floe size categories
201 :
202 : real (kind=dbl_kind), intent(in) :: &
203 : dt, & ! time step ! LCOV_EXCL_LINE
204 : aice, & ! ice area fraction ! LCOV_EXCL_LINE
205 : vice ! ice volume per unit area
206 :
207 : real (kind=dbl_kind), dimension(ncat), intent(in) :: &
208 : aicen ! ice area fraction (categories)
209 :
210 : real(kind=dbl_kind), dimension(:), intent(in) :: &
211 : floe_rad_l, & ! fsd size lower bound in m (radius) ! LCOV_EXCL_LINE
212 : floe_rad_c ! fsd size bin centre in m (radius)
213 :
214 : real (kind=dbl_kind), dimension (:), intent(in) :: &
215 : wavefreq, & ! wave frequencies (s^-1) ! LCOV_EXCL_LINE
216 : dwavefreq ! wave frequency bin widths (s^-1)
217 :
218 : real (kind=dbl_kind), dimension(:), intent(in) :: &
219 : wave_spectrum ! ocean surface wave spectrum as a function of frequency
220 : ! power spectral density of surface elevation, E(f) (units m^2 s)
221 :
222 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
223 : trcrn ! tracer array
224 :
225 : real (kind=dbl_kind), dimension(:), intent(out) :: &
226 : d_afsd_wave ! change in fsd due to waves
227 :
228 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
229 0 : d_afsdn_wave ! change in fsd due to waves, per category
230 :
231 : !autodocument_end
232 : ! local variables
233 : integer (kind=int_kind) :: &
234 : n, k, & ! LCOV_EXCL_LINE
235 : nsubt ! number of subcycles
236 :
237 : real (kind=dbl_kind), dimension (nfsd, nfsd) :: &
238 0 : frac
239 :
240 : real (kind=dbl_kind) :: &
241 : hbar , & ! mean ice thickness ! LCOV_EXCL_LINE
242 : elapsed_t , & ! elapsed subcycling time ! LCOV_EXCL_LINE
243 : subdt , & ! subcycling time step ! LCOV_EXCL_LINE
244 0 : cons_error ! area conservation error
245 :
246 : real (kind=dbl_kind), dimension (nfsd) :: &
247 : fracture_hist, & ! fracture histogram ! LCOV_EXCL_LINE
248 : afsd_init , & ! tracer array ! LCOV_EXCL_LINE
249 : afsd_tmp , & ! tracer array ! LCOV_EXCL_LINE
250 0 : d_afsd_tmp ! change
251 :
252 : character(len=*),parameter :: &
253 : subname='(icepack_step_wavefracture)'
254 :
255 : !------------------------------------
256 :
257 : ! initialize
258 0 : d_afsd_wave (:) = c0
259 0 : d_afsdn_wave (:,:) = c0
260 0 : fracture_hist (:) = c0
261 :
262 : ! if all ice is not in first floe size category
263 0 : if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then
264 :
265 :
266 : ! do not try to fracture for minimal ice concentration or zero wave spectrum
267 0 : if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then
268 :
269 0 : hbar = vice / aice
270 :
271 : ! calculate fracture histogram
272 : call wave_frac(nfsd, nfreq, wave_spec_type, &
273 : floe_rad_l, floe_rad_c, & ! LCOV_EXCL_LINE
274 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
275 0 : hbar, wave_spectrum, fracture_hist)
276 :
277 0 : if (icepack_warnings_aborted(subname)) return
278 :
279 : ! if fracture occurs
280 0 : if (MAXVAL(fracture_hist) > puny) then
281 : ! protect against small numerical errors
282 0 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
283 0 : if (icepack_warnings_aborted(subname)) return
284 :
285 0 : do n = 1, ncat
286 :
287 0 : afsd_init(:) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
288 :
289 : ! if there is ice, and a FSD, and not all ice is the smallest floe size
290 0 : if ((aicen(n) > puny) .and. (SUM(afsd_init(:)) > puny) &
291 0 : .and. (afsd_init(1) < c1)) then
292 :
293 0 : afsd_tmp = afsd_init
294 :
295 : ! frac does not vary within subcycle
296 0 : frac(:,:) = c0
297 0 : do k = 2, nfsd
298 0 : frac(k,1:k-1) = fracture_hist(1:k-1)
299 : end do
300 0 : do k = 1, nfsd
301 0 : if (SUM(frac(k,:)) > c0) frac(k,:) = frac(k,:)/SUM(frac(k,:))
302 : end do
303 :
304 : ! adaptive sub-timestep
305 0 : elapsed_t = c0
306 0 : cons_error = c0
307 0 : nsubt = 0
308 0 : DO WHILE (elapsed_t < dt)
309 0 : nsubt = nsubt + 1
310 :
311 : ! if all floes in smallest category already, exit
312 0 : if (afsd_tmp(1).ge.c1-puny) EXIT
313 :
314 : ! calculate d_afsd using current afstd
315 0 : d_afsd_tmp = get_dafsd_wave(nfsd, afsd_tmp, fracture_hist, frac)
316 :
317 : ! check in case wave fracture struggles to converge
318 0 : if (nsubt>100) then
319 0 : write(warnstr,*) subname, &
320 0 : 'warning: step_wavefracture struggling to converge'
321 0 : call icepack_warnings_add(warnstr)
322 : endif
323 :
324 : ! required timestep
325 0 : subdt = get_subdt_fsd(nfsd, afsd_tmp, d_afsd_tmp)
326 0 : subdt = MIN(subdt, dt)
327 :
328 : ! update afsd
329 0 : afsd_tmp = afsd_tmp + subdt * d_afsd_tmp(:)
330 :
331 : ! check conservation and negatives
332 0 : if (MINVAL(afsd_tmp) < -puny) then
333 0 : write(warnstr,*) subname, 'wb, <0 loop'
334 0 : call icepack_warnings_add(warnstr)
335 : endif
336 0 : if (MAXVAL(afsd_tmp) > c1+puny) then
337 0 : write(warnstr,*) subname, 'wb, >1 loop'
338 0 : call icepack_warnings_add(warnstr)
339 : endif
340 :
341 : ! update time
342 0 : elapsed_t = elapsed_t + subdt
343 :
344 : END DO ! elapsed_t < dt
345 :
346 : ! In some cases---particularly for strong fracturing---the equation
347 : ! for wave fracture does not quite conserve area.
348 : ! With the dummy wave forcing, this happens < 2% of the time (in
349 : ! 1997) and is always less than 10^-7.
350 : ! Simply renormalizing may cause the first floe size
351 : ! category to reduce, which is not physically allowed
352 : ! to happen. So we adjust here
353 0 : cons_error = SUM(afsd_tmp) - c1
354 :
355 : ! area loss: add to first category
356 0 : if (cons_error.lt.c0) then
357 0 : afsd_tmp(1) = afsd_tmp(1) - cons_error
358 : else
359 : ! area gain: take it from the largest possible category
360 0 : do k = nfsd, 1, -1
361 0 : if (afsd_tmp(k).gt.cons_error) then
362 0 : afsd_tmp(k) = afsd_tmp(k) - cons_error
363 0 : EXIT
364 : end if
365 : end do
366 : end if
367 :
368 : ! update trcrn
369 0 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp)
370 0 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
371 0 : if (icepack_warnings_aborted(subname)) return
372 :
373 : ! for diagnostics
374 0 : d_afsdn_wave(:,n) = afsd_tmp(:) - afsd_init(:)
375 0 : d_afsd_wave (:) = d_afsd_wave(:) + aicen(n)*d_afsdn_wave(:,n)
376 : endif ! aicen > puny
377 : enddo ! n
378 : endif ! fracture hist > 0
379 :
380 : endif ! aice > p01
381 : endif ! all small floes
382 :
383 : end subroutine icepack_step_wavefracture
384 :
385 : !=======================================================================
386 : !
387 : ! Calculates functions to describe the change in the FSD when waves
388 : ! fracture ice, given a wave spectrum (1D frequency, nfreq (default 25)
389 : ! frequency bins)
390 : !
391 : ! We calculate extrema and if these are successive maximum,
392 : ! minimum, maximum or vice versa, and have strain greater than a
393 : ! critical strain, break ice and create new floes with lengths equal
394 : ! to these distances. Based on MatLab code written by Chris Horvat,
395 : ! from Horvat & Tziperman (2015).
396 : !
397 : ! Note that a realization of sea surface height requires a random phase.
398 : !
399 : ! authors: 2018 Lettie Roach, NIWA/VUW
400 :
401 0 : subroutine wave_frac(nfsd, nfreq, wave_spec_type, &
402 : floe_rad_l, floe_rad_c, & ! LCOV_EXCL_LINE
403 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
404 0 : hbar, spec_efreq, frac_local)
405 :
406 : integer (kind=int_kind), intent(in) :: &
407 : nfsd, & ! number of floe size categories ! LCOV_EXCL_LINE
408 : nfreq ! number of wave frequency categories
409 :
410 : character (len=char_len), intent(in) :: &
411 : wave_spec_type ! type of wave spectrum forcing
412 :
413 : real (kind=dbl_kind), intent(in) :: &
414 : hbar ! mean ice thickness (m)
415 :
416 : real(kind=dbl_kind), dimension(:), intent(in) :: &
417 : floe_rad_l, & ! fsd size lower bound in m (radius) ! LCOV_EXCL_LINE
418 : floe_rad_c ! fsd size bin centre in m (radius)
419 :
420 : real (kind=dbl_kind), dimension (:), intent(in) :: &
421 : wavefreq, & ! wave frequencies (s^-1) ! LCOV_EXCL_LINE
422 : dwavefreq, & ! wave frequency bin widths (s^-1) ! LCOV_EXCL_LINE
423 : spec_efreq ! wave spectrum (m^2 s)
424 :
425 : real (kind=dbl_kind), dimension (nfsd), intent(out) :: &
426 : frac_local ! fracturing histogram
427 :
428 : ! local variables
429 :
430 : integer (kind=int_kind) :: j, k, iter, loop_max_iter
431 :
432 : real (kind=dbl_kind) :: &
433 0 : fracerror ! difference between successive histograms
434 :
435 : real (kind=dbl_kind), parameter :: &
436 : errortol = 6.5e-4 ! tolerance in error between successive histograms
437 :
438 : real (kind=dbl_kind), dimension(nfreq) :: &
439 : lambda, & ! wavelengths (m) ! LCOV_EXCL_LINE
440 : spec_coeff, & ! LCOV_EXCL_LINE
441 0 : phi, rand_array, summand
442 :
443 : real (kind=dbl_kind), dimension(nx) :: &
444 0 : fraclengths
445 :
446 : real (kind=dbl_kind), dimension(nx) :: &
447 : X, & ! spatial domain (m) ! LCOV_EXCL_LINE
448 0 : eta ! sea surface height field (m)
449 :
450 : real (kind=dbl_kind), dimension(nfsd) :: &
451 : frachistogram, & ! histogram ! LCOV_EXCL_LINE
452 0 : prev_frac_local ! previous histogram
453 :
454 : character(len=*),parameter :: &
455 : subname='(wave_frac)'
456 :
457 :
458 0 : if (trim(wave_spec_type).eq.'random') then
459 : ! run wave fracture to convergence
460 0 : loop_max_iter = max_no_iter
461 : else
462 0 : loop_max_iter = 1
463 : end if
464 :
465 : ! spatial domain
466 0 : do j = 1, nx
467 0 : X(j)= j*dx
468 : end do
469 :
470 : ! dispersion relation
471 0 : lambda (:) = gravit/(c2*pi*wavefreq (:)**2)
472 :
473 : ! spectral coefficients
474 0 : spec_coeff = sqrt(c2*spec_efreq*dwavefreq)
475 :
476 : ! initialize frac lengths
477 0 : fraclengths(:) = c0
478 0 : prev_frac_local(:) = c0
479 0 : frachistogram(:) = c0
480 0 : fracerror = bignum
481 :
482 : ! loop while fracerror greater than error tolerance
483 0 : iter = 0
484 0 : do while (iter < loop_max_iter .and. fracerror > errortol)
485 0 : iter = iter + 1
486 :
487 : ! Phase for each Fourier component may be constant or
488 : ! a random phase that varies in each i loop
489 : ! See documentation for discussion
490 0 : if (trim(wave_spec_type).eq.'random') then
491 0 : call RANDOM_NUMBER(rand_array)
492 0 : if (icepack_warnings_aborted(subname)) return
493 : else
494 0 : rand_array(:) = p5
495 : endif
496 0 : phi = c2*pi*rand_array
497 :
498 0 : do j = 1, nx
499 : ! SSH field in space (sum over wavelengths, no attenuation)
500 0 : summand = spec_coeff*COS(2*pi*X(j)/lambda+phi)
501 0 : eta(j) = SUM(summand)
502 : end do
503 :
504 0 : fraclengths(:) = c0
505 0 : if ((SUM(ABS(eta)) > puny).and.(hbar > puny)) then
506 0 : call get_fraclengths(X, eta, fraclengths, hbar)
507 0 : if (icepack_warnings_aborted(subname)) return
508 : end if
509 :
510 : ! convert from diameter to radii
511 0 : fraclengths(:) = fraclengths(:)/c2
512 :
513 0 : if (ALL(fraclengths.lt.floe_rad_l(1))) then
514 0 : frac_local(:) = c0
515 : else
516 : ! bin into FS cats
517 : ! accumulate the frac histogram each iteration
518 0 : do j = 1, size(fraclengths)
519 0 : if (fraclengths(j).gt.floe_rad_l(1)) then
520 0 : do k = 1, nfsd-1
521 0 : if ((fraclengths(j) >= floe_rad_l(k)) .and. &
522 0 : (fraclengths(j) < floe_rad_l(k+1))) then
523 0 : frachistogram(k) = frachistogram(k) + 1
524 : end if
525 : end do
526 0 : if (fraclengths(j)>floe_rad_l(nfsd)) frachistogram(nfsd) = frachistogram(nfsd) + 1
527 : end if
528 : end do
529 :
530 0 : do k = 1, nfsd
531 0 : frac_local(k) = floe_rad_c(k)*frachistogram(k)
532 : end do
533 :
534 : ! normalize
535 0 : if (SUM(frac_local) /= c0) frac_local(:) = frac_local(:) / SUM(frac_local(:))
536 :
537 : end if
538 :
539 : ! wave fracture run to convergence
540 0 : if (trim(wave_spec_type).eq.'random') then
541 :
542 : ! check avg frac local against previous iteration
543 0 : fracerror = SUM(ABS(frac_local - prev_frac_local))/nfsd
544 :
545 : ! save histogram for next iteration
546 0 : prev_frac_local = frac_local
547 :
548 : end if
549 :
550 : END DO
551 :
552 0 : if (iter >= max_no_iter) then
553 0 : write(warnstr,*) subname,'warning: wave_frac struggling to converge'
554 0 : call icepack_warnings_add(warnstr)
555 : endif
556 :
557 : end subroutine wave_frac
558 :
559 : !===========================================================================
560 : !
561 : ! Given the (attenuated) sea surface height, find the strain across triplets
562 : ! of max, min, max or min, max, min (local extrema within 10m).
563 : ! If this strain is greater than the critical strain, ice can fracture
564 : ! and new floes are formed with sizes equal to the distances between
565 : ! extrema. Based on MatLab code written by Chris Horvat,
566 : ! from Horvat & Tziperman (2015).
567 : !
568 : ! authors: 2016 Lettie Roach, NIWA/VUW
569 : !
570 0 : subroutine get_fraclengths(X, eta, fraclengths, hbar)
571 :
572 : real (kind=dbl_kind), intent(in) :: &
573 : hbar ! mean thickness (m)
574 :
575 : real (kind=dbl_kind), intent(in), dimension (nx) :: &
576 : X, & ! spatial domain (m) ! LCOV_EXCL_LINE
577 : eta ! sea surface height field (m)
578 :
579 : real (kind=dbl_kind), intent(inout), dimension (nx) :: &
580 : fraclengths ! The distances between fracture points
581 : ! Size cannot be greater than nx.
582 : ! In practice, will be much less
583 :
584 : ! local variables
585 : integer (kind=int_kind) :: &
586 : spcing, & ! distance over which to search for extrema on each side of point ! LCOV_EXCL_LINE
587 : j, k, & ! indices to iterate over domain ! LCOV_EXCL_LINE
588 : first, last, & ! indices over which to search for extrema ! LCOV_EXCL_LINE
589 : j_neg, & ! nearest extrema backwards ! LCOV_EXCL_LINE
590 : j_pos, & ! nearest extrema forwards ! LCOV_EXCL_LINE
591 : n_above ! number of points where strain is above critical
592 :
593 : real (kind=dbl_kind), dimension(nx) :: &
594 : fracdistances, & ! distances in space where fracture has occurred ! LCOV_EXCL_LINE
595 0 : strain ! the strain between triplets of extrema
596 :
597 : logical (kind=log_kind), dimension(nx) :: &
598 : is_max, is_min,& ! arrays to hold whether each point is a local max or min ! LCOV_EXCL_LINE
599 : is_extremum, & ! or extremum ! LCOV_EXCL_LINE
600 : is_triplet ! or triplet of extrema
601 :
602 : real (kind=dbl_kind) :: &
603 : denominator, & ! denominator in strain equation ! LCOV_EXCL_LINE
604 : delta, & ! difference in x between current and prev extrema ! LCOV_EXCL_LINE
605 0 : delta_pos ! difference in x between next and current extrema
606 :
607 : integer (kind=int_kind), dimension(1) :: &
608 : maxj, minj ! indices of local max and min
609 :
610 : ! ------- equivalent of peakfinder2
611 : ! given eta and spcing, compute extremelocs in ascending order
612 0 : spcing = nint(threshold/dx)
613 :
614 0 : is_max = .false.
615 0 : is_min = .false.
616 0 : is_extremum = .false.
617 0 : is_triplet = .false.
618 0 : strain = c0
619 0 : j_neg = 0
620 0 : j_pos = 0
621 0 : fraclengths(:) = c0
622 :
623 : ! search for local max and min within spacing
624 : ! on either side of each point
625 :
626 0 : do j = 1, nx
627 :
628 : ! indices within which to search for local max and min
629 0 : first = MAX(1,j-spcing)
630 0 : last = MIN(nx,j+spcing)
631 :
632 : ! location of max and min within spacing
633 0 : maxj = MAXLOC(eta(first:last))
634 0 : minj = MINLOC(eta(first:last))
635 :
636 : ! current j is the max or the min, save it
637 0 : if (maxj(1)+first-1 == j) is_max(j) = .true.
638 0 : if (minj(1)+first-1 == j) is_min(j) = .true.
639 :
640 : ! save whether max or min in one array
641 0 : if (is_min(j).or.is_max(j)) is_extremum(j) = .true.
642 : end do
643 :
644 : ! loop over points
645 : ! nothing can happen at the first or last
646 0 : do j = 2, nx-1
647 0 : if (is_extremum(j)) then
648 0 : if (j == 2) then
649 0 : if (is_extremum(1)) j_neg = 1
650 : else
651 0 : do k = j-1, 1, -1
652 0 : if (is_extremum(k)) then
653 0 : j_neg = k
654 0 : EXIT
655 : end if
656 : end do
657 : end if
658 :
659 0 : do k = j+1, nx
660 0 : if (is_extremum(k)) then
661 0 : j_pos = k
662 0 : EXIT
663 : end if
664 : end do
665 :
666 : ! find triplets of max and min
667 0 : if ((j_neg > 0).and.(j_pos > 0)) then
668 0 : if (is_max(j_neg).and.is_min(j).and.is_max(j_pos)) &
669 0 : is_triplet(j) = .true.
670 0 : if (is_min(j_neg).and.is_max(j).and.is_min(j_pos)) &
671 0 : is_triplet(j) = .true.
672 : end if
673 :
674 : ! calculate strain
675 0 : if (is_triplet(j)) then
676 :
677 : ! finite differences
678 0 : delta_pos = X(j_pos) - X(j )
679 0 : delta = X(j ) - X(j_neg)
680 :
681 : ! This equation differs from HT2015 by a factor 2 in numerator
682 : ! and eta(j_pos). This is the correct form of the equation.
683 :
684 0 : denominator = delta*delta_pos*(delta+delta_pos)
685 :
686 0 : if (denominator.ne.c0) &
687 : strain(j) = ABS(hbar*(eta(j_neg)* delta_pos & ! LCOV_EXCL_LINE
688 : - eta(j )*(delta_pos+delta) & ! LCOV_EXCL_LINE
689 : + eta(j_pos)* delta) & ! LCOV_EXCL_LINE
690 0 : / denominator)
691 :
692 : end if ! is triplet
693 : end if ! is extremum
694 :
695 : end do ! loop over j
696 :
697 0 : n_above = COUNT(strain > straincrit)
698 0 : fracdistances(:) = c0
699 :
700 : ! only do if we have some strains exceeding strain crit
701 0 : if (n_above>0) then
702 :
703 0 : k = 0
704 0 : do j = 1, nx
705 0 : if (strain(j) > straincrit) then
706 0 : k = k + 1
707 0 : fracdistances(k) = X(j)
708 : end if
709 : end do
710 :
711 0 : do j = 1, n_above-1
712 0 : fraclengths(j) = fracdistances(j+1) - fracdistances(j)
713 : end do
714 :
715 :
716 : end if ! n_above
717 :
718 0 : end subroutine get_fraclengths
719 :
720 : !=======================================================================
721 :
722 : end module icepack_wavefracspec
723 :
724 : !=======================================================================
725 :
726 :
|
