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, ncat, nfsd
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 154 : subroutine icepack_init_wave(nfreq, &
73 154 : wave_spectrum_profile, & ! LCOV_EXCL_LINE
74 154 : 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 : wave_spectrum_data ! default values for nfreq profile
91 :
92 : ! set for 25 frequencies
93 :
94 154 : wave_spectrum_data = c0
95 :
96 : ! FOR TESTING ONLY - do not use for actual runs!!
97 154 : wave_spectrum_data(1) = 0.00015429197810590267
98 154 : wave_spectrum_data(2) = 0.002913531381636858
99 154 : wave_spectrum_data(3) = 0.02312942035496235
100 154 : wave_spectrum_data(4) = 0.07201970368623734
101 154 : wave_spectrum_data(5) = 0.06766948103904724
102 154 : wave_spectrum_data(6) = 0.005527883302420378
103 154 : wave_spectrum_data(7) = 3.326293881400488e-05
104 154 : wave_spectrum_data(8) = 6.815936703929992e-10
105 154 : wave_spectrum_data(9) = 2.419401186610744e-20
106 :
107 4004 : do k = 1, nfreq
108 4004 : 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 : 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 4004 : 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 4004 : dwavefreq(:) = wavefreq(:)*(SQRT(1.1_dbl_kind) - SQRT(c1/1.1_dbl_kind))
122 :
123 154 : 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 176435 : function get_dafsd_wave(afsd_init, fracture_hist, frac) &
132 176435 : result(d_afsd)
133 :
134 : real (kind=dbl_kind), dimension (:), intent(in) :: &
135 : afsd_init, fracture_hist
136 :
137 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
138 : frac
139 :
140 : ! output
141 : real (kind=dbl_kind), dimension (nfsd) :: &
142 : d_afsd
143 :
144 : ! local variables
145 : real (kind=dbl_kind), dimension (nfsd) :: &
146 352870 : loss, gain, omega
147 :
148 : integer (kind=int_kind) :: k
149 :
150 : character(len=*),parameter :: subname='(get_dafsd_wave)'
151 :
152 2293655 : do k = 1, nfsd
153 : ! fracture_hist is already normalized
154 13938365 : omega(k) = afsd_init(k)*SUM(fracture_hist(1:k-1))
155 : end do
156 :
157 2293655 : loss = omega
158 :
159 2293655 : do k =1,nfsd
160 27700295 : gain(k) = SUM(omega*frac(:,k))
161 : end do
162 :
163 2293655 : d_afsd(:) = gain(:) - loss(:)
164 :
165 2293655 : if (SUM(d_afsd(:)) > puny) then
166 0 : write(warnstr,*) subname, 'area not conserved, waves'
167 0 : call icepack_warnings_add(warnstr)
168 : endif
169 :
170 2293655 : WHERE (ABS(d_afsd).lt.puny) d_afsd = c0
171 :
172 176435 : end function get_dafsd_wave
173 :
174 : !=======================================================================
175 : !autodocument_start icepack_step_wavefracture
176 : !
177 : ! Given fracture histogram computed from local wave spectrum, evolve
178 : ! the floe size distribution
179 : !
180 : ! authors: 2018 Lettie Roach, NIWA/VUW
181 : !
182 22595640 : subroutine icepack_step_wavefracture(wave_spec_type, &
183 : dt, nfreq, & ! LCOV_EXCL_LINE
184 22595640 : aice, vice, aicen, & ! LCOV_EXCL_LINE
185 22595640 : wave_spectrum, wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
186 22595640 : trcrn, d_afsd_wave)
187 :
188 :
189 : character (len=char_len), intent(in) :: &
190 : wave_spec_type ! type of wave spectrum forcing
191 :
192 : integer (kind=int_kind), intent(in) :: &
193 : nfreq ! number of wave frequency categories
194 :
195 : real (kind=dbl_kind), intent(in) :: &
196 : dt, & ! time step ! LCOV_EXCL_LINE
197 : aice, & ! ice area fraction ! LCOV_EXCL_LINE
198 : vice ! ice volume per unit area
199 :
200 : real (kind=dbl_kind), dimension(ncat), intent(in) :: &
201 : aicen ! ice area fraction (categories)
202 :
203 : real (kind=dbl_kind), dimension (:), intent(in) :: &
204 : wavefreq, & ! wave frequencies (s^-1) ! LCOV_EXCL_LINE
205 : dwavefreq ! wave frequency bin widths (s^-1)
206 :
207 : real (kind=dbl_kind), dimension(:), intent(in) :: &
208 : wave_spectrum ! ocean surface wave spectrum as a function of frequency
209 : ! power spectral density of surface elevation, E(f) (units m^2 s)
210 :
211 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
212 : trcrn ! tracer array
213 :
214 : real (kind=dbl_kind), dimension(:), intent(out) :: &
215 : d_afsd_wave ! change in fsd due to waves
216 :
217 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
218 45191280 : d_afsdn_wave ! change in fsd due to waves, per category
219 :
220 : !autodocument_end
221 : ! local variables
222 : integer (kind=int_kind) :: &
223 : n, k, & ! LCOV_EXCL_LINE
224 : nsubt ! number of subcycles
225 :
226 : real (kind=dbl_kind), dimension (nfsd, nfsd) :: &
227 45191280 : frac
228 :
229 : real (kind=dbl_kind) :: &
230 : hbar , & ! mean ice thickness ! LCOV_EXCL_LINE
231 : elapsed_t , & ! elapsed subcycling time ! LCOV_EXCL_LINE
232 : subdt , & ! subcycling time step ! LCOV_EXCL_LINE
233 : cons_error ! area conservation error
234 :
235 : real (kind=dbl_kind), dimension (nfsd) :: &
236 45191280 : fracture_hist, & ! fracture histogram ! LCOV_EXCL_LINE
237 22595640 : afsd_init , & ! tracer array ! LCOV_EXCL_LINE
238 45191280 : afsd_tmp , & ! tracer array ! LCOV_EXCL_LINE
239 22595640 : d_afsd_tmp ! change
240 :
241 : real (kind=dbl_kind) :: &
242 : local_sig_ht
243 :
244 : character(len=*),parameter :: &
245 : subname='(icepack_step_wavefracture)'
246 :
247 : !------------------------------------
248 :
249 : ! initialize
250 293743320 : d_afsd_wave (:) = c0
251 1491312240 : d_afsdn_wave (:,:) = c0
252 293743320 : fracture_hist (:) = c0
253 :
254 : ! if all ice is not in first floe size category
255 22636476 : if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then
256 :
257 587274714 : local_sig_ht = c4*SQRT(SUM(wave_spectrum(:)*dwavefreq(:)))
258 : ! do not try to fracture for minimal ice concentration or zero wave spectrum
259 : ! if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then
260 22587489 : if ((aice > p01).and.(local_sig_ht>0.1_dbl_kind)) then
261 :
262 11417 : hbar = vice / aice
263 :
264 : ! calculate fracture histogram
265 : call wave_frac(nfreq, wave_spec_type, &
266 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
267 11417 : hbar, wave_spectrum, fracture_hist)
268 :
269 11417 : if (icepack_warnings_aborted(subname)) return
270 :
271 : ! if fracture occurs
272 159838 : if (MAXVAL(fracture_hist) > puny) then
273 : ! protect against small numerical errors
274 9329 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
275 9329 : if (icepack_warnings_aborted(subname)) return
276 :
277 55974 : do n = 1, ncat
278 :
279 606385 : afsd_init(:) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
280 :
281 : ! if there is ice, and a FSD, and not all ice is the smallest floe size
282 : if ((aicen(n) > puny) .and. (SUM(afsd_init(:)) > puny) &
283 615714 : .and. (afsd_init(1) < c1)) then
284 :
285 581113 : afsd_tmp = afsd_init
286 :
287 : ! frac does not vary within subcycle
288 7018057 : frac(:,:) = c0
289 536412 : do k = 2, nfsd
290 3486678 : frac(k,1:k-1) = fracture_hist(1:k-1)
291 : end do
292 581113 : do k = 1, nfsd
293 18218881 : if (SUM(frac(k,:)) > c0) frac(k,:) = frac(k,:)/SUM(frac(k,:))
294 : end do
295 :
296 : ! adaptive sub-timestep
297 44701 : elapsed_t = c0
298 44701 : cons_error = c0
299 44701 : nsubt = 0
300 221136 : DO WHILE (elapsed_t < dt)
301 209511 : nsubt = nsubt + 1
302 :
303 : ! if all floes in smallest category already, exit
304 209511 : if (afsd_tmp(1).ge.c1-puny) EXIT
305 :
306 : ! calculate d_afsd using current afstd
307 176435 : d_afsd_tmp = get_dafsd_wave(afsd_tmp, fracture_hist, frac)
308 :
309 : ! check in case wave fracture struggles to converge
310 176435 : if (nsubt>100) then
311 0 : write(warnstr,*) subname, &
312 0 : 'warning: step_wavefracture struggling to converge'
313 0 : call icepack_warnings_add(warnstr)
314 : endif
315 :
316 : ! required timestep
317 176435 : subdt = get_subdt_fsd(afsd_tmp, d_afsd_tmp)
318 176435 : subdt = MIN(subdt, dt)
319 :
320 : ! update afsd
321 2293655 : afsd_tmp = afsd_tmp + subdt * d_afsd_tmp(:)
322 :
323 : ! check conservation and negatives
324 2470090 : if (MINVAL(afsd_tmp) < -puny) then
325 0 : write(warnstr,*) subname, 'wb, <0 loop'
326 0 : call icepack_warnings_add(warnstr)
327 : endif
328 2470090 : if (MAXVAL(afsd_tmp) > c1+puny) then
329 0 : write(warnstr,*) subname, 'wb, >1 loop'
330 0 : call icepack_warnings_add(warnstr)
331 : endif
332 :
333 : ! update time
334 176435 : elapsed_t = elapsed_t + subdt
335 :
336 : END DO ! elapsed_t < dt
337 :
338 : ! In some cases---particularly for strong fracturing---the equation
339 : ! for wave fracture does not quite conserve area.
340 : ! With the dummy wave forcing, this happens < 2% of the time (in
341 : ! 1997) and is always less than 10^-7.
342 : ! Simply renormalizing may cause the first floe size
343 : ! category to reduce, which is not physically allowed
344 : ! to happen. So we adjust here
345 581113 : cons_error = SUM(afsd_tmp) - c1
346 :
347 : ! area loss: add to first category
348 44701 : if (cons_error.lt.c0) then
349 7536 : afsd_tmp(1) = afsd_tmp(1) - cons_error
350 : else
351 : ! area gain: take it from the largest possible category
352 393643 : do k = nfsd, 1, -1
353 393643 : if (afsd_tmp(k).gt.cons_error) then
354 37165 : afsd_tmp(k) = afsd_tmp(k) - cons_error
355 37165 : EXIT
356 : end if
357 : end do
358 : end if
359 :
360 : ! update trcrn
361 1117525 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp)
362 44701 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
363 44701 : if (icepack_warnings_aborted(subname)) return
364 :
365 : ! for diagnostics
366 581113 : d_afsdn_wave(:,n) = afsd_tmp(:) - afsd_init(:)
367 581113 : d_afsd_wave (:) = d_afsd_wave(:) + aicen(n)*d_afsdn_wave(:,n)
368 : endif ! aicen > puny
369 : enddo ! n
370 : endif ! fracture hist > 0
371 :
372 : endif ! aice > p01
373 : endif ! all small floes
374 :
375 : end subroutine icepack_step_wavefracture
376 :
377 : !=======================================================================
378 : !
379 : ! Calculates functions to describe the change in the FSD when waves
380 : ! fracture ice, given a wave spectrum (1D frequency, nfreq (default 25)
381 : ! frequency bins)
382 : !
383 : ! We calculate extrema and if these are successive maximum,
384 : ! minimum, maximum or vice versa, and have strain greater than a
385 : ! critical strain, break ice and create new floes with lengths equal
386 : ! to these distances. Based on MatLab code written by Chris Horvat,
387 : ! from Horvat & Tziperman (2015).
388 : !
389 : ! Note that a realization of sea surface height requires a random phase.
390 : !
391 : ! authors: 2018 Lettie Roach, NIWA/VUW
392 :
393 11417 : subroutine wave_frac(nfreq, wave_spec_type, &
394 22834 : wavefreq, dwavefreq, & ! LCOV_EXCL_LINE
395 11417 : hbar, spec_efreq, frac_local)
396 :
397 : integer (kind=int_kind), intent(in) :: &
398 : nfreq ! number of wave frequency categories
399 :
400 : character (len=char_len), intent(in) :: &
401 : wave_spec_type ! type of wave spectrum forcing
402 :
403 : real (kind=dbl_kind), intent(in) :: &
404 : hbar ! mean ice thickness (m)
405 :
406 : real (kind=dbl_kind), dimension (:), intent(in) :: &
407 : wavefreq, & ! wave frequencies (s^-1) ! LCOV_EXCL_LINE
408 : dwavefreq, & ! wave frequency bin widths (s^-1) ! LCOV_EXCL_LINE
409 : spec_efreq ! wave spectrum (m^2 s)
410 :
411 : real (kind=dbl_kind), dimension (nfsd), intent(out) :: &
412 : frac_local ! fracturing histogram
413 :
414 : ! local variables
415 :
416 : integer (kind=int_kind) :: j, k, iter, loop_max_iter
417 :
418 : real (kind=dbl_kind) :: &
419 : fracerror ! difference between successive histograms
420 :
421 : real (kind=dbl_kind), parameter :: &
422 : errortol = 6.5e-4 ! tolerance in error between successive histograms
423 :
424 : real (kind=dbl_kind), dimension(nfreq) :: &
425 22834 : lambda, & ! wavelengths (m) ! LCOV_EXCL_LINE
426 22834 : spec_coeff, & ! LCOV_EXCL_LINE
427 22834 : phi, rand_array, summand
428 :
429 : real (kind=dbl_kind), dimension(nx) :: &
430 : fraclengths
431 :
432 : real (kind=dbl_kind), dimension(nx) :: &
433 : X, & ! spatial domain (m) ! LCOV_EXCL_LINE
434 : eta ! sea surface height field (m)
435 :
436 : real (kind=dbl_kind), dimension(nfsd) :: &
437 11417 : frachistogram, & ! histogram ! LCOV_EXCL_LINE
438 11417 : prev_frac_local ! previous histogram
439 :
440 : character(len=*),parameter :: &
441 : subname='(wave_frac)'
442 :
443 :
444 11417 : if (trim(wave_spec_type).eq.'random') then
445 : ! run wave fracture to convergence
446 0 : loop_max_iter = max_no_iter
447 : else
448 11417 : loop_max_iter = 1
449 : end if
450 :
451 : ! spatial domain
452 114181417 : do j = 1, nx
453 114181417 : X(j)= j*dx
454 : end do
455 :
456 : ! dispersion relation
457 296842 : lambda (:) = gravit/(c2*pi*wavefreq (:)**2)
458 :
459 : ! spectral coefficients
460 296842 : spec_coeff = sqrt(c2*spec_efreq*dwavefreq)
461 :
462 : ! initialize frac lengths
463 11417 : fraclengths(:) = c0
464 148421 : prev_frac_local(:) = c0
465 148421 : frachistogram(:) = c0
466 11417 : fracerror = bignum
467 :
468 : ! loop while fracerror greater than error tolerance
469 11417 : iter = 0
470 22834 : do while (iter < loop_max_iter .and. fracerror > errortol)
471 11417 : iter = iter + 1
472 :
473 : ! Phase for each Fourier component may be constant or
474 : ! a random phase that varies in each i loop
475 : ! See documentation for discussion
476 11417 : if (trim(wave_spec_type).eq.'random') then
477 0 : call RANDOM_NUMBER(rand_array)
478 0 : if (icepack_warnings_aborted(subname)) return
479 : else
480 296842 : rand_array(:) = p5
481 : endif
482 296842 : phi = c2*pi*rand_array
483 :
484 114181417 : do j = 1, nx
485 : ! SSH field in space (sum over wavelengths, no attenuation)
486 2968420000 : summand = spec_coeff*COS(2*pi*X(j)/lambda+phi)
487 2968431417 : eta(j) = SUM(summand)
488 : end do
489 :
490 11417 : fraclengths(:) = c0
491 114181417 : if ((SUM(ABS(eta)) > puny).and.(hbar > puny)) then
492 11417 : call get_fraclengths(X, eta, fraclengths, hbar)
493 11417 : if (icepack_warnings_aborted(subname)) return
494 : end if
495 :
496 : ! convert from diameter to radii
497 114181417 : fraclengths(:) = fraclengths(:)/c2
498 :
499 20891417 : if (ALL(fraclengths.lt.floe_rad_l(1))) then
500 27144 : frac_local(:) = c0
501 : else
502 : ! bin into FS cats
503 : ! accumulate the frac histogram each iteration
504 93299329 : do j = 1, size(fraclengths)
505 93299329 : if (fraclengths(j).gt.floe_rad_l(1)) then
506 39743520 : do k = 1, nfsd-1
507 36431560 : if ((fraclengths(j) >= floe_rad_l(k)) .and. &
508 3311960 : (fraclengths(j) < floe_rad_l(k+1))) then
509 3311716 : frachistogram(k) = frachistogram(k) + 1
510 : end if
511 : end do
512 3311960 : if (fraclengths(j)>floe_rad_l(nfsd)) frachistogram(nfsd) = frachistogram(nfsd) + 1
513 : end if
514 : end do
515 :
516 121277 : do k = 1, nfsd
517 121277 : frac_local(k) = floe_rad_c(k)*frachistogram(k)
518 : end do
519 :
520 : ! normalize
521 345173 : if (SUM(frac_local) /= c0) frac_local(:) = frac_local(:) / SUM(frac_local(:))
522 :
523 : end if
524 :
525 : ! wave fracture run to convergence
526 11417 : if (trim(wave_spec_type).eq.'random') then
527 :
528 : ! check avg frac local against previous iteration
529 0 : fracerror = SUM(ABS(frac_local - prev_frac_local))/nfsd
530 :
531 : ! save histogram for next iteration
532 0 : prev_frac_local = frac_local
533 :
534 : end if
535 :
536 : END DO
537 :
538 11417 : if (iter >= max_no_iter) then
539 0 : write(warnstr,*) subname,'warning: wave_frac struggling to converge'
540 0 : call icepack_warnings_add(warnstr)
541 : endif
542 :
543 : end subroutine wave_frac
544 :
545 : !===========================================================================
546 : !
547 : ! Given the (attenuated) sea surface height, find the strain across triplets
548 : ! of max, min, max or min, max, min (local extrema within 10m).
549 : ! If this strain is greater than the critical strain, ice can fracture
550 : ! and new floes are formed with sizes equal to the distances between
551 : ! extrema. Based on MatLab code written by Chris Horvat,
552 : ! from Horvat & Tziperman (2015).
553 : !
554 : ! authors: 2016 Lettie Roach, NIWA/VUW
555 : !
556 11417 : subroutine get_fraclengths(X, eta, fraclengths, hbar)
557 :
558 : real (kind=dbl_kind), intent(in) :: &
559 : hbar ! mean thickness (m)
560 :
561 : real (kind=dbl_kind), intent(in), dimension (nx) :: &
562 : X, & ! spatial domain (m) ! LCOV_EXCL_LINE
563 : eta ! sea surface height field (m)
564 :
565 : real (kind=dbl_kind), intent(inout), dimension (nx) :: &
566 : fraclengths ! The distances between fracture points
567 : ! Size cannot be greater than nx.
568 : ! In practice, will be much less
569 :
570 : ! local variables
571 : integer (kind=int_kind) :: &
572 : spcing, & ! distance over which to search for extrema on each side of point ! LCOV_EXCL_LINE
573 : j, k, & ! indices to iterate over domain ! LCOV_EXCL_LINE
574 : first, last, & ! indices over which to search for extrema ! LCOV_EXCL_LINE
575 : j_neg, & ! nearest extrema backwards ! LCOV_EXCL_LINE
576 : j_pos, & ! nearest extrema forwards ! LCOV_EXCL_LINE
577 : n_above ! number of points where strain is above critical
578 :
579 : real (kind=dbl_kind), dimension(nx) :: &
580 : fracdistances, & ! distances in space where fracture has occurred ! LCOV_EXCL_LINE
581 : strain ! the strain between triplets of extrema
582 :
583 : logical (kind=log_kind), dimension(nx) :: &
584 : is_max, is_min,& ! arrays to hold whether each point is a local max or min ! LCOV_EXCL_LINE
585 : is_extremum, & ! or extremum ! LCOV_EXCL_LINE
586 : is_triplet ! or triplet of extrema
587 :
588 : real (kind=dbl_kind) :: &
589 : denominator, & ! denominator in strain equation ! LCOV_EXCL_LINE
590 : delta, & ! difference in x between current and prev extrema ! LCOV_EXCL_LINE
591 : delta_pos ! difference in x between next and current extrema
592 :
593 : integer (kind=int_kind), dimension(1) :: &
594 : maxj, minj ! indices of local max and min
595 :
596 : ! ------- equivalent of peakfinder2
597 : ! given eta and spcing, compute extremelocs in ascending order
598 11417 : spcing = nint(threshold/dx)
599 :
600 11417 : is_max = .false.
601 11417 : is_min = .false.
602 11417 : is_extremum = .false.
603 11417 : is_triplet = .false.
604 11417 : strain = c0
605 11417 : j_neg = 0
606 11417 : j_pos = 0
607 11417 : fraclengths(:) = c0
608 :
609 : ! search for local max and min within spacing
610 : ! on either side of each point
611 :
612 114181417 : do j = 1, nx
613 :
614 : ! indices within which to search for local max and min
615 114170000 : first = MAX(1,j-spcing)
616 114170000 : last = MIN(nx,j+spcing)
617 :
618 : ! location of max and min within spacing
619 114170000 : maxj = MAXLOC(eta(first:last))
620 114170000 : minj = MINLOC(eta(first:last))
621 :
622 : ! current j is the max or the min, save it
623 114170000 : if (maxj(1)+first-1 == j) is_max(j) = .true.
624 114170000 : if (minj(1)+first-1 == j) is_min(j) = .true.
625 :
626 : ! save whether max or min in one array
627 114181417 : if (is_min(j).or.is_max(j)) is_extremum(j) = .true.
628 : end do
629 :
630 : ! loop over points
631 : ! nothing can happen at the first or last
632 114158583 : do j = 2, nx-1
633 114158583 : if (is_extremum(j)) then
634 4076547 : if (j == 2) then
635 0 : if (is_extremum(1)) j_neg = 1
636 : else
637 113846393 : do k = j-1, 1, -1
638 113846393 : if (is_extremum(k)) then
639 4076547 : j_neg = k
640 4076547 : EXIT
641 : end if
642 : end do
643 : end if
644 :
645 112955371 : do k = j+1, nx
646 112955371 : if (is_extremum(k)) then
647 4070738 : j_pos = k
648 4070738 : EXIT
649 : end if
650 : end do
651 :
652 : ! find triplets of max and min
653 4076547 : if ((j_neg > 0).and.(j_pos > 0)) then
654 4076547 : if (is_max(j_neg).and.is_min(j).and.is_max(j_pos)) &
655 1760070 : is_triplet(j) = .true.
656 4076547 : if (is_min(j_neg).and.is_max(j).and.is_min(j_pos)) &
657 1763554 : is_triplet(j) = .true.
658 : end if
659 :
660 : ! calculate strain
661 4076547 : if (is_triplet(j)) then
662 :
663 : ! finite differences
664 3523624 : delta_pos = X(j_pos) - X(j )
665 3523624 : delta = X(j ) - X(j_neg)
666 :
667 : ! This equation differs from HT2015 by a factor 2 in numerator
668 : ! and eta(j_pos). This is the correct form of the equation.
669 :
670 3523624 : denominator = delta*delta_pos*(delta+delta_pos)
671 :
672 3523624 : if (denominator.ne.c0) &
673 : strain(j) = ABS(hbar*(eta(j_neg)* delta_pos & ! LCOV_EXCL_LINE
674 : - eta(j )*(delta_pos+delta) & ! LCOV_EXCL_LINE
675 : + eta(j_pos)* delta) & ! LCOV_EXCL_LINE
676 3523624 : / denominator)
677 :
678 : end if ! is triplet
679 : end if ! is extremum
680 :
681 : end do ! loop over j
682 :
683 114181417 : n_above = COUNT(strain > straincrit)
684 11417 : fracdistances(:) = c0
685 :
686 : ! only do if we have some strains exceeding strain crit
687 11417 : if (n_above>0) then
688 :
689 9334 : k = 0
690 93349334 : do j = 1, nx
691 93349334 : if (strain(j) > straincrit) then
692 3321294 : k = k + 1
693 3321294 : fracdistances(k) = X(j)
694 : end if
695 : end do
696 :
697 3321294 : do j = 1, n_above-1
698 3321294 : fraclengths(j) = fracdistances(j+1) - fracdistances(j)
699 : end do
700 :
701 :
702 : end if ! n_above
703 :
704 11417 : end subroutine get_fraclengths
705 :
706 : !=======================================================================
707 :
708 : end module icepack_wavefracspec
709 :
710 : !=======================================================================
711 :
712 :
|
