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)
44 : straincrit = 3.e-5_dbl_kind, & ! critical strain
45 : D = 1.e4_dbl_kind, & ! domain size
46 : dx = c1, & ! domain spacing
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 24132 : subroutine icepack_init_wave(nfreq, &
73 24132 : wave_spectrum_profile, &
74 24132 : 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
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 24132 : wave_spectrum_data = c0
95 :
96 : ! FOR TESTING ONLY - do not use for actual runs!!
97 24132 : wave_spectrum_data(1) = 0.00015429197810590267
98 24132 : wave_spectrum_data(2) = 0.002913531381636858
99 24132 : wave_spectrum_data(3) = 0.02312942035496235
100 24132 : wave_spectrum_data(4) = 0.07201970368623734
101 24132 : wave_spectrum_data(5) = 0.06766948103904724
102 24132 : wave_spectrum_data(6) = 0.005527883302420378
103 24132 : wave_spectrum_data(7) = 3.326293881400488e-05
104 24132 : wave_spectrum_data(8) = 6.815936703929992e-10
105 24132 : wave_spectrum_data(9) = 2.419401186610744e-20
106 :
107 627432 : do k = 1, nfreq
108 627432 : 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, &
116 : 0.10681032, 0.11749136, 0.1292405, 0.14216454, 0.15638101, &
117 : 0.17201911, 0.18922101, 0.20814312, 0.22895744, 0.25185317, &
118 627432 : 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 627432 : dwavefreq(:) = wavefreq(:)*(SQRT(1.1_dbl_kind) - SQRT(c1/1.1_dbl_kind))
122 :
123 24132 : 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(afsd_init, fracture_hist, frac) &
132 0 : 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 0 : loss, gain, omega
147 :
148 : integer (kind=int_kind) :: k
149 :
150 : character(len=*),parameter :: subname='(get_dafsd_wave)'
151 :
152 0 : do k = 1, nfsd
153 : ! fracture_hist is already normalized
154 0 : omega(k) = afsd_init(k)*SUM(fracture_hist(1:k-1))
155 : end do
156 :
157 0 : loss = omega
158 :
159 0 : do k =1,nfsd
160 0 : gain(k) = SUM(omega*frac(:,k))
161 : end do
162 :
163 0 : d_afsd(:) = gain(:) - loss(:)
164 :
165 0 : 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 0 : WHERE (ABS(d_afsd).lt.puny) d_afsd = c0
171 :
172 0 : 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 96528 : subroutine icepack_step_wavefracture(wave_spec_type, &
183 : dt, nfreq, &
184 96528 : aice, vice, aicen, &
185 96528 : wave_spectrum, wavefreq, dwavefreq, &
186 96528 : 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
197 : aice, & ! ice area fraction
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)
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 193056 : 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, &
224 : nsubt ! number of subcycles
225 :
226 : real (kind=dbl_kind), dimension (nfsd, nfsd) :: &
227 193056 : frac
228 :
229 : real (kind=dbl_kind) :: &
230 : hbar , & ! mean ice thickness
231 : elapsed_t , & ! elapsed subcycling time
232 : subdt , & ! subcycling time step
233 : cons_error ! area conservation error
234 :
235 : real (kind=dbl_kind), dimension (nfsd) :: &
236 193056 : fracture_hist, & ! fracture histogram
237 96528 : afsd_init , & ! tracer array
238 193056 : afsd_tmp , & ! tracer array
239 96528 : 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 1254864 : d_afsd_wave (:) = c0
251 6370848 : d_afsdn_wave (:,:) = c0
252 1254864 : fracture_hist (:) = c0
253 :
254 : ! if all ice is not in first floe size category
255 96530 : if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then
256 :
257 2509728 : 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 96528 : if ((aice > p01).and.(local_sig_ht>0.1_dbl_kind)) then
261 :
262 72386 : hbar = vice / aice
263 :
264 : ! calculate fracture histogram
265 : call wave_frac(nfreq, wave_spec_type, &
266 : wavefreq, dwavefreq, &
267 72386 : hbar, wave_spectrum, fracture_hist)
268 :
269 72386 : if (icepack_warnings_aborted(subname)) return
270 :
271 : ! if fracture occurs
272 1013404 : if (MAXVAL(fracture_hist) > puny) then
273 : ! protect against small numerical errors
274 0 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
275 0 : if (icepack_warnings_aborted(subname)) return
276 :
277 0 : do n = 1, ncat
278 :
279 0 : 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 0 : .and. (afsd_init(1) < c1)) then
284 :
285 0 : afsd_tmp = afsd_init
286 :
287 : ! frac does not vary within subcycle
288 0 : frac(:,:) = c0
289 0 : do k = 2, nfsd
290 0 : frac(k,1:k-1) = fracture_hist(1:k-1)
291 : end do
292 0 : do k = 1, nfsd
293 0 : if (SUM(frac(k,:)) > c0) frac(k,:) = frac(k,:)/SUM(frac(k,:))
294 : end do
295 :
296 : ! adaptive sub-timestep
297 0 : elapsed_t = c0
298 0 : cons_error = c0
299 0 : nsubt = 0
300 0 : DO WHILE (elapsed_t < dt)
301 0 : nsubt = nsubt + 1
302 :
303 : ! if all floes in smallest category already, exit
304 0 : if (afsd_tmp(1).ge.c1-puny) EXIT
305 :
306 : ! calculate d_afsd using current afstd
307 0 : d_afsd_tmp = get_dafsd_wave(afsd_tmp, fracture_hist, frac)
308 :
309 : ! check in case wave fracture struggles to converge
310 0 : 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 0 : subdt = get_subdt_fsd(afsd_tmp, d_afsd_tmp)
318 0 : subdt = MIN(subdt, dt)
319 :
320 : ! update afsd
321 0 : afsd_tmp = afsd_tmp + subdt * d_afsd_tmp(:)
322 :
323 : ! check conservation and negatives
324 0 : if (MINVAL(afsd_tmp) < -puny) then
325 0 : write(warnstr,*) subname, 'wb, <0 loop'
326 0 : call icepack_warnings_add(warnstr)
327 : endif
328 0 : 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 0 : 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 0 : cons_error = SUM(afsd_tmp) - c1
346 :
347 : ! area loss: add to first category
348 0 : if (cons_error.lt.c0) then
349 0 : afsd_tmp(1) = afsd_tmp(1) - cons_error
350 : else
351 : ! area gain: take it from the largest possible category
352 0 : do k = nfsd, 1, -1
353 0 : if (afsd_tmp(k).gt.cons_error) then
354 0 : afsd_tmp(k) = afsd_tmp(k) - cons_error
355 0 : EXIT
356 : end if
357 : end do
358 : end if
359 :
360 : ! update trcrn
361 0 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp)
362 0 : call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
363 0 : if (icepack_warnings_aborted(subname)) return
364 :
365 : ! for diagnostics
366 0 : d_afsdn_wave(:,n) = afsd_tmp(:) - afsd_init(:)
367 0 : 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 72386 : subroutine wave_frac(nfreq, wave_spec_type, &
394 144772 : wavefreq, dwavefreq, &
395 72386 : 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)
408 : dwavefreq, & ! wave frequency bin widths (s^-1)
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 144772 : lambda, & ! wavelengths (m)
426 144772 : spec_coeff, &
427 144772 : 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)
434 : eta ! sea surface height field (m)
435 :
436 : real (kind=dbl_kind), dimension(nfsd) :: &
437 72386 : frachistogram, & ! histogram
438 72386 : prev_frac_local ! previous histogram
439 :
440 : character(len=*),parameter :: &
441 : subname='(wave_frac)'
442 :
443 :
444 72386 : 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 72386 : loop_max_iter = 1
449 : end if
450 :
451 : ! spatial domain
452 723932386 : do j = 1, nx
453 723932386 : X(j)= j*dx
454 : end do
455 :
456 : ! dispersion relation
457 1882036 : lambda (:) = gravit/(c2*pi*wavefreq (:)**2)
458 :
459 : ! spectral coefficients
460 1882036 : spec_coeff = sqrt(c2*spec_efreq*dwavefreq)
461 :
462 : ! initialize frac lengths
463 72386 : fraclengths(:) = c0
464 941018 : prev_frac_local(:) = c0
465 941018 : frachistogram(:) = c0
466 72386 : fracerror = bignum
467 :
468 : ! loop while fracerror greater than error tolerance
469 72386 : iter = 0
470 144772 : do while (iter < loop_max_iter .and. fracerror > errortol)
471 72386 : 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 72386 : 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 1882036 : rand_array(:) = p5
481 : endif
482 1882036 : phi = c2*pi*rand_array
483 :
484 723932386 : do j = 1, nx
485 : ! SSH field in space (sum over wavelengths, no attenuation)
486 18820360000 : summand = spec_coeff*COS(2*pi*X(j)/lambda+phi)
487 18820432386 : eta(j) = SUM(summand)
488 : end do
489 :
490 72386 : fraclengths(:) = c0
491 723932386 : if ((SUM(ABS(eta)) > puny).and.(hbar > puny)) then
492 72386 : call get_fraclengths(X, eta, fraclengths, hbar)
493 72386 : if (icepack_warnings_aborted(subname)) return
494 : end if
495 :
496 : ! convert from diameter to radii
497 723932386 : fraclengths(:) = fraclengths(:)/c2
498 :
499 723932386 : if (ALL(fraclengths.lt.floe_rad_l(1))) then
500 941018 : frac_local(:) = c0
501 : else
502 : ! bin into FS cats
503 : ! accumulate the frac histogram each iteration
504 0 : do j = 1, size(fraclengths)
505 0 : if (fraclengths(j).gt.floe_rad_l(1)) then
506 0 : do k = 1, nfsd-1
507 0 : if ((fraclengths(j) >= floe_rad_l(k)) .and. &
508 0 : (fraclengths(j) < floe_rad_l(k+1))) then
509 0 : frachistogram(k) = frachistogram(k) + 1
510 : end if
511 : end do
512 0 : if (fraclengths(j)>floe_rad_l(nfsd)) frachistogram(nfsd) = frachistogram(nfsd) + 1
513 : end if
514 : end do
515 :
516 0 : do k = 1, nfsd
517 0 : frac_local(k) = floe_rad_c(k)*frachistogram(k)
518 : end do
519 :
520 : ! normalize
521 0 : if (SUM(frac_local) /= c0) frac_local(:) = frac_local(:) / SUM(frac_local(:))
522 :
523 : end if
524 :
525 : ! wave fracture run to convergence
526 72386 : 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 72386 : 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 72386 : 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)
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
573 : j, k, & ! indices to iterate over domain
574 : first, last, & ! indices over which to search for extrema
575 : j_neg, & ! nearest extrema backwards
576 : j_pos, & ! nearest extrema forwards
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
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
585 : is_extremum, & ! or extremum
586 : is_triplet ! or triplet of extrema
587 :
588 : real (kind=dbl_kind) :: &
589 : denominator, & ! denominator in strain equation
590 : delta, & ! difference in x between current and prev extrema
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 72386 : spcing = nint(threshold/dx)
599 :
600 72386 : is_max = .false.
601 72386 : is_min = .false.
602 72386 : is_extremum = .false.
603 72386 : is_triplet = .false.
604 72386 : strain = c0
605 72386 : j_neg = 0
606 72386 : j_pos = 0
607 72386 : fraclengths(:) = c0
608 :
609 : ! search for local max and min within spacing
610 : ! on either side of each point
611 :
612 723932386 : do j = 1, nx
613 :
614 : ! indices within which to search for local max and min
615 723860000 : first = MAX(1,j-spcing)
616 723860000 : last = MIN(nx,j+spcing)
617 :
618 : ! location of max and min within spacing
619 723860000 : maxj = MAXLOC(eta(first:last))
620 723860000 : minj = MINLOC(eta(first:last))
621 :
622 : ! current j is the max or the min, save it
623 723860000 : if (maxj(1)+first-1 == j) is_max(j) = .true.
624 723860000 : if (minj(1)+first-1 == j) is_min(j) = .true.
625 :
626 : ! save whether max or min in one array
627 723932386 : 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 723787614 : do j = 2, nx-1
633 723787614 : if (is_extremum(j)) then
634 3329756 : if (j == 2) then
635 0 : if (is_extremum(1)) j_neg = 1
636 : else
637 716476628 : do k = j-1, 1, -1
638 716476628 : if (is_extremum(k)) then
639 3329756 : j_neg = k
640 3329756 : EXIT
641 : end if
642 : end do
643 : end if
644 :
645 706921676 : do k = j+1, nx
646 706921676 : if (is_extremum(k)) then
647 3329756 : j_pos = k
648 3329756 : EXIT
649 : end if
650 : end do
651 :
652 : ! find triplets of max and min
653 3329756 : if ((j_neg > 0).and.(j_pos > 0)) then
654 3329756 : if (is_max(j_neg).and.is_min(j).and.is_max(j_pos)) &
655 1664878 : is_triplet(j) = .true.
656 3329756 : if (is_min(j_neg).and.is_max(j).and.is_min(j_pos)) &
657 1664878 : is_triplet(j) = .true.
658 : end if
659 :
660 : ! calculate strain
661 3329756 : if (is_triplet(j)) then
662 :
663 : ! finite differences
664 3329756 : delta_pos = X(j_pos) - X(j )
665 3329756 : 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 3329756 : denominator = delta*delta_pos*(delta+delta_pos)
671 :
672 3329756 : if (denominator.ne.c0) &
673 : strain(j) = ABS(hbar*(eta(j_neg)* delta_pos &
674 : - eta(j )*(delta_pos+delta) &
675 : + eta(j_pos)* delta) &
676 3329756 : / denominator)
677 :
678 : end if ! is triplet
679 : end if ! is extremum
680 :
681 : end do ! loop over j
682 :
683 723932386 : n_above = COUNT(strain > straincrit)
684 72386 : fracdistances(:) = c0
685 :
686 : ! only do if we have some strains exceeding strain crit
687 72386 : if (n_above>0) then
688 :
689 0 : k = 0
690 0 : do j = 1, nx
691 0 : if (strain(j) > straincrit) then
692 0 : k = k + 1
693 0 : fracdistances(k) = X(j)
694 : end if
695 : end do
696 :
697 0 : do j = 1, n_above-1
698 0 : fraclengths(j) = fracdistances(j+1) - fracdistances(j)
699 : end do
700 :
701 :
702 : end if ! n_above
703 :
704 72386 : end subroutine get_fraclengths
705 :
706 : !=======================================================================
707 :
708 : end module icepack_wavefracspec
709 :
710 : !=======================================================================
711 :
712 :
|
