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)
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 884760 : 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 0 : 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(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 61488 : subroutine icepack_step_wavefracture(wave_spec_type, &
186 : dt, ncat, nfsd, &
187 : nfreq, &
188 96528 : aice, vice, aicen, &
189 96528 : floe_rad_l, floe_rad_c, &
190 193056 : wave_spectrum, wavefreq, dwavefreq, &
191 96528 : 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
199 : ncat, & ! number of thickness categories
200 : nfsd ! number of floe size categories
201 :
202 : real (kind=dbl_kind), intent(in) :: &
203 : dt, & ! time step
204 : aice, & ! ice area fraction
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)
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)
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 2435616 : 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, t, &
235 : nsubt ! number of subcycles
236 :
237 : real (kind=dbl_kind), dimension(nfsd,ncat) :: &
238 : afsdn ! floe size and thickness distribution
239 :
240 : real (kind=dbl_kind), dimension (nfsd, nfsd) :: &
241 5624256 : frac
242 :
243 : real (kind=dbl_kind) :: &
244 35040 : hbar , & ! mean ice thickness
245 35040 : elapsed_t , & ! elapsed subcycling time
246 35040 : subdt , & ! subcycling time step
247 35040 : cons_error ! area conservation error
248 :
249 : real (kind=dbl_kind), dimension (nfsd) :: &
250 578496 : fracture_hist, & ! fracture histogram
251 648576 : afsd_init , & ! tracer array
252 578496 : afsd_tmp , & ! tracer array
253 587088 : d_afsd_tmp ! change
254 :
255 : character(len=*),parameter :: &
256 : subname='(icepack_step_wavefracture)'
257 :
258 : !------------------------------------
259 :
260 : ! initialize
261 1254864 : d_afsd_wave (:) = c0
262 6370848 : d_afsdn_wave (:,:) = c0
263 1254864 : fracture_hist (:) = c0
264 :
265 : ! if all ice is not in first floe size category
266 96535 : if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then
267 :
268 :
269 : ! do not try to fracture for minimal ice concentration or zero wave spectrum
270 2606256 : if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then
271 :
272 67558 : hbar = vice / aice
273 :
274 : ! calculate fracture histogram
275 : call wave_frac(nfsd, nfreq, wave_spec_type, &
276 0 : floe_rad_l, floe_rad_c, &
277 0 : wavefreq, dwavefreq, &
278 67558 : hbar, wave_spectrum, fracture_hist)
279 :
280 67558 : if (icepack_warnings_aborted(subname)) return
281 :
282 : ! if fracture occurs
283 945812 : if (MAXVAL(fracture_hist) > puny) then
284 : ! protect against small numerical errors
285 0 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
286 0 : if (icepack_warnings_aborted(subname)) return
287 :
288 0 : do n = 1, ncat
289 :
290 0 : afsd_init(:) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
291 :
292 : ! if there is ice, and a FSD, and not all ice is the smallest floe size
293 0 : if ((aicen(n) > puny) .and. (SUM(afsd_init(:)) > puny) &
294 0 : .and. (afsd_init(1) < c1)) then
295 :
296 0 : afsd_tmp = afsd_init
297 :
298 : ! frac does not vary within subcycle
299 0 : frac(:,:) = c0
300 0 : do k = 2, nfsd
301 0 : frac(k,1:k-1) = fracture_hist(1:k-1)
302 : end do
303 0 : do k = 1, nfsd
304 0 : if (SUM(frac(k,:)) > c0) frac(k,:) = frac(k,:)/SUM(frac(k,:))
305 : end do
306 :
307 : ! adaptive sub-timestep
308 0 : elapsed_t = c0
309 0 : cons_error = c0
310 0 : nsubt = 0
311 0 : DO WHILE (elapsed_t < dt)
312 0 : nsubt = nsubt + 1
313 :
314 : ! if all floes in smallest category already, exit
315 0 : if (afsd_tmp(1).ge.c1-puny) EXIT
316 :
317 : ! calculate d_afsd using current afstd
318 0 : d_afsd_tmp = get_dafsd_wave(nfsd, afsd_tmp, fracture_hist, frac)
319 :
320 : ! check in case wave fracture struggles to converge
321 0 : if (nsubt>100) then
322 0 : write(warnstr,*) subname, &
323 0 : 'warning: step_wavefracture struggling to converge'
324 0 : call icepack_warnings_add(warnstr)
325 : endif
326 :
327 : ! required timestep
328 0 : subdt = get_subdt_fsd(nfsd, afsd_tmp, d_afsd_tmp)
329 0 : subdt = MIN(subdt, dt)
330 :
331 : ! update afsd
332 0 : afsd_tmp = afsd_tmp + subdt * d_afsd_tmp(:)
333 :
334 : ! check conservation and negatives
335 0 : if (MINVAL(afsd_tmp) < -puny) then
336 0 : write(warnstr,*) subname, 'wb, <0 loop'
337 0 : call icepack_warnings_add(warnstr)
338 : endif
339 0 : if (MAXVAL(afsd_tmp) > c1+puny) then
340 0 : write(warnstr,*) subname, 'wb, >1 loop'
341 0 : call icepack_warnings_add(warnstr)
342 : endif
343 :
344 : ! update time
345 0 : elapsed_t = elapsed_t + subdt
346 :
347 : END DO ! elapsed_t < dt
348 :
349 : ! In some cases---particularly for strong fracturing---the equation
350 : ! for wave fracture does not quite conserve area.
351 : ! With the dummy wave forcing, this happens < 2% of the time (in
352 : ! 1997) and is always less than 10^-7.
353 : ! Simply renormalizing may cause the first floe size
354 : ! category to reduce, which is not physically allowed
355 : ! to happen. So we adjust here
356 0 : cons_error = SUM(afsd_tmp) - c1
357 :
358 : ! area loss: add to first category
359 0 : if (cons_error.lt.c0) then
360 0 : afsd_tmp(1) = afsd_tmp(1) - cons_error
361 : else
362 : ! area gain: take it from the largest possible category
363 0 : do k = nfsd, 1, -1
364 0 : if (afsd_tmp(k).gt.cons_error) then
365 0 : afsd_tmp(k) = afsd_tmp(k) - cons_error
366 0 : EXIT
367 : end if
368 : end do
369 : end if
370 :
371 : ! update trcrn
372 0 : trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp)
373 0 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
374 0 : if (icepack_warnings_aborted(subname)) return
375 :
376 : ! for diagnostics
377 0 : d_afsdn_wave(:,n) = afsd_tmp(:) - afsd_init(:)
378 0 : d_afsd_wave (:) = d_afsd_wave(:) + aicen(n)*d_afsdn_wave(:,n)
379 : endif ! aicen > puny
380 : enddo ! n
381 : endif ! fracture hist > 0
382 :
383 : endif ! aice > p01
384 : endif ! all small floes
385 :
386 : end subroutine icepack_step_wavefracture
387 :
388 : !=======================================================================
389 : !
390 : ! Calculates functions to describe the change in the FSD when waves
391 : ! fracture ice, given a wave spectrum (1D frequency, nfreq (default 25)
392 : ! frequency bins)
393 : !
394 : ! We calculate extrema and if these are successive maximum,
395 : ! minimum, maximum or vice versa, and have strain greater than a
396 : ! critical strain, break ice and create new floes with lengths equal
397 : ! to these distances. Based on MatLab code written by Chris Horvat,
398 : ! from Horvat & Tziperman (2015).
399 : !
400 : ! Note that a realization of sea surface height requires a random phase.
401 : !
402 : ! authors: 2018 Lettie Roach, NIWA/VUW
403 :
404 42991 : subroutine wave_frac(nfsd, nfreq, wave_spec_type, &
405 135116 : floe_rad_l, floe_rad_c, &
406 135116 : wavefreq, dwavefreq, &
407 135116 : hbar, spec_efreq, frac_local)
408 :
409 : integer (kind=int_kind), intent(in) :: &
410 : nfsd, & ! number of floe size categories
411 : nfreq ! number of wave frequency categories
412 :
413 : character (len=char_len), intent(in) :: &
414 : wave_spec_type ! type of wave spectrum forcing
415 :
416 : real (kind=dbl_kind), intent(in) :: &
417 : hbar ! mean ice thickness (m)
418 :
419 : real(kind=dbl_kind), dimension(:), intent(in) :: &
420 : floe_rad_l, & ! fsd size lower bound in m (radius)
421 : floe_rad_c ! fsd size bin centre in m (radius)
422 :
423 : real (kind=dbl_kind), dimension (:), intent(in) :: &
424 : wavefreq, & ! wave frequencies (s^-1)
425 : dwavefreq, & ! wave frequency bin widths (s^-1)
426 : spec_efreq ! wave spectrum (m^2 s)
427 :
428 : real (kind=dbl_kind), dimension (nfsd), intent(out) :: &
429 : frac_local ! fracturing histogram
430 :
431 : ! local variables
432 :
433 : integer (kind=int_kind) :: i, j, k, iter, loop_max_iter
434 :
435 : real (kind=dbl_kind) :: &
436 24567 : fracerror ! difference between successive histograms
437 :
438 : real (kind=dbl_kind), parameter :: &
439 : errortol = 6.5e-4 ! tolerance in error between successive histograms
440 :
441 : real (kind=dbl_kind), dimension(nfreq) :: &
442 724724 : lambda, & ! wavelengths (m)
443 724724 : spec_coeff, &
444 2174172 : phi, rand_array, summand
445 :
446 : real (kind=dbl_kind), dimension(nx) :: &
447 : fraclengths
448 :
449 : real (kind=dbl_kind), dimension(nx) :: &
450 : X, & ! spatial domain (m)
451 : eta ! sea surface height field (m)
452 :
453 : real (kind=dbl_kind), dimension(nfsd) :: &
454 454487 : frachistogram, & ! histogram
455 411496 : prev_frac_local ! previous histogram
456 :
457 : character(len=*),parameter :: &
458 : subname='(wave_frac)'
459 :
460 :
461 67558 : if (trim(wave_spec_type).eq.'random') then
462 : ! run wave fracture to convergence
463 0 : loop_max_iter = max_no_iter
464 : else
465 67558 : loop_max_iter = 1
466 : end if
467 :
468 : ! spatial domain
469 675647558 : do j = 1, nx
470 675647558 : X(j)= j*dx
471 : end do
472 :
473 : ! dispersion relation
474 1756508 : lambda (:) = gravit/(c2*pi*wavefreq (:)**2)
475 :
476 : ! spectral coefficients
477 1756508 : spec_coeff = sqrt(c2*spec_efreq*dwavefreq)
478 :
479 : ! initialize frac lengths
480 67558 : fraclengths(:) = c0
481 878254 : prev_frac_local(:) = c0
482 878254 : frachistogram(:) = c0
483 67558 : fracerror = bignum
484 :
485 : ! loop while fracerror greater than error tolerance
486 67558 : iter = 0
487 135116 : do while (iter < loop_max_iter .and. fracerror > errortol)
488 67558 : iter = iter + 1
489 :
490 : ! Phase for each Fourier component may be constant or
491 : ! a random phase that varies in each i loop
492 : ! See documentation for discussion
493 67558 : if (trim(wave_spec_type).eq.'random') then
494 0 : call RANDOM_NUMBER(rand_array)
495 0 : if (icepack_warnings_aborted(subname)) return
496 : else
497 1756508 : rand_array(:) = p5
498 : endif
499 1756508 : phi = c2*pi*rand_array
500 :
501 675647558 : do j = 1, nx
502 : ! SSH field in space (sum over wavelengths, no attenuation)
503 17565080000 : summand = spec_coeff*COS(2*pi*X(j)/lambda+phi)
504 17565147558 : eta(j) = SUM(summand)
505 : end do
506 :
507 67558 : fraclengths(:) = c0
508 675647558 : if ((SUM(ABS(eta)) > puny).and.(hbar > puny)) then
509 67558 : call get_fraclengths(X, eta, fraclengths, hbar)
510 67558 : if (icepack_warnings_aborted(subname)) return
511 : end if
512 :
513 : ! convert from diameter to radii
514 675647558 : fraclengths(:) = fraclengths(:)/c2
515 :
516 675647558 : if (ALL(fraclengths.lt.floe_rad_l(1))) then
517 878254 : frac_local(:) = c0
518 : else
519 : ! bin into FS cats
520 : ! accumulate the frac histogram each iteration
521 0 : do j = 1, size(fraclengths)
522 0 : if (fraclengths(j).gt.floe_rad_l(1)) then
523 0 : do k = 1, nfsd-1
524 0 : if ((fraclengths(j) >= floe_rad_l(k)) .and. &
525 0 : (fraclengths(j) < floe_rad_l(k+1))) then
526 0 : frachistogram(k) = frachistogram(k) + 1
527 : end if
528 : end do
529 0 : if (fraclengths(j)>floe_rad_l(nfsd)) frachistogram(nfsd) = frachistogram(nfsd) + 1
530 : end if
531 : end do
532 :
533 0 : do k = 1, nfsd
534 0 : frac_local(k) = floe_rad_c(k)*frachistogram(k)
535 : end do
536 :
537 : ! normalize
538 0 : if (SUM(frac_local) /= c0) frac_local(:) = frac_local(:) / SUM(frac_local(:))
539 :
540 : end if
541 :
542 : ! wave fracture run to convergence
543 135116 : if (trim(wave_spec_type).eq.'random') then
544 :
545 : ! check avg frac local against previous iteration
546 0 : fracerror = SUM(ABS(frac_local - prev_frac_local))/nfsd
547 :
548 : ! save histogram for next iteration
549 0 : prev_frac_local = frac_local
550 :
551 : end if
552 :
553 : END DO
554 :
555 67558 : if (iter >= max_no_iter) then
556 0 : write(warnstr,*) subname,'warning: wave_frac struggling to converge'
557 0 : call icepack_warnings_add(warnstr)
558 : endif
559 :
560 : end subroutine wave_frac
561 :
562 : !===========================================================================
563 : !
564 : ! Given the (attenuated) sea surface height, find the strain across triplets
565 : ! of max, min, max or min, max, min (local extrema within 10m).
566 : ! If this strain is greater than the critical strain, ice can fracture
567 : ! and new floes are formed with sizes equal to the distances between
568 : ! extrema. Based on MatLab code written by Chris Horvat,
569 : ! from Horvat & Tziperman (2015).
570 : !
571 : ! authors: 2016 Lettie Roach, NIWA/VUW
572 : !
573 67558 : subroutine get_fraclengths(X, eta, fraclengths, hbar)
574 :
575 : real (kind=dbl_kind), intent(in) :: &
576 : hbar ! mean thickness (m)
577 :
578 : real (kind=dbl_kind), intent(in), dimension (nx) :: &
579 : X, & ! spatial domain (m)
580 : eta ! sea surface height field (m)
581 :
582 : real (kind=dbl_kind), intent(inout), dimension (nx) :: &
583 : fraclengths ! The distances between fracture points
584 : ! Size cannot be greater than nx.
585 : ! In practice, will be much less
586 :
587 : ! local variables
588 : integer (kind=int_kind) :: &
589 : spcing, & ! distance over which to search for extrema on each side of point
590 : j, k, & ! indices to iterate over domain
591 : first, last, & ! indices over which to search for extrema
592 : j_neg, & ! nearest extrema backwards
593 : j_pos, & ! nearest extrema forwards
594 : n_above ! number of points where strain is above critical
595 :
596 : real (kind=dbl_kind), dimension(nx) :: &
597 : fracdistances, & ! distances in space where fracture has occurred
598 : strain ! the strain between triplets of extrema
599 :
600 : logical (kind=log_kind), dimension(nx) :: &
601 : is_max, is_min,& ! arrays to hold whether each point is a local max or min
602 : is_extremum, & ! or extremum
603 : is_triplet ! or triplet of extrema
604 :
605 : real (kind=dbl_kind) :: &
606 24567 : denominator, & ! denominator in strain equation
607 24567 : delta, & ! difference in x between current and prev extrema
608 24567 : delta_pos ! difference in x between next and current extrema
609 :
610 : integer (kind=int_kind), dimension(1) :: &
611 : maxj, minj ! indices of local max and min
612 :
613 : ! ------- equivalent of peakfinder2
614 : ! given eta and spcing, compute extremelocs in ascending order
615 67558 : spcing = nint(threshold/dx)
616 :
617 67558 : is_max = .false.
618 67558 : is_min = .false.
619 67558 : is_extremum = .false.
620 67558 : is_triplet = .false.
621 67558 : strain = c0
622 67558 : j_neg = 0
623 67558 : j_pos = 0
624 67558 : fraclengths(:) = c0
625 :
626 : ! search for local max and min within spacing
627 : ! on either side of each point
628 :
629 675647558 : do j = 1, nx
630 :
631 : ! indices within which to search for local max and min
632 675580000 : first = MAX(1,j-spcing)
633 675580000 : last = MIN(nx,j+spcing)
634 :
635 : ! location of max and min within spacing
636 675580000 : maxj = MAXLOC(eta(first:last))
637 675580000 : minj = MINLOC(eta(first:last))
638 :
639 : ! current j is the max or the min, save it
640 675580000 : if (maxj(1)+first-1 == j) is_max(j) = .true.
641 675580000 : if (minj(1)+first-1 == j) is_min(j) = .true.
642 :
643 : ! save whether max or min in one array
644 675647558 : if (is_min(j).or.is_max(j)) is_extremum(j) = .true.
645 : end do
646 :
647 : ! loop over points
648 : ! nothing can happen at the first or last
649 675512442 : do j = 2, nx-1
650 675512442 : if (is_extremum(j)) then
651 3107668 : if (j == 2) then
652 0 : if (is_extremum(1)) j_neg = 1
653 : else
654 668689084 : do k = j-1, 1, -1
655 668689084 : if (is_extremum(k)) then
656 3107668 : j_neg = k
657 3107668 : EXIT
658 : end if
659 : end do
660 : end if
661 :
662 659771428 : do k = j+1, nx
663 659771428 : if (is_extremum(k)) then
664 3107668 : j_pos = k
665 3107668 : EXIT
666 : end if
667 : end do
668 :
669 : ! find triplets of max and min
670 3107668 : if ((j_neg > 0).and.(j_pos > 0)) then
671 3107668 : if (is_max(j_neg).and.is_min(j).and.is_max(j_pos)) &
672 1553834 : is_triplet(j) = .true.
673 3107668 : if (is_min(j_neg).and.is_max(j).and.is_min(j_pos)) &
674 1553834 : is_triplet(j) = .true.
675 : end if
676 :
677 : ! calculate strain
678 3107668 : if (is_triplet(j)) then
679 :
680 : ! finite differences
681 3107668 : delta_pos = X(j_pos) - X(j )
682 3107668 : delta = X(j ) - X(j_neg)
683 :
684 : ! This equation differs from HT2015 by a factor 2 in numerator
685 : ! and eta(j_pos). This is the correct form of the equation.
686 :
687 3107668 : denominator = delta*delta_pos*(delta+delta_pos)
688 :
689 3107668 : if (denominator.ne.c0) &
690 1130082 : strain(j) = ABS(hbar*(eta(j_neg)* delta_pos &
691 1130082 : - eta(j )*(delta_pos+delta) &
692 1130082 : + eta(j_pos)* delta) &
693 3107668 : / denominator)
694 :
695 : end if ! is triplet
696 : end if ! is extremum
697 :
698 : end do ! loop over j
699 :
700 675647558 : n_above = COUNT(strain > straincrit)
701 67558 : fracdistances(:) = c0
702 :
703 : ! only do if we have some strains exceeding strain crit
704 67558 : if (n_above>0) then
705 :
706 0 : k = 0
707 0 : do j = 1, nx
708 0 : if (strain(j) > straincrit) then
709 0 : k = k + 1
710 0 : fracdistances(k) = X(j)
711 : end if
712 : end do
713 :
714 0 : do j = 1, n_above-1
715 0 : fraclengths(j) = fracdistances(j+1) - fracdistances(j)
716 : end do
717 :
718 :
719 : end if ! n_above
720 :
721 67558 : end subroutine get_fraclengths
722 :
723 : !=======================================================================
724 :
725 : end module icepack_wavefracspec
726 :
727 : !=======================================================================
728 :
729 :
|
