24 implicit none ;
private
26 #include <MOM_memory.h>
39 real,
allocatable,
dimension(:,:,:) :: w_strct
41 real,
allocatable,
dimension(:,:,:) :: u_strct
43 real,
allocatable,
dimension(:,:,:) :: w_profile
47 real,
allocatable,
dimension(:,:,:) :: uavg_profile
50 real,
allocatable,
dimension(:,:,:) :: z_depths
52 real,
allocatable,
dimension(:,:,:) :: n2
54 integer,
allocatable,
dimension(:,:):: num_intfaces
56 real :: int_tide_source_x
58 real :: int_tide_source_y
91 subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halos)
95 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
98 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: cn
100 integer,
intent(in) :: modenum
101 real,
intent(in) :: freq
104 real,
dimension(SZI_(G),SZJ_(G)), &
105 optional,
intent(in) :: en
106 logical,
optional,
intent(in) :: full_halos
109 real,
dimension(SZK_(G)+1) :: &
116 real,
dimension(SZK_(G)) :: &
119 real,
dimension(SZK_(G),SZI_(G)) :: &
124 real,
dimension(SZK_(G)) :: &
130 real,
dimension(SZI_(G),SZJ_(G)) :: &
135 real,
dimension(SZI_(G)) :: &
144 real,
parameter :: tol1 = 0.0001, tol2 = 0.001
145 real,
pointer,
dimension(:,:,:) :: t => null(), s => null()
148 integer :: kf(szi_(g))
149 integer,
parameter :: max_itt = 1
151 real,
parameter :: a_int = 0.5
157 real,
dimension(SZK_(G)+1) :: w_strct, u_strct, w_profile, uavg_profile, z_int, n2
160 real,
dimension(SZK_(G)+1) :: w_strct2, u_strct2
162 real,
dimension(SZK_(G)) :: dz
172 real,
dimension(SZK_(G)-1) :: lam_z
174 real,
dimension(SZK_(G)-1) :: a_diag, b_diag, c_diag
177 real,
dimension(SZK_(G)-1) :: e_guess
178 real,
dimension(SZK_(G)-1) :: e_itt
181 integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
183 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
187 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_wave_structure: "// &
188 "Module must be initialized before it is used.")
191 if (
present(full_halos))
then ;
if (full_halos)
then
192 is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
197 s => tv%S ; t => tv%T
198 g_rho0 = us%L_T_to_m_s**2 * gv%g_Earth / gv%Rho0
199 cg_subro = 1e-100*us%m_s_to_L_T
200 use_eos =
associated(tv%eqn_of_state)
202 h_to_pres = gv%Z_to_H*gv%H_to_Pa
205 min_h_frac = tol1 / real(nz)
211 do i=is,ie ; htot(i,j) = 0.0 ;
enddo
212 do k=1,nz ;
do i=is,ie ; htot(i,j) = htot(i,j) + h(i,j,k)*gv%H_to_Z ;
enddo ;
enddo
215 hmin(i) = htot(i,j)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
216 hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
219 do k=1,nz ;
do i=is,ie
220 if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i)))
then
221 hf(kf(i),i) = h_here(i)
222 tf(kf(i),i) = hxt_here(i) / h_here(i)
223 sf(kf(i),i) = hxs_here(i) / h_here(i)
227 h_here(i) = h(i,j,k)*gv%H_to_Z
228 hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
229 hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
231 h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
232 hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
233 hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
236 do i=is,ie ;
if (h_here(i) > 0.0)
then
237 hf(kf(i),i) = h_here(i)
238 tf(kf(i),i) = hxt_here(i) / h_here(i)
239 sf(kf(i),i) = hxs_here(i) / h_here(i)
242 do k=1,nz ;
do i=is,ie
243 if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i)))
then
244 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
248 h_here(i) = h(i,j,k)*gv%H_to_Z
249 hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
251 h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
252 hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
255 do i=is,ie ;
if (h_here(i) > 0.0)
then
256 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
262 do i=is,ie ;
if (cn(i,j)>0.0)
then
264 ig = i + g%idg_offset ; jg = j + g%jdg_offset
267 if (g%mask2dT(i,j) > 0.5)
then
269 lam = 1/(us%L_T_to_m_s**2 * cn(i,j)**2)
275 pres(k) = pres(k-1) + h_to_pres*hf(k-1,i)
276 t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
277 s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
280 kf(i)-1, tv%eqn_of_state, scale=us%kg_m3_to_R)
286 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
287 max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
288 drho_ds(k)*(sf(k,i)-sf(k-1,i)))
293 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
294 max(0.0,rf(k,i)-rf(k-1,i))
300 if (drxh_sum >= 0.0)
then
305 hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
307 if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
308 (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum)
then
310 i_hnew = 1.0 / (hc(kc) + hf(k,i))
311 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
312 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
313 hc(kc) = (hc(kc) + hf(k,i))
318 if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
319 (hc(k2) + hc(k2-1)) < tol2*drxh_sum)
then
321 i_hnew = 1.0 / (hc(kc) + hc(kc-1))
322 tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
323 sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
324 hc(kc-1) = (hc(kc) + hc(kc-1))
331 drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
332 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
337 gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
338 drho_ds(k)*(sc(k)-sc(k-1)))
343 hc(1) = hf(1,i) ; rc(1) = rf(1,i)
345 if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum)
then
347 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
348 hc(kc) = (hc(kc) + hf(k,i))
353 if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum)
then
355 rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
356 hc(kc-1) = (hc(kc) + hc(kc-1))
363 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
368 gprime(k) = g_rho0 * (rc(k)-rc(k-1))
380 if (kc >= modenum + 1)
then
386 igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
387 z_int(k) = z_int(k-1) + hc(k-1)
388 n2(k) = us%m_to_Z**2*gprime(k)/(0.5*(hc(k)+hc(k-1)))
391 n2(1) = n2(2) ; n2(kc+1) = n2(kc)
393 z_int(kc+1) = z_int(kc)+hc(kc)
395 if (abs(z_int(kc+1)-htot(i,j)) > 1.e-14*htot(i,j))
then
396 call mom_error(fatal,
"wave_structure: mismatch in total depths")
409 gp_unscaled = us%m_to_Z*gprime(k)
410 lam_z(row) = lam*gp_unscaled
411 a_diag(row) = gp_unscaled*(-igu(k))
412 b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
413 c_diag(row) = gp_unscaled*(-igl(k))
414 if (isnan(lam_z(row)))
then ; print *,
"Wave_structure: lam_z(row) is NAN" ;
endif
415 if (isnan(a_diag(row)))
then ; print *,
"Wave_structure: a(k) is NAN" ;
endif
416 if (isnan(b_diag(row)))
then ; print *,
"Wave_structure: b(k) is NAN" ;
endif
417 if (isnan(c_diag(row)))
then ; print *,
"Wave_structure: c(k) is NAN" ;
endif
421 gp_unscaled = us%m_to_Z*gprime(k)
422 lam_z(row) = lam*gp_unscaled
424 b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
425 c_diag(row) = gp_unscaled*(-igl(k))
428 gp_unscaled = us%m_to_Z*gprime(k)
429 lam_z(row) = lam*gp_unscaled
430 a_diag(row) = gp_unscaled*(-igu(k))
431 b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
435 e_guess(1:kc-1) = sin((z_int(2:kc)/htot(i,j)) *pi)
436 e_guess(1:kc-1) = e_guess(1:kc-1)/sqrt(sum(e_guess(1:kc-1)**2))
440 call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), &
441 -lam_z(1:kc-1),e_guess(1:kc-1),
"TDMA_H",e_itt)
442 e_guess(1:kc-1) = e_itt(1:kc-1) / sqrt(sum(e_itt(1:kc-1)**2))
444 w_strct(2:kc) = e_guess(1:kc-1)
449 ig_stop = 0 ; jg_stop = 0
450 if (isnan(sum(w_strct(1:kc+1))))
then
451 print *,
"Wave_structure: w_strct has a NAN at ig=", ig,
", jg=", jg
452 if (i<g%isc .or. i>g%iec .or. j<g%jsc .or. j>g%jec)
then
453 print *,
"This is occuring at a halo point."
455 ig_stop = ig ; jg_stop = jg
464 dz(k) = us%Z_to_m*hc(k)
465 w2avg = w2avg + 0.5*(w_strct(k)**2+w_strct(k+1)**2)*dz(k)
468 w2avg = w2avg / htot(i,j)
469 w_strct(:) = w_strct(:) / sqrt(htot(i,j)*w2avg*i_a_int)
473 u_strct(k) = 0.5*((w_strct(k-1) - w_strct(k) )/dz(k-1) + &
474 (w_strct(k) - w_strct(k+1))/dz(k))
476 u_strct(1) = (w_strct(1) - w_strct(2) )/dz(1)
477 u_strct(nzm) = (w_strct(nzm-1)- w_strct(nzm))/dz(nzm-1)
480 f2 = g%CoriolisBu(i,j)**2
483 kmag2 = us%m_to_L**2 * (freq**2 - f2) / (cn(i,j)**2 + cg_subro**2)
486 int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_n2w2 = 0.0
487 u_strct2(:) = u_strct(1:nzm)**2
488 w_strct2(:) = w_strct(1:nzm)**2
491 int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(k)+u_strct2(k+1)) * us%m_to_Z*dz(k)
492 int_w2 = int_w2 + 0.5*(w_strct2(k)+w_strct2(k+1)) * us%m_to_Z*dz(k)
493 int_n2w2 = int_n2w2 + 0.5*(w_strct2(k)*n2(k)+w_strct2(k+1)*n2(k+1)) * us%m_to_Z*dz(k)
497 if (
present(en) .and. (freq**2*kmag2 > 0.0))
then
499 ke_term = 0.25*gv%Rho0*( ((freq**2 + f2) / (freq**2*kmag2))*int_dwdz2 + int_w2 )
500 pe_term = 0.25*gv%Rho0*( int_n2w2 / (us%s_to_T*freq)**2 )
501 if (en(i,j) >= 0.0)
then
502 w0 = sqrt( en(i,j) / (ke_term + pe_term) )
504 call mom_error(warning,
"wave_structure: En < 0.0; setting to W0 to 0.0")
505 print *,
"En(i,j)=", en(i,j),
" at ig=", ig,
", jg=", jg
509 w_profile(:) = w0*w_strct(:)
512 uavg_profile(:) = us%Z_to_L*abs(w0*u_strct(:)) * sqrt((freq**2 + f2) / (2.0*freq**2*kmag2))
516 uavg_profile(:) = 0.0
520 cs%w_strct(i,j,1:nzm) = w_strct(:)
521 cs%u_strct(i,j,1:nzm) = u_strct(:)
522 cs%W_profile(i,j,1:nzm) = w_profile(:)
523 cs%Uavg_profile(i,j,1:nzm)= uavg_profile(:)
524 cs%z_depths(i,j,1:nzm) = us%Z_to_m*z_int(:)
525 cs%N2(i,j,1:nzm) = n2(:)
526 cs%num_intfaces(i,j) = nzm
530 cs%w_strct(i,j,1:nzm) = 0.0
531 cs%u_strct(i,j,1:nzm) = 0.0
532 cs%W_profile(i,j,1:nzm) = 0.0
533 cs%Uavg_profile(i,j,1:nzm)= 0.0
534 cs%z_depths(i,j,1:nzm) = 0.0
535 cs%N2(i,j,1:nzm) = 0.0
536 cs%num_intfaces(i,j) = nzm
546 cs%w_strct(i,j,1:nzm) = 0.0
547 cs%u_strct(i,j,1:nzm) = 0.0
548 cs%W_profile(i,j,1:nzm) = 0.0
549 cs%Uavg_profile(i,j,1:nzm)= 0.0
550 cs%z_depths(i,j,1:nzm) = 0.0
551 cs%N2(i,j,1:nzm) = 0.0
552 cs%num_intfaces(i,j) = nzm
562 real,
dimension(:),
intent(in) :: a
563 real,
dimension(:),
intent(in) :: b
564 real,
dimension(:),
intent(in) :: c
565 real,
dimension(:),
intent(in) :: h
572 real,
dimension(:),
intent(in) :: y
573 character(len=*),
intent(in) :: method
574 real,
dimension(:),
intent(out) :: x
579 real,
allocatable,
dimension(:) :: c_prime, y_prime, q, alpha
581 real :: Q_prime, beta
586 allocate(c_prime(nrow))
587 allocate(y_prime(nrow))
589 allocate(alpha(nrow))
593 if (method ==
'TDMA_T')
then
596 c_prime(:) = 0.0 ; y_prime(:) = 0.0
597 c_prime(1) = c(1)/b(1) ; y_prime(1) = y(1)/b(1)
601 c_prime(k) = c(k)/(b(k)-a(k)*c_prime(k-1))
605 y_prime(k) = (y(k)-a(k)*y_prime(k-1))/(b(k)-a(k)*c_prime(k-1))
608 x(nrow) = y_prime(nrow)
612 x(k) = y_prime(k)-c_prime(k)*x(k+1)
631 elseif (method ==
'TDMA_H')
then
641 if (abs(a(k+1)-c(k)) > 1.e-10*(abs(a(k+1))+abs(c(k))))
then
642 call mom_error(warning,
"tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H")
648 alpha(nrow) = b(nrow)-h(nrow)-alpha(nrow-1)
651 y_prime(:) = 0.0 ; q(:) = 0.0
652 y_prime(1) = beta*y(1) ; q(1) = beta*alpha(1)
657 beta = 1/(h(k)+alpha(k-1)*q_prime+alpha(k))
658 if (isnan(beta))
then ; print *,
"Tridiag_solver: beta is NAN" ;
endif
660 y_prime(k) = beta*(y(k)+alpha(k-1)*y_prime(k-1))
661 q_prime = beta*(h(k)+alpha(k-1)*q_prime)
663 if ((h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow)) == 0.0)
then
664 call mom_error(fatal,
"Tridiag_solver: this system is not stable.")
667 beta = 1/(h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow))
669 y_prime(nrow) = beta*(y(nrow)+alpha(nrow-1)*y_prime(nrow-1))
670 x(nrow) = y_prime(nrow)
673 x(k) = y_prime(k)+q(k)*x(k+1)
679 deallocate(c_prime,y_prime,q,alpha)
686 type(time_type),
intent(in) :: time
690 type(
diag_ctrl),
target,
intent(in) :: diag
695 #include "version_variable.h"
696 character(len=40) :: mdl =
"MOM_wave_structure"
697 integer :: isd, ied, jsd, jed, nz
699 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
701 if (
associated(cs))
then
702 call mom_error(warning,
"wave_structure_init called with an "// &
703 "associated control structure.")
705 else ;
allocate(cs) ;
endif
707 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
708 "X Location of generation site for internal tide", default=1.)
709 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
710 "Y Location of generation site for internal tide", default=1.)
716 allocate(cs%w_strct(isd:ied,jsd:jed,nz+1))
717 allocate(cs%u_strct(isd:ied,jsd:jed,nz+1))
718 allocate(cs%W_profile(isd:ied,jsd:jed,nz+1))
719 allocate(cs%Uavg_profile(isd:ied,jsd:jed,nz+1))
720 allocate(cs%z_depths(isd:ied,jsd:jed,nz+1))
721 allocate(cs%N2(isd:ied,jsd:jed,nz+1))
722 allocate(cs%num_intfaces(isd:ied,jsd:jed))