16 implicit none ;
private
18 #include <MOM_memory.h>
20 public wave_speed, wave_speeds, wave_speed_init, wave_speed_set_param
29 logical :: use_ebt_mode = .false.
33 real :: mono_n2_column_fraction = 0.
37 real :: mono_n2_depth = -1.
49 subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
50 mono_N2_column_fraction, mono_N2_depth, modal_structure)
54 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
57 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: cg1
59 logical,
optional,
intent(in) :: full_halos
61 logical,
optional,
intent(in) :: use_ebt_mode
63 real,
optional,
intent(in) :: mono_n2_column_fraction
66 real,
optional,
intent(in) :: mono_n2_depth
69 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
70 optional,
intent(out) :: modal_structure
73 real,
dimension(SZK_(G)+1) :: &
77 real,
dimension(SZK_(G)) :: &
80 real,
dimension(SZK_(G),SZI_(G)) :: &
82 real,
dimension(SZK_(G)) :: &
84 real det, ddet, detkm1, detkm2, ddetkm1, ddetkm2
85 real :: lam, dlam, lam0
88 real,
dimension(SZI_(G)) :: &
90 h_here, hxt_here, hxs_here, hxr_here
92 real :: i_hnew, drxh_sum
94 real,
parameter :: tol1 = 0.0001, tol2 = 0.001
95 real,
pointer,
dimension(:,:,:) :: t => null(), s => null()
97 real :: rescale, i_rescale
98 integer :: kf(szi_(g))
99 integer,
parameter :: max_itt = 10
100 real :: lam_it(max_itt), det_it(max_itt), ddet_it(max_itt)
104 integer :: i, j, k, k2, itt, is, ie, js, je, nz
105 real :: hw, gp, sum_hc, n2min
106 logical :: l_use_ebt_mode, calc_modal_structure
107 real :: l_mono_n2_column_fraction, l_mono_n2_depth
108 real :: mode_struct(szk_(g)), ms_min, ms_max, ms_sq
110 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
112 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_wave_speed: "// &
113 "Module must be initialized before it is used.")
114 if (
present(full_halos))
then ;
if (full_halos)
then
115 is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
118 l2_to_z2 = us%m_to_Z**2
120 l_use_ebt_mode = cs%use_ebt_mode
121 if (
present(use_ebt_mode)) l_use_ebt_mode = use_ebt_mode
122 l_mono_n2_column_fraction = cs%mono_N2_column_fraction
123 if (
present(mono_n2_column_fraction)) l_mono_n2_column_fraction = mono_n2_column_fraction
124 l_mono_n2_depth = us%m_to_Z*cs%mono_N2_depth
125 if (
present(mono_n2_depth)) l_mono_n2_depth = us%m_to_Z*mono_n2_depth
126 calc_modal_structure = l_use_ebt_mode
127 if (
present(modal_structure)) calc_modal_structure = .true.
128 if (calc_modal_structure)
then
129 do k=1,nz;
do j=js,je;
do i=is,ie
130 modal_structure(i,j,k) = 0.0
131 enddo ;
enddo ;
enddo
134 s => tv%S ; t => tv%T
135 g_rho0 = us%L_T_to_m_s**2 * gv%g_Earth / gv%Rho0
136 z_to_pa = gv%Z_to_H * gv%H_to_Pa
137 use_eos =
associated(tv%eqn_of_state)
139 rescale = 1024.0**4 ; i_rescale = 1.0/rescale
141 min_h_frac = tol1 / real(nz)
157 do i=is,ie ; htot(i) = 0.0 ;
enddo
158 do k=1,nz ;
do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ;
enddo ;
enddo
161 hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
162 hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
165 do k=1,nz ;
do i=is,ie
166 if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i)))
then
167 hf(kf(i),i) = h_here(i)
168 tf(kf(i),i) = hxt_here(i) / h_here(i)
169 sf(kf(i),i) = hxs_here(i) / h_here(i)
173 h_here(i) = h(i,j,k)*gv%H_to_Z
174 hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
175 hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
177 h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
178 hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
179 hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
182 do i=is,ie ;
if (h_here(i) > 0.0)
then
183 hf(kf(i),i) = h_here(i)
184 tf(kf(i),i) = hxt_here(i) / h_here(i)
185 sf(kf(i),i) = hxs_here(i) / h_here(i)
188 do k=1,nz ;
do i=is,ie
189 if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i)))
then
190 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
194 h_here(i) = h(i,j,k)*gv%H_to_Z
195 hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
197 h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
198 hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
201 do i=is,ie ;
if (h_here(i) > 0.0)
then
202 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
208 do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
212 pres(k) = pres(k-1) + z_to_pa*hf(k-1,i)
213 t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
214 s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
217 kf(i)-1, tv%eqn_of_state)
223 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
224 max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
225 drho_ds(k)*(sf(k,i)-sf(k-1,i)))
230 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
231 max(0.0,rf(k,i)-rf(k-1,i))
235 if (calc_modal_structure)
then
241 if (drxh_sum <= 0.0)
then
248 hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
250 if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
251 (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum)
then
253 i_hnew = 1.0 / (hc(kc) + hf(k,i))
254 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
255 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
256 hc(kc) = (hc(kc) + hf(k,i))
261 if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
262 (hc(k2) + hc(k2-1)) < tol2*drxh_sum)
then
264 i_hnew = 1.0 / (hc(kc) + hc(kc-1))
265 tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
266 sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
267 hc(kc-1) = (hc(kc) + hc(kc-1))
274 drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
275 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
280 gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
281 drho_ds(k)*(sc(k)-sc(k-1)))
286 hc(1) = hf(1,i) ; rc(1) = rf(1,i)
288 if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum)
then
290 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
291 hc(kc) = (hc(kc) + hf(k,i))
296 if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum)
then
298 rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
299 hc(kc-1) = (hc(kc) + hc(kc-1))
306 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
311 gprime(k) = g_rho0 * (rc(k)-rc(k-1))
320 if (l_use_ebt_mode)
then
323 n2min = l2_to_z2*gprime(2)/hc(1)
325 hw = 0.5*(hc(k-1)+hc(k))
327 if (l_mono_n2_column_fraction>0. .or. l_mono_n2_depth>=0.)
then
328 if ( ((g%bathyT(i,j)-sum_hc < l_mono_n2_column_fraction*g%bathyT(i,j)) .or. &
329 ((l_mono_n2_depth >= 0.) .and. (sum_hc > l_mono_n2_depth))) .and. &
330 (l2_to_z2*gp > n2min*hw) )
then
333 gp = us%Z_to_m**2 * (n2min*hw)
335 n2min = l2_to_z2 * gp/hw
338 igu(k) = 1.0/(gp*hc(k))
339 igl(k-1) = 1.0/(gp*hc(k-1))
340 speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
341 sum_hc = sum_hc + hc(k)
347 igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
348 speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
352 if (calc_modal_structure)
then
353 mode_struct(1:kc) = 1.
357 if (calc_modal_structure)
then
358 lam0 = 0.5 / speed2_tot ; lam = lam0
360 lam0 = 1.0 / speed2_tot ; lam = lam0
365 if (l_use_ebt_mode)
then
374 detkm1 = 1.0 ; ddetkm1 = 0.0
375 det = ( igl(1)-lam) ; ddet = -1.0
377 detkm2 = detkm1; ddetkm2 = ddetkm1
378 detkm1 = det; ddetkm1 = ddet
379 det = (igu(2)+igl(2)-lam)*detkm1 - (igu(2)*igl(1))*detkm2
380 ddet = (igu(2)+igl(2)-lam)*ddetkm1 - (igu(2)*igl(1))*ddetkm2 - detkm1
392 detkm1 = 1.0 ; ddetkm1 = 0.0
393 det = (igu(2)+igl(2)-lam) ; ddet = -1.0
400 detkm2 = detkm1; ddetkm2 = ddetkm1
401 detkm1 = det; ddetkm1 = ddet
403 det = (igu(k)+igl(k)-lam)*detkm1 - (igu(k)*igl(k-1))*detkm2
404 ddet = (igu(k)+igl(k)-lam)*ddetkm1 - (igu(k)*igl(k-1))*ddetkm2 - detkm1
407 if (abs(det) > rescale)
then
408 det = i_rescale*det ; detkm1 = i_rescale*detkm1
409 ddet = i_rescale*ddet ; ddetkm1 = i_rescale*ddetkm1
413 det_it(itt) = det ; ddet_it(itt) = ddet
415 if ((ddet >= 0.0) .or. (-det > -0.5*lam*ddet))
then
426 if (calc_modal_structure)
then
427 call tdma6(kc, -igu, igu+igl, -igl, lam, mode_struct)
428 ms_min = mode_struct(1)
429 ms_max = mode_struct(1)
430 ms_sq = mode_struct(1)**2
432 ms_min = min(ms_min, mode_struct(k))
433 ms_max = max(ms_max, mode_struct(k))
434 ms_sq = ms_sq + mode_struct(k)**2
436 if (ms_min<0. .and. ms_max>0.)
then
437 lam = 0.5 * ( lam - dlam )
439 mode_struct(1:kc) = abs(mode_struct(1:kc)) / sqrt( ms_sq )
441 mode_struct(1:kc) = mode_struct(1:kc) / sqrt( ms_sq )
445 if (abs(dlam) < tol2*lam)
exit
449 if (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)
451 if (
present(modal_structure))
then
452 if (mode_struct(1)/=0.)
then
453 mode_struct(1:kc) = mode_struct(1:kc) / mode_struct(1)
459 call remapping_core_h(cs%remapping_CS, kc, gv%Z_to_H*hc(:), mode_struct, &
460 nz, h(i,j,:), modal_structure(i,j,:), 1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
464 if (
present(modal_structure)) modal_structure(i,j,:) = 0.
469 if (
present(modal_structure)) modal_structure(i,j,:) = 0.
473 end subroutine wave_speed
477 subroutine tdma6(n, a, b, c, lam, y)
478 integer,
intent(in) :: n
479 real,
dimension(n),
intent(in) :: a
480 real,
dimension(n),
intent(in) :: b
481 real,
dimension(n),
intent(in) :: c
482 real,
intent(in) :: lam
483 real,
dimension(n),
intent(inout) :: y
486 real :: beta(n), yy(n), lambda
488 beta(1) = b(1) - lambda
489 if (beta(1)==0.)
then
491 lambda = (1. + 1.e-5) * lambda
492 beta(1) = b(1) - lambda
494 beta(1) = 1. / beta(1)
497 beta(k) = ( b(k) - lambda ) - a(k) * c(k-1) * beta(k-1)
498 if (beta(k)==0.)
then
500 lambda = (1. + 1.e-5) * lambda
501 beta(1) = 1. / ( b(1) - lambda )
503 beta(l) = 1. / ( ( b(l) - lambda ) - a(l) * c(l-1) * beta(l-1) )
504 yy(l) = y(l) - a(l) * yy(l-1) * beta(l-1)
507 beta(k) = 1. / beta(k)
509 yy(k) = y(k) - a(k) * yy(k-1) * beta(k-1)
511 y(n) = yy(n) * beta(n)
513 y(k) = ( yy(k) - c(k) * y(k+1) ) * beta(k)
518 subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
522 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
524 integer,
intent(in) :: nmodes
525 real,
dimension(G%isd:G%ied,G%jsd:G%jed,nmodes),
intent(out) :: cn
527 logical,
optional,
intent(in) :: full_halos
530 real,
dimension(SZK_(G)+1) :: &
532 pres, t_int, s_int, &
534 real,
dimension(SZK_(G)) :: &
537 real,
dimension(SZK_(G)-1) :: &
538 a_diag, b_diag, c_diag
541 real,
dimension(SZK_(G),SZI_(G)) :: &
543 real,
dimension(SZK_(G)) :: &
545 real,
parameter :: c1_thresh = 0.01
556 real :: ddet_l,ddet_r
557 real :: det_sub,ddet_sub
560 real,
dimension(nmodes) :: &
563 integer :: nrootsfound
566 real,
dimension(SZI_(G)) :: &
568 h_here, hxt_here, hxs_here, hxr_here
571 real,
parameter :: reduct_factor = 0.5
573 real :: i_hnew, drxh_sum
574 real,
parameter :: tol1 = 0.0001, tol2 = 0.001
575 real,
pointer,
dimension(:,:,:) :: t => null(), s => null()
577 integer :: kf(szi_(g))
578 integer,
parameter :: max_itt = 10
580 real,
dimension(SZK_(G)+1) :: z_int, n2
582 integer,
parameter :: sub_it_max = 4
585 logical :: sub_rootfound
587 integer :: sub, sub_it
588 integer :: i, j, k, k2, itt, is, ie, js, je, nz, row, iint, m, ig, jg
590 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
592 if (
present(cs))
then
593 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_wave_speed: "// &
594 "Module must be initialized before it is used.")
597 if (
present(full_halos))
then ;
if (full_halos)
then
598 is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
601 s => tv%S ; t => tv%T
602 g_rho0 = us%L_T_to_m_s**2 * gv%g_Earth / gv%Rho0
603 use_eos =
associated(tv%eqn_of_state)
604 z_to_pa = gv%Z_to_H * gv%H_to_Pa
606 min_h_frac = tol1 / real(nz)
613 do i=is,ie ; htot(i) = 0.0 ;
enddo
614 do k=1,nz ;
do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ;
enddo ;
enddo
617 hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
618 hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
621 do k=1,nz ;
do i=is,ie
622 if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i)))
then
623 hf(kf(i),i) = h_here(i)
624 tf(kf(i),i) = hxt_here(i) / h_here(i)
625 sf(kf(i),i) = hxs_here(i) / h_here(i)
629 h_here(i) = h(i,j,k)*gv%H_to_Z
630 hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
631 hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
633 h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
634 hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
635 hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
638 do i=is,ie ;
if (h_here(i) > 0.0)
then
639 hf(kf(i),i) = h_here(i)
640 tf(kf(i),i) = hxt_here(i) / h_here(i)
641 sf(kf(i),i) = hxs_here(i) / h_here(i)
644 do k=1,nz ;
do i=is,ie
645 if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i)))
then
646 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
650 h_here(i) = h(i,j,k)*gv%H_to_Z
651 hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
653 h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
654 hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
657 do i=is,ie ;
if (h_here(i) > 0.0)
then
658 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
665 if (g%mask2dT(i,j) > 0.5)
then
669 pres(k) = pres(k-1) + z_to_pa*hf(k-1,i)
670 t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
671 s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
674 kf(i)-1, tv%eqn_of_state)
680 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
681 max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
682 drho_ds(k)*(sf(k,i)-sf(k-1,i)))
687 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
688 max(0.0,rf(k,i)-rf(k-1,i))
693 if (drxh_sum <= 0.0)
then
700 hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
702 if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
703 (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum)
then
705 i_hnew = 1.0 / (hc(kc) + hf(k,i))
706 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
707 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
708 hc(kc) = (hc(kc) + hf(k,i))
713 if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
714 (hc(k2) + hc(k2-1)) < tol2*drxh_sum)
then
716 i_hnew = 1.0 / (hc(kc) + hc(kc-1))
717 tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
718 sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
719 hc(kc-1) = (hc(kc) + hc(kc-1))
726 drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
727 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
732 gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
733 drho_ds(k)*(sc(k)-sc(k-1)))
738 hc(1) = hf(1,i) ; rc(1) = rf(1,i)
740 if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum)
then
742 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
743 hc(kc) = (hc(kc) + hf(k,i))
748 if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum)
then
750 rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
751 hc(kc-1) = (hc(kc) + hc(kc-1))
758 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
763 gprime(k) = g_rho0 * (rc(k)-rc(k-1))
768 ig = i + g%idg_offset ; jg = j + g%jdg_offset
779 igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
780 z_int(k) = z_int(k-1) + hc(k-1)
781 n2(k) = us%m_to_Z**2*gprime(k)/(0.5*(hc(k)+hc(k-1)))
782 speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
785 n2(1) = n2(2) ; n2(kc+1) = n2(kc)
787 z_int(kc+1) = z_int(kc)+hc(kc)
789 if (abs(z_int(kc+1)-htot(i)) > 1.e-12*htot(i))
then
790 call mom_error(fatal,
"wave_structure: mismatch in total depths")
797 a_diag(row) = (-igu(k))
798 b_diag(row) = (igu(k)+igl(k))
799 c_diag(row) = (-igl(k))
804 b_diag(row) = (igu(k)+igl(k))
805 c_diag(row) = (-igl(k))
808 a_diag(row) = (-igu(k))
809 b_diag(row) = (igu(k)+igl(k))
815 lam_1 = 1.0 / speed2_tot
820 call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
821 nrows,lam_1,det,ddet)
824 if ((ddet >= 0.0) .or. (-det > -0.5*lam_1*ddet))
then
832 if (abs(dlam) < tol2*lam_1)
then
834 if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1)
842 if (nmodes>1 .and. kc>=nmodes+1 .and. cn(i,j,1)>c1_thresh)
then
845 lammin = lam_1*(1.0 + tol2)
847 speed2_min = (reduct_factor*cn(i,j,1)/real(nmodes))**2
848 lammax = 1.0 / speed2_min
852 numint = nint((lammax - lammin)/laminc)
858 call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
859 nrows,lammin,det_l,ddet_l)
862 xr = lammin + laminc * iint
864 call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
865 nrows,xr,det_r,ddet_r)
866 if (det_l*det_r < 0.0)
then
867 if (det_l*ddet_l < 0.0)
then
868 nrootsfound = nrootsfound + 1
869 xbl(nrootsfound) = xl
870 xbr(nrootsfound) = xr
879 sub_rootfound = .false.
880 do sub_it=1,sub_it_max
884 xl_sub = xl + laminc/(nsub)*sub
885 call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
886 nrows,xl_sub,det_sub,ddet_sub)
887 if (det_sub*det_r < 0.0)
then
888 if (det_sub*ddet_sub < 0.0)
then
889 sub_rootfound = .true.
890 nrootsfound = nrootsfound + 1
891 xbl(nrootsfound) = xl_sub
892 xbr(nrootsfound) = xr
897 if (sub_rootfound)
exit
900 if (sub_it == sub_it_max)
then
901 call mom_error(warning,
"wave_speed: root not found "// &
902 " after sub_it_max subdivisions of original"// &
909 if (nrootsfound >= nmodes-1)
then
912 elseif (iint == numint)
then
915 cn(i,j,nrootsfound+2:nmodes) = 0.0
928 call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
929 nrows,lam_n,det,ddet)
933 if (abs(dlam) < tol2*lam_1)
then
935 if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n)
941 cn(i,j,2:nmodes) = 0.0
953 end subroutine wave_speeds
956 subroutine tridiag_det(a,b,c,nrows,lam,det_out,ddet_out)
957 real,
dimension(:),
intent(in) :: a
958 real,
dimension(:),
intent(in) :: b
959 real,
dimension(:),
intent(in) :: c
960 integer,
intent(in) :: nrows
961 real,
intent(in) :: lam
962 real,
intent(out):: det_out
963 real,
intent(out):: ddet_out
965 real,
dimension(nrows) :: det
966 real,
dimension(nrows) :: ddet
967 real,
parameter:: rescale = 1024.0**4
971 if (
size(b) /= nrows)
call mom_error(warning,
"Diagonal b must be same length as nrows.")
972 if (
size(a) /= nrows)
call mom_error(warning,
"Diagonal a must be same length as nrows.")
973 if (
size(c) /= nrows)
call mom_error(warning,
"Diagonal c must be same length as nrows.")
975 i_rescale = 1.0/rescale
977 det(1) = 1.0 ; ddet(1) = 0.0
978 det(2) = b(2)-lam ; ddet(2) = -1.0
980 det(n) = (b(n)-lam)*det(n-1) - (a(n)*c(n-1))*det(n-2)
981 ddet(n) = (b(n)-lam)*ddet(n-1) - (a(n)*c(n-1))*ddet(n-2) - det(n-1)
983 if (abs(det(n)) > rescale)
then
984 det(n) = i_rescale*det(n) ; det(n-1) = i_rescale*det(n-1)
985 ddet(n) = i_rescale*ddet(n) ; ddet(n-1) = i_rescale*ddet(n-1)
989 ddet_out = ddet(nrows)
991 end subroutine tridiag_det
994 subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth)
996 logical,
optional,
intent(in) :: use_ebt_mode
998 real,
optional,
intent(in) :: mono_n2_column_fraction
1001 real,
optional,
intent(in) :: mono_n2_depth
1005 #include "version_variable.h"
1006 character(len=40) :: mdl =
"MOM_wave_speed"
1008 if (
associated(cs))
then
1009 call mom_error(warning,
"wave_speed_init called with an "// &
1010 "associated control structure.")
1012 else ;
allocate(cs) ;
endif
1017 call wave_speed_set_param(cs, use_ebt_mode=use_ebt_mode, mono_n2_column_fraction=mono_n2_column_fraction)
1019 call initialize_remapping(cs%remapping_CS,
'PLM', boundary_extrapolation=.false.)
1021 end subroutine wave_speed_init
1024 subroutine wave_speed_set_param(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth)
1026 logical,
optional,
intent(in) :: use_ebt_mode
1028 real,
optional,
intent(in) :: mono_n2_column_fraction
1031 real,
optional,
intent(in) :: mono_n2_depth
1035 if (.not.
associated(cs))
call mom_error(fatal, &
1036 "wave_speed_set_param called with an associated control structure.")
1038 if (
present(use_ebt_mode)) cs%use_ebt_mode = use_ebt_mode
1039 if (
present(mono_n2_column_fraction)) cs%mono_N2_column_fraction = mono_n2_column_fraction
1040 if (
present(mono_n2_depth)) cs%mono_N2_depth = mono_n2_depth
1042 end subroutine wave_speed_set_param