Calls subroutines in this file that are needed to refract, propagate, and dissipate energy density of the internal tide.
152 type(ocean_grid_type),
intent(inout) :: G
153 type(verticalGrid_type),
intent(in) :: GV
154 type(unit_scale_type),
intent(in) :: US
155 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
157 type(thermo_var_ptrs),
intent(in) :: tv
159 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: TKE_itidal_input
161 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: vel_btTide
163 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: Nb
164 real,
intent(in) :: dt
166 type(int_tide_CS),
pointer :: CS
168 real,
dimension(SZI_(G),SZJ_(G),CS%nMode), &
171 real,
dimension(SZI_(G),SZJ_(G),2) :: &
173 real,
dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
176 real,
dimension(SZI_(G),SZJB_(G)) :: &
179 real,
dimension(SZI_(G),SZJ_(G)) :: &
181 tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_Froude_loss, tot_allprocesses_loss, &
184 itidal_loss_mode, allprocesses_loss_mode
186 real :: frac_per_sector, f2, I_rho0, I_D_here, freq2, Kmag2
187 real :: c_phase, loss_rate, Fr2_max
188 real,
parameter :: cn_subRO = 1e-100
189 real :: En_new, En_check
190 real :: En_initial, Delta_E_check
191 real :: TKE_Froude_loss_check, TKE_Froude_loss_tot
192 character(len=160) :: mesg
193 integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm
194 integer :: id_g, jd_g
195 type(group_pass_type),
save :: pass_test, pass_En
197 if (.not.
associated(cs))
return
198 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
199 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nangle = cs%NAngle
200 i_rho0 = 1.0 / gv%Rho0
210 if (cs%energized_angle <= 0)
then
211 frac_per_sector = 1.0 / real(cs%nAngle * cs%nMode * cs%nFreq)
212 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
213 f2 = 0.25*us%s_to_T**2*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
214 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
215 if (cs%frequency(fr)**2 > f2) &
216 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + &
217 dt*frac_per_sector*(1-cs%q_itides)*tke_itidal_input(i,j)
218 enddo ;
enddo ;
enddo ;
enddo ;
enddo
219 elseif (cs%energized_angle <= cs%nAngle)
then
220 frac_per_sector = 1.0 / real(cs%nMode * cs%nFreq)
221 a = cs%energized_angle
222 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do j=js,je ;
do i=is,ie
223 f2 = 0.25*us%s_to_T**2*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
224 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
225 if (cs%frequency(fr)**2 > f2) &
226 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + &
227 dt*frac_per_sector**(1-cs%q_itides)*tke_itidal_input(i,j)
228 enddo ;
enddo ;
enddo ;
enddo
230 call mom_error(warning,
"Internal tide energy is being put into a angular "//&
231 "band that does not exist.")
235 do j=jsd,jed ;
do i=isd,ied ; test(i,j,1) = 1.0 ; test(i,j,2) = 0.0 ;
enddo ;
enddo
236 do m=1,cs%nMode ;
do fr=1,cs%nFreq
237 call create_group_pass(pass_en, cs%En(:,:,:,fr,m), g%domain)
239 call create_group_pass(pass_test, test(:,:,1), test(:,:,2), g%domain, stagger=agrid)
240 call start_group_pass(pass_test, g%domain)
243 do m=1,cs%nMode ;
do fr=1,cs%nFreq
244 call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, g, us, cs%nAngle, cs%use_PPMang)
248 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
249 do j=js,je ;
do i=is,ie
250 if (cs%En(i,j,a,fr,m)<0.0)
then
251 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
252 write(mesg,*)
'After first refraction: En<0.0 at ig=', id_g,
', jg=', jd_g, &
253 'En=',cs%En(i,j,a,fr,m)
254 call mom_error(warning,
"propagate_int_tide: "//trim(mesg))
255 cs%En(i,j,a,fr,m) = 0.0
259 enddo ;
enddo ;
enddo
261 call do_group_pass(pass_en, g%domain)
263 call complete_group_pass(pass_test, g%domain)
266 call correct_halo_rotation(cs%En, test, g, cs%nAngle)
269 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
270 call propagate(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), dt, g, us, cs, cs%NAngle)
274 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
275 do j=js,je ;
do i=is,ie
276 if (cs%En(i,j,a,fr,m)<0.0)
then
277 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
278 if (abs(cs%En(i,j,a,fr,m))>1.0)
then
279 write(mesg,*)
'After propagation: En<0.0 at ig=', id_g,
', jg=', jd_g, &
280 'En=', cs%En(i,j,a,fr,m)
281 call mom_error(warning,
"propagate_int_tide: "//trim(mesg))
284 cs%En(i,j,a,fr,m) = 0.0
287 enddo ;
enddo ;
enddo
290 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
291 call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, g, us, cs%NAngle, cs%use_PPMang)
295 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
296 do j=js,je ;
do i=is,ie
297 if (cs%En(i,j,a,fr,m)<0.0)
then
298 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
299 write(mesg,*)
'After second refraction: En<0.0 at ig=', id_g,
', jg=', jd_g, &
300 'En=',cs%En(i,j,a,fr,m)
301 call mom_error(warning,
"propagate_int_tide: "//trim(mesg))
302 cs%En(i,j,a,fr,m) = 0.0
306 enddo ;
enddo ;
enddo
309 if (cs%apply_background_drag .or. cs%apply_bottom_drag &
310 .or. cs%apply_wave_drag .or. cs%apply_Froude_drag &
311 .or. (cs%id_tot_En > 0))
then
313 tot_en_mode(:,:,:,:) = 0.0
314 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
315 do j=jsd,jed ;
do i=isd,ied ;
do a=1,cs%nAngle
316 tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
317 tot_en_mode(i,j,fr,m) = tot_en_mode(i,j,fr,m) + cs%En(i,j,a,fr,m)
318 enddo ;
enddo ;
enddo
323 if (cs%apply_background_drag)
then
324 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do a=1,cs%nAngle ;
do j=jsd,jed ;
do i=isd,ied
327 cs%TKE_leak_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * cs%decay_rate
328 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt *cs%decay_rate)
329 enddo ;
enddo ;
enddo ;
enddo ;
enddo
332 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
333 do j=js,je ;
do i=is,ie
334 if (cs%En(i,j,a,fr,m)<0.0)
then
335 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
336 write(mesg,*)
'After leak loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
337 'En=',cs%En(i,j,a,fr,m)
338 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
339 cs%En(i,j,a,fr,m) = 0.0
343 enddo ;
enddo ;
enddo
346 if (cs%apply_bottom_drag)
then
347 do j=jsd,jed ;
do i=isd,ied
349 i_d_here = 1.0 / (us%Z_to_m*max(g%bathyT(i,j), 1.0*us%m_to_Z))
350 drag_scale(i,j) = cs%cdrag * sqrt(max(0.0, vel_bttide(i,j)**2 + &
351 tot_en(i,j) * i_rho0 * i_d_here)) * i_d_here
353 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do a=1,cs%nAngle ;
do j=jsd,jed ;
do i=isd,ied
356 cs%TKE_quad_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * drag_scale(i,j)
357 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt *drag_scale(i,j))
358 enddo ;
enddo ;
enddo ;
enddo ;
enddo
361 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
362 do j=js,je ;
do i=is,ie
363 if (cs%En(i,j,a,fr,m)<0.0)
then
364 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
365 write(mesg,*)
'After bottom loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
366 'En=',cs%En(i,j,a,fr,m)
367 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
368 cs%En(i,j,a,fr,m) = 0.0
373 enddo ;
enddo ;
enddo
378 if (cs%apply_wave_drag .or. cs%apply_Froude_drag)
then
379 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
381 call wave_structure(h, tv, g, gv, us, cn(:,:,m), m, cs%frequency(fr), &
382 cs%wave_structure_CSp, tot_en_mode(:,:,fr,m), full_halos=.true.)
384 do j=jsd,jed ;
do i=isd,ied
385 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
386 nzm = cs%wave_structure_CSp%num_intfaces(i,j)
387 ub(i,j,fr,m) = cs%wave_structure_CSp%Uavg_profile(i,j,nzm)
388 umax(i,j,fr,m) = maxval(cs%wave_structure_CSp%Uavg_profile(i,j,1:nzm))
393 if (cs%apply_wave_drag)
then
395 call itidal_lowmode_loss(g, us, cs, nb, ub, cs%En, cs%TKE_itidal_loss_fixed, &
396 cs%TKE_itidal_loss, dt, full_halos=.false.)
399 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
400 do j=js,je ;
do i=is,ie
401 if (cs%En(i,j,a,fr,m)<0.0)
then
402 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
403 write(mesg,*)
'After wave drag loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
404 'En=',cs%En(i,j,a,fr,m)
405 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
406 cs%En(i,j,a,fr,m) = 0.0
410 enddo ;
enddo ;
enddo
413 if (cs%apply_Froude_drag)
then
415 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
416 freq2 = cs%frequency(fr)**2
417 do j=jsd,jed ;
do i=isd,ied
418 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
420 f2 = 0.25*us%s_to_T**2*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
421 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
422 kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subro**2)
424 if (kmag2 > 0.0)
then
425 c_phase = sqrt(freq2/kmag2)
426 nzm = cs%wave_structure_CSp%num_intfaces(i,j)
427 fr2_max = (umax(i,j,fr,m)/c_phase)**2
429 if (fr2_max > 1.0)
then
430 en_initial = sum(cs%En(i,j,:,fr,m))
432 loss_rate = (1/fr2_max - 1.0)/dt
435 cs%TKE_Froude_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * abs(loss_rate)
437 en_new = cs%En(i,j,a,fr,m)/fr2_max
438 en_check = cs%En(i,j,a,fr,m) - cs%TKE_Froude_loss(i,j,a,fr,m)*dt
440 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m)/fr2_max
442 if (abs(en_new - en_check) > 1e-10)
then
443 call mom_error(warning,
"MOM_internal_tides: something is wrong with Fr-breaking.", &
445 write(mesg,*)
"En_new=", en_new ,
"En_check=", en_check
446 call mom_error(warning,
"MOM_internal_tides: "//trim(mesg), all_print=.true.)
450 delta_e_check = en_initial - sum(cs%En(i,j,:,fr,m))
451 tke_froude_loss_check = abs(delta_e_check)/dt
452 tke_froude_loss_tot = sum(cs%TKE_Froude_loss(i,j,:,fr,m))
453 if (abs(tke_froude_loss_check - tke_froude_loss_tot) > 1e-10)
then
454 call mom_error(warning,
"MOM_internal_tides: something is wrong with Fr energy update.", &
456 write(mesg,*)
"TKE_Froude_loss_check=", tke_froude_loss_check, &
457 "TKE_Froude_loss_tot=", tke_froude_loss_tot
458 call mom_error(warning,
"MOM_internal_tides: "//trim(mesg), all_print=.true.)
462 cs%cp(i,j,fr,m) = c_phase
467 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
468 do j=js,je ;
do i=is,ie
469 if (cs%En(i,j,a,fr,m)<0.0)
then
470 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
471 write(mesg,*)
'After Froude loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
472 'En=',cs%En(i,j,a,fr,m)
473 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
474 cs%En(i,j,a,fr,m) = 0.0
479 enddo ;
enddo ;
enddo
482 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
483 call sum_en(g,cs,cs%En(:,:,:,fr,m),
'prop_int_tide')
487 if (query_averaging_enabled(cs%diag))
then
489 if (cs%id_tot_En > 0)
call post_data(cs%id_tot_En, tot_en, cs%diag)
490 if (cs%id_itide_drag > 0)
call post_data(cs%id_itide_drag, drag_scale, cs%diag)
491 if (cs%id_TKE_itidal_input > 0)
call post_data(cs%id_TKE_itidal_input, &
492 tke_itidal_input, cs%diag)
495 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_En_mode(fr,m) > 0)
then
497 do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
498 tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
499 enddo ;
enddo ;
enddo
500 call post_data(cs%id_En_mode(fr,m), tot_en, cs%diag)
501 endif ;
enddo ;
enddo
504 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_En_ang_mode(fr,m) > 0)
then
505 call post_data(cs%id_En_ang_mode(fr,m), cs%En(:,:,:,fr,m) , cs%diag)
506 endif ;
enddo ;
enddo
509 tot_leak_loss(:,:) = 0.0
510 tot_quad_loss(:,:) = 0.0
511 tot_itidal_loss(:,:) = 0.0
512 tot_froude_loss(:,:) = 0.0
513 tot_allprocesses_loss(:,:) = 0.0
514 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
515 tot_leak_loss(i,j) = tot_leak_loss(i,j) + cs%TKE_leak_loss(i,j,a,fr,m)
516 tot_quad_loss(i,j) = tot_quad_loss(i,j) + cs%TKE_quad_loss(i,j,a,fr,m)
517 tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
518 tot_froude_loss(i,j) = tot_froude_loss(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
519 enddo ;
enddo ;
enddo ;
enddo ;
enddo
520 do j=js,je ;
do i=is,ie
521 tot_allprocesses_loss(i,j) = tot_leak_loss(i,j) + tot_quad_loss(i,j) + &
522 tot_itidal_loss(i,j) + tot_froude_loss(i,j)
524 cs%tot_leak_loss = tot_leak_loss
525 cs%tot_quad_loss = tot_quad_loss
526 cs%tot_itidal_loss = tot_itidal_loss
527 cs%tot_Froude_loss = tot_froude_loss
528 cs%tot_allprocesses_loss = tot_allprocesses_loss
529 if (cs%id_tot_leak_loss > 0)
then
530 call post_data(cs%id_tot_leak_loss, tot_leak_loss, cs%diag)
532 if (cs%id_tot_quad_loss > 0)
then
533 call post_data(cs%id_tot_quad_loss, tot_quad_loss, cs%diag)
535 if (cs%id_tot_itidal_loss > 0)
then
536 call post_data(cs%id_tot_itidal_loss, tot_itidal_loss, cs%diag)
538 if (cs%id_tot_Froude_loss > 0)
then
539 call post_data(cs%id_tot_Froude_loss, tot_froude_loss, cs%diag)
541 if (cs%id_tot_allprocesses_loss > 0)
then
542 call post_data(cs%id_tot_allprocesses_loss, tot_allprocesses_loss, cs%diag)
546 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
547 if (cs%id_itidal_loss_mode(fr,m) > 0 .or. cs%id_allprocesses_loss_mode(fr,m) > 0)
then
548 itidal_loss_mode(:,:) = 0.0
549 allprocesses_loss_mode(:,:) = 0.0
550 do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
551 itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
552 allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + &
553 cs%TKE_leak_loss(i,j,a,fr,m) + cs%TKE_quad_loss(i,j,a,fr,m) + &
554 cs%TKE_itidal_loss(i,j,a,fr,m) + cs%TKE_Froude_loss(i,j,a,fr,m)
555 enddo ;
enddo ;
enddo
556 call post_data(cs%id_itidal_loss_mode(fr,m), itidal_loss_mode, cs%diag)
557 call post_data(cs%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, cs%diag)
558 endif ;
enddo ;
enddo
561 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_itidal_loss_ang_mode(fr,m) > 0)
then
562 call post_data(cs%id_itidal_loss_ang_mode(fr,m), cs%TKE_itidal_loss(:,:,:,fr,m) , cs%diag)
563 endif ;
enddo ;
enddo
566 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_Ub_mode(fr,m) > 0)
then
567 call post_data(cs%id_Ub_mode(fr,m), ub(:,:,fr,m), cs%diag)
568 endif ;
enddo ;
enddo
571 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_cp_mode(fr,m) > 0)
then
572 call post_data(cs%id_cp_mode(fr,m), cs%cp(:,:,fr,m), cs%diag)
573 endif ;
enddo ;
enddo