21 use mom_time_manager,
only : time_type, time_type_to_real,
operator(+),
operator(/),
operator(-)
29 implicit none ;
private
31 #include <MOM_memory.h>
39 logical :: do_int_tides
42 integer :: nangle = 24
43 integer :: energized_angle = -1
53 real,
allocatable,
dimension(:,:) :: refl_angle
56 real :: nullangle = -999.9
57 real,
allocatable,
dimension(:,:) :: refl_pref
60 logical,
allocatable,
dimension(:,:) :: refl_pref_logical
63 logical,
allocatable,
dimension(:,:) :: refl_dbl
67 real,
allocatable,
dimension(:,:,:,:) :: cp
69 real,
allocatable,
dimension(:,:,:,:,:) :: tke_leak_loss
71 real,
allocatable,
dimension(:,:,:,:,:) :: tke_quad_loss
73 real,
allocatable,
dimension(:,:,:,:,:) :: tke_froude_loss
75 real,
allocatable,
dimension(:,:) :: tke_itidal_loss_fixed
79 real,
allocatable,
dimension(:,:,:,:,:) :: tke_itidal_loss
81 real,
allocatable,
dimension(:,:) :: tot_leak_loss
83 real,
allocatable,
dimension(:,:) :: tot_quad_loss
85 real,
allocatable,
dimension(:,:) :: tot_itidal_loss
87 real,
allocatable,
dimension(:,:) :: tot_froude_loss
89 real,
allocatable,
dimension(:,:) :: tot_allprocesses_loss
93 type(time_type),
pointer :: time => null()
94 character(len=200) :: inputdir
98 logical :: apply_background_drag
100 logical :: apply_bottom_drag
102 logical :: apply_wave_drag
104 logical :: apply_froude_drag
106 real,
dimension(:,:,:,:,:),
pointer :: en => null()
109 real,
dimension(:,:,:),
pointer :: en_restart => null()
111 real,
allocatable,
dimension(:) :: frequency
120 integer :: id_tot_en = -1, id_tke_itidal_input = -1, id_itide_drag = -1
121 integer :: id_refl_pref = -1, id_refl_ang = -1, id_land_mask = -1
122 integer :: id_dx_cv = -1, id_dy_cu = -1
124 integer :: id_tot_leak_loss = -1, id_tot_quad_loss = -1, id_tot_itidal_loss = -1
125 integer :: id_tot_froude_loss = -1, id_tot_allprocesses_loss = -1
127 integer,
allocatable,
dimension(:,:) :: &
129 id_itidal_loss_mode, &
130 id_allprocesses_loss_mode, &
134 integer,
allocatable,
dimension(:,:) :: &
136 id_itidal_loss_ang_mode
144 integer :: ish, ieh, jsh, jeh
157 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
161 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: tke_itidal_input
163 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: vel_bttide
165 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: nb
166 real,
intent(in) :: dt
170 real,
dimension(SZI_(G),SZJ_(G),CS%nMode), &
174 real,
dimension(SZI_(G),SZJ_(G),2) :: &
176 real,
dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
180 real,
dimension(SZI_(G),SZJB_(G)) :: &
183 real,
dimension(SZI_(G),SZJ_(G)) :: &
185 tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_froude_loss, tot_allprocesses_loss, &
188 itidal_loss_mode, allprocesses_loss_mode
190 real :: frac_per_sector, f2, kmag2
198 real :: en_new, en_check
199 real :: en_initial, delta_e_check
200 real :: tke_froude_loss_check, tke_froude_loss_tot
201 character(len=160) :: mesg
202 integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nangle, nzm
203 integer :: id_g, jd_g
204 type(group_pass_type),
save :: pass_test, pass_en
206 if (.not.
associated(cs))
return
207 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
208 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nangle = cs%NAngle
209 i_rho0 = 1.0 / gv%Rho0
210 cn_subro = 1e-100*us%m_s_to_L_T
220 if (cs%energized_angle <= 0)
then
221 frac_per_sector = 1.0 / real(cs%nAngle * cs%nMode * cs%nFreq)
222 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
223 f2 = 0.25*((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) + dt*frac_per_sector*(1.0-cs%q_itides) * &
227 tke_itidal_input(i,j)
228 enddo ;
enddo ;
enddo ;
enddo ;
enddo
229 elseif (cs%energized_angle <= cs%nAngle)
then
230 frac_per_sector = 1.0 / real(cs%nMode * cs%nFreq)
231 a = cs%energized_angle
232 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do j=js,je ;
do i=is,ie
233 f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
234 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
235 if (cs%frequency(fr)**2 > f2) &
236 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-cs%q_itides) * &
237 tke_itidal_input(i,j)
238 enddo ;
enddo ;
enddo ;
enddo
240 call mom_error(warning,
"Internal tide energy is being put into a angular "//&
241 "band that does not exist.")
245 do j=jsd,jed ;
do i=isd,ied ; test(i,j,1) = 1.0 ; test(i,j,2) = 0.0 ;
enddo ;
enddo
246 do m=1,cs%nMode ;
do fr=1,cs%nFreq
249 call create_group_pass(pass_test, test(:,:,1), test(:,:,2), g%domain, stagger=agrid)
253 do m=1,cs%nMode ;
do fr=1,cs%nFreq
254 call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, &
255 g, us, cs%nAngle, cs%use_PPMang)
259 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
260 do j=js,je ;
do i=is,ie
261 if (cs%En(i,j,a,fr,m)<0.0)
then
262 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
263 write(mesg,*)
'After first refraction: En<0.0 at ig=', id_g,
', jg=', jd_g, &
264 'En=',cs%En(i,j,a,fr,m)
265 call mom_error(warning,
"propagate_int_tide: "//trim(mesg))
266 cs%En(i,j,a,fr,m) = 0.0
270 enddo ;
enddo ;
enddo
272 call do_group_pass(pass_en, g%domain)
280 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
281 call propagate(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), dt, &
282 g, us, cs, cs%NAngle)
286 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
287 do j=js,je ;
do i=is,ie
288 if (cs%En(i,j,a,fr,m)<0.0)
then
289 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
290 if (abs(cs%En(i,j,a,fr,m))>1.0)
then
291 write(mesg,*)
'After propagation: En<0.0 at ig=', id_g,
', jg=', jd_g, &
292 'En=', cs%En(i,j,a,fr,m)
293 call mom_error(warning,
"propagate_int_tide: "//trim(mesg))
296 cs%En(i,j,a,fr,m) = 0.0
299 enddo ;
enddo ;
enddo
302 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
303 call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, &
304 g, us, cs%NAngle, cs%use_PPMang)
308 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
309 do j=js,je ;
do i=is,ie
310 if (cs%En(i,j,a,fr,m)<0.0)
then
311 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
312 write(mesg,*)
'After second refraction: En<0.0 at ig=', id_g,
', jg=', jd_g, &
313 'En=',cs%En(i,j,a,fr,m)
314 call mom_error(warning,
"propagate_int_tide: "//trim(mesg))
315 cs%En(i,j,a,fr,m) = 0.0
319 enddo ;
enddo ;
enddo
322 if (cs%apply_background_drag .or. cs%apply_bottom_drag &
323 .or. cs%apply_wave_drag .or. cs%apply_Froude_drag &
324 .or. (cs%id_tot_En > 0))
then
326 tot_en_mode(:,:,:,:) = 0.0
327 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
328 do j=jsd,jed ;
do i=isd,ied ;
do a=1,cs%nAngle
329 tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
330 tot_en_mode(i,j,fr,m) = tot_en_mode(i,j,fr,m) + cs%En(i,j,a,fr,m)
331 enddo ;
enddo ;
enddo
336 if (cs%apply_background_drag)
then
337 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do a=1,cs%nAngle ;
do j=jsd,jed ;
do i=isd,ied
340 cs%TKE_leak_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * cs%decay_rate
341 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt * cs%decay_rate)
342 enddo ;
enddo ;
enddo ;
enddo ;
enddo
345 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
346 do j=js,je ;
do i=is,ie
347 if (cs%En(i,j,a,fr,m)<0.0)
then
348 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
349 write(mesg,*)
'After leak loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
350 'En=',cs%En(i,j,a,fr,m)
351 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
352 cs%En(i,j,a,fr,m) = 0.0
356 enddo ;
enddo ;
enddo
359 if (cs%apply_bottom_drag)
then
360 do j=jsd,jed ;
do i=isd,ied
362 i_d_here = 1.0 / (max(g%bathyT(i,j), 1.0*us%m_to_Z))
363 drag_scale(i,j) = cs%cdrag * sqrt(max(0.0, us%L_to_Z**2*vel_bttide(i,j)**2 + &
364 tot_en(i,j) * i_rho0 * i_d_here)) * i_d_here
366 do m=1,cs%nMode ;
do fr=1,cs%nFreq ;
do a=1,cs%nAngle ;
do j=jsd,jed ;
do i=isd,ied
369 cs%TKE_quad_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * drag_scale(i,j)
370 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt * drag_scale(i,j))
371 enddo ;
enddo ;
enddo ;
enddo ;
enddo
374 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
375 do j=js,je ;
do i=is,ie
376 if (cs%En(i,j,a,fr,m)<0.0)
then
377 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
378 write(mesg,*)
'After bottom loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
379 'En=',cs%En(i,j,a,fr,m)
380 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
381 cs%En(i,j,a,fr,m) = 0.0
386 enddo ;
enddo ;
enddo
391 if (cs%apply_wave_drag .or. cs%apply_Froude_drag)
then
392 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
394 call wave_structure(h, tv, g, gv, us, cn(:,:,m), m, cs%frequency(fr), &
395 cs%wave_structure_CSp, tot_en_mode(:,:,fr,m), full_halos=.true.)
397 do j=jsd,jed ;
do i=isd,ied
398 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
399 nzm = cs%wave_structure_CSp%num_intfaces(i,j)
400 ub(i,j,fr,m) = cs%wave_structure_CSp%Uavg_profile(i,j,nzm)
401 umax(i,j,fr,m) = maxval(cs%wave_structure_CSp%Uavg_profile(i,j,1:nzm))
406 if (cs%apply_wave_drag)
then
409 cs%TKE_itidal_loss, dt, full_halos=.false.)
412 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
413 do j=js,je ;
do i=is,ie
414 if (cs%En(i,j,a,fr,m)<0.0)
then
415 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
416 write(mesg,*)
'After wave drag loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
417 'En=',cs%En(i,j,a,fr,m)
418 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
419 cs%En(i,j,a,fr,m) = 0.0
423 enddo ;
enddo ;
enddo
426 if (cs%apply_Froude_drag)
then
428 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
429 freq2 = cs%frequency(fr)**2
430 do j=jsd,jed ;
do i=isd,ied
431 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
433 f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
434 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
435 kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subro**2)
437 if (kmag2 > 0.0)
then
438 c_phase = sqrt(freq2/kmag2)
439 nzm = cs%wave_structure_CSp%num_intfaces(i,j)
440 fr2_max = (umax(i,j,fr,m) / c_phase)**2
442 if (fr2_max > 1.0)
then
443 en_initial = sum(cs%En(i,j,:,fr,m))
445 loss_rate = (1.0 - fr2_max) / (fr2_max * dt)
448 cs%TKE_Froude_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * abs(loss_rate)
450 en_new = cs%En(i,j,a,fr,m)/fr2_max
451 en_check = cs%En(i,j,a,fr,m) - cs%TKE_Froude_loss(i,j,a,fr,m)*dt
453 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m)/fr2_max
455 if (abs(en_new - en_check) > 1e-10*us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**2)
then
456 call mom_error(warning,
"MOM_internal_tides: something is wrong with Fr-breaking.", &
458 write(mesg,*)
"En_new=", en_new ,
"En_check=", en_check
459 call mom_error(warning,
"MOM_internal_tides: "//trim(mesg), all_print=.true.)
463 delta_e_check = en_initial - sum(cs%En(i,j,:,fr,m))
464 tke_froude_loss_check = abs(delta_e_check)/dt
465 tke_froude_loss_tot = sum(cs%TKE_Froude_loss(i,j,:,fr,m))
466 if (abs(tke_froude_loss_check - tke_froude_loss_tot) > 1e-10)
then
467 call mom_error(warning,
"MOM_internal_tides: something is wrong with Fr energy update.", &
469 write(mesg,*)
"TKE_Froude_loss_check=", tke_froude_loss_check, &
470 "TKE_Froude_loss_tot=", tke_froude_loss_tot
471 call mom_error(warning,
"MOM_internal_tides: "//trim(mesg), all_print=.true.)
475 cs%cp(i,j,fr,m) = c_phase
480 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle
481 do j=js,je ;
do i=is,ie
482 if (cs%En(i,j,a,fr,m)<0.0)
then
483 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
484 write(mesg,*)
'After Froude loss: En<0.0 at ig=', id_g,
', jg=', jd_g, &
485 'En=',cs%En(i,j,a,fr,m)
486 call mom_error(warning,
"propagate_int_tide: "//trim(mesg), all_print=.true.)
487 cs%En(i,j,a,fr,m) = 0.0
492 enddo ;
enddo ;
enddo
495 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
496 call sum_en(g,cs,cs%En(:,:,:,fr,m),
'prop_int_tide')
500 if (query_averaging_enabled(cs%diag))
then
502 if (cs%id_tot_En > 0)
call post_data(cs%id_tot_En, tot_en, cs%diag)
503 if (cs%id_itide_drag > 0)
call post_data(cs%id_itide_drag, drag_scale, cs%diag)
504 if (cs%id_TKE_itidal_input > 0)
call post_data(cs%id_TKE_itidal_input, &
505 tke_itidal_input, cs%diag)
508 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_En_mode(fr,m) > 0)
then
510 do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
511 tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
512 enddo ;
enddo ;
enddo
513 call post_data(cs%id_En_mode(fr,m), tot_en, cs%diag)
514 endif ;
enddo ;
enddo
517 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_En_ang_mode(fr,m) > 0)
then
518 call post_data(cs%id_En_ang_mode(fr,m), cs%En(:,:,:,fr,m) , cs%diag)
519 endif ;
enddo ;
enddo
522 tot_leak_loss(:,:) = 0.0
523 tot_quad_loss(:,:) = 0.0
524 tot_itidal_loss(:,:) = 0.0
525 tot_froude_loss(:,:) = 0.0
526 tot_allprocesses_loss(:,:) = 0.0
527 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
528 tot_leak_loss(i,j) = tot_leak_loss(i,j) + cs%TKE_leak_loss(i,j,a,fr,m)
529 tot_quad_loss(i,j) = tot_quad_loss(i,j) + cs%TKE_quad_loss(i,j,a,fr,m)
530 tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
531 tot_froude_loss(i,j) = tot_froude_loss(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
532 enddo ;
enddo ;
enddo ;
enddo ;
enddo
533 do j=js,je ;
do i=is,ie
534 tot_allprocesses_loss(i,j) = tot_leak_loss(i,j) + tot_quad_loss(i,j) + &
535 tot_itidal_loss(i,j) + tot_froude_loss(i,j)
537 cs%tot_leak_loss = tot_leak_loss
538 cs%tot_quad_loss = tot_quad_loss
539 cs%tot_itidal_loss = tot_itidal_loss
540 cs%tot_Froude_loss = tot_froude_loss
541 cs%tot_allprocesses_loss = tot_allprocesses_loss
542 if (cs%id_tot_leak_loss > 0)
then
543 call post_data(cs%id_tot_leak_loss, tot_leak_loss, cs%diag)
545 if (cs%id_tot_quad_loss > 0)
then
546 call post_data(cs%id_tot_quad_loss, tot_quad_loss, cs%diag)
548 if (cs%id_tot_itidal_loss > 0)
then
549 call post_data(cs%id_tot_itidal_loss, tot_itidal_loss, cs%diag)
551 if (cs%id_tot_Froude_loss > 0)
then
552 call post_data(cs%id_tot_Froude_loss, tot_froude_loss, cs%diag)
554 if (cs%id_tot_allprocesses_loss > 0)
then
555 call post_data(cs%id_tot_allprocesses_loss, tot_allprocesses_loss, cs%diag)
559 do m=1,cs%NMode ;
do fr=1,cs%Nfreq
560 if (cs%id_itidal_loss_mode(fr,m) > 0 .or. cs%id_allprocesses_loss_mode(fr,m) > 0)
then
561 itidal_loss_mode(:,:) = 0.0
562 allprocesses_loss_mode(:,:) = 0.0
563 do a=1,cs%nAngle ;
do j=js,je ;
do i=is,ie
564 itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
565 allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + &
566 cs%TKE_leak_loss(i,j,a,fr,m) + cs%TKE_quad_loss(i,j,a,fr,m) + &
567 cs%TKE_itidal_loss(i,j,a,fr,m) + cs%TKE_Froude_loss(i,j,a,fr,m)
568 enddo ;
enddo ;
enddo
569 call post_data(cs%id_itidal_loss_mode(fr,m), itidal_loss_mode, cs%diag)
570 call post_data(cs%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, cs%diag)
571 endif ;
enddo ;
enddo
574 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_itidal_loss_ang_mode(fr,m) > 0)
then
575 call post_data(cs%id_itidal_loss_ang_mode(fr,m), cs%TKE_itidal_loss(:,:,:,fr,m) , cs%diag)
576 endif ;
enddo ;
enddo
579 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_Ub_mode(fr,m) > 0)
then
580 call post_data(cs%id_Ub_mode(fr,m), ub(:,:,fr,m), cs%diag)
581 endif ;
enddo ;
enddo
584 do m=1,cs%NMode ;
do fr=1,cs%Nfreq ;
if (cs%id_cp_mode(fr,m) > 0)
then
585 call post_data(cs%id_cp_mode(fr,m), cs%cp(:,:,fr,m), cs%diag)
586 endif ;
enddo ;
enddo
593 subroutine sum_en(G, CS, En, label)
597 real,
dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle), &
599 character(len=*),
intent(in) :: label
602 real :: tmpForSumming
612 en_sum = en_sum + tmpforsumming
635 subroutine itidal_lowmode_loss(G, US, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos)
640 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
642 real,
dimension(G%isd:G%ied,G%jsd:G%jed,CS%nFreq,CS%nMode), &
645 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
646 intent(in) :: TKE_loss_fixed
648 real,
dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
650 real,
dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
651 intent(out) :: TKE_loss
653 real,
intent(in) :: dt
654 logical,
optional,
intent(in) :: full_halos
657 integer :: j,i,m,fr,a, is, ie, js, je
660 real :: TKE_sum_check
661 real :: frac_per_sector
667 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
669 q_itides = cs%q_itides
670 en_negl = 1e-30*us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**2
672 if (
present(full_halos))
then ;
if (full_halos)
then
673 is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
676 do j=js,je ;
do i=is,ie ;
do m=1,cs%nMode ;
do fr=1,cs%nFreq
681 en_tot = en_tot + en(i,j,a,fr,m)
685 tke_loss_tot = q_itides * tke_loss_fixed(i,j) * nb(i,j) * ub(i,j,fr,m)**2
689 if (en_tot > 0.0)
then
691 frac_per_sector = en(i,j,a,fr,m)/en_tot
692 tke_loss(i,j,a,fr,m) = frac_per_sector*tke_loss_tot
693 loss_rate = tke_loss(i,j,a,fr,m) / (en(i,j,a,fr,m) + en_negl)
694 en(i,j,a,fr,m) = en(i,j,a,fr,m) / (1.0 + dt*loss_rate)
698 tke_loss(i,j,:,fr,m) = 0.0
719 enddo ;
enddo ;
enddo ;
enddo
728 integer,
intent(in) :: i
729 integer,
intent(in) :: j
733 character(len=*),
intent(in) :: mechanism
734 real,
intent(out) :: tke_loss_sum
737 if (mechanism ==
'LeakDrag') tke_loss_sum = cs%tot_leak_loss(i,j)
738 if (mechanism ==
'QuadDrag') tke_loss_sum = cs%tot_quad_loss(i,j)
739 if (mechanism ==
'WaveDrag') tke_loss_sum = cs%tot_itidal_loss(i,j)
740 if (mechanism ==
'Froude') tke_loss_sum = cs%tot_Froude_loss(i,j)
745 subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang)
747 integer,
intent(in) :: NAngle
749 real,
dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
753 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
755 real,
intent(in) :: freq
756 real,
intent(in) :: dt
758 logical,
intent(in) :: use_PPMang
761 integer,
parameter :: stencil = 2
762 real,
dimension(SZI_(G),1-stencil:NAngle+stencil) :: &
764 real,
dimension(1-stencil:NAngle+stencil) :: &
766 real,
dimension(SZI_(G)) :: &
767 Dk_Dt_Kmag, Dl_Dt_Kmag
768 real,
dimension(SZI_(G),0:nAngle) :: &
770 real,
dimension(SZI_(G),SZJ_(G),1-stencil:NAngle+stencil) :: &
777 real :: Angle_size, dt_Angle_size, angle
778 real :: Ifreq, Kmag2, I_Kmag
780 integer :: is, ie, js, je, asd, aed, na
783 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na =
size(en,3)
784 asd = 1-stencil ; aed = nangle+stencil
787 cn_subro = 1e-100*us%m_s_to_L_T
788 angle_size = (8.0*atan(1.0)) / (real(nangle))
789 dt_angle_size = dt / angle_size
792 angle = (real(a) - 0.5) * angle_size
793 cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
800 do a=1,na ;
do i=is,ie
801 en2d(i,a) = en(i,j,a)
803 do a=asd,0 ;
do i=is,ie
804 en2d(i,a) = en2d(i,a+nangle)
805 en2d(i,nangle+stencil+a) = en2d(i,stencil+a)
810 f2 = 0.25* ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
811 (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
812 favg = 0.25*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
813 (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j)))
814 df_dx = 0.5*((g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)) - &
815 (g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1))) * g%IdxT(i,j)
816 dlncn_dx = 0.5*( g%IdxCu(i,j) * (cn(i+1,j) - cn(i,j)) / &
817 (0.5*(cn(i+1,j) + cn(i,j)) + cn_subro) + &
818 g%IdxCu(i-1,j) * (cn(i,j) - cn(i-1,j)) / &
819 (0.5*(cn(i,j) + cn(i-1,j)) + cn_subro) )
820 df_dy = 0.5*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) - &
821 (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1))) * g%IdyT(i,j)
822 dlncn_dy = 0.5*( g%IdyCv(i,j) * (cn(i,j+1) - cn(i,j)) / &
823 (0.5*(cn(i,j+1) + cn(i,j)) + cn_subro) + &
824 g%IdyCv(i,j-1) * (cn(i,j) - cn(i,j-1)) / &
825 (0.5*(cn(i,j) + cn(i,j-1)) + cn_subro) )
826 kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subro**2)
827 if (kmag2 > 0.0)
then
828 i_kmag = 1.0 / sqrt(kmag2)
829 dk_dt_kmag(i) = -ifreq * (favg*df_dx + (freq**2 - f2) * dlncn_dx) * i_kmag
830 dl_dt_kmag(i) = -ifreq * (favg*df_dy + (freq**2 - f2) * dlncn_dy) * i_kmag
838 do a=asd,aed ;
do i=is,ie
839 cfl_ang(i,j,a) = (cos_angle(a) * dl_dt_kmag(i) - sin_angle(a) * dk_dt_kmag(i)) * dt_angle_size
840 if (abs(cfl_ang(i,j,a)) > 1.0)
then
841 call mom_error(warning,
"refract: CFL exceeds 1.", .true.)
842 if (cfl_ang(i,j,a) > 0.0)
then ; cfl_ang(i,j,a) = 1.0 ;
else ; cfl_ang(i,j,a) = -1.0 ;
endif
847 if (.not.use_ppmang)
then
849 do a=0,na ;
do i=is,ie
850 if (cfl_ang(i,j,a) > 0.0)
then
851 flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a)
853 flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a+1)
864 do a=1,na ;
do i=is,ie
868 en(i,j,a) = en2d(i,a) + (flux_e(i,a-1) - flux_e(i,a))
876 integer,
intent(in) :: NAngle
878 real,
intent(in) :: dt
879 integer,
intent(in) :: halo_ang
880 real,
dimension(1-halo_ang:NAngle+halo_ang), &
883 real,
dimension(1-halo_ang:NAngle+halo_ang), &
884 intent(in) :: CFL_ang
885 real,
dimension(0:NAngle),
intent(out) :: Flux_En
894 real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6
897 angle_size = (8.0*atan(1.0)) / (real(nangle))
898 i_angle_size = 1 / angle_size
902 u_ang = cfl_ang(a)*angle_size*i_dt
903 if (u_ang >= 0.0)
then
905 ep = en2d(a+1)*i_angle_size
906 ec = en2d(a) *i_angle_size
907 em = en2d(a-1)*i_angle_size
908 al = ( 5.*ec + ( 2.*em - ep ) )/6.
909 al = max( min(ec,em), al) ; al = min( max(ec,em), al)
910 ar = ( 5.*ec + ( 2.*ep - em ) )/6.
911 ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar)
912 da = ar - al ; ma = 0.5*( ar + al )
913 if ((ep-ec)*(ec-em) <= 0.)
then
915 elseif ( da*(ec-ma) > (da*da)/6. )
then
917 elseif ( da*(ec-ma) < - (da*da)/6. )
then
920 a6 = 6.*ec - 3. * (ar + al)
922 flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
925 flux_en(a) = dt * flux
929 ep = en2d(a+2)*i_angle_size
930 ec = en2d(a+1)*i_angle_size
931 em = en2d(a) *i_angle_size
932 al = ( 5.*ec + ( 2.*em - ep ) )/6.
933 al = max( min(ec,em), al) ; al = min( max(ec,em), al)
934 ar = ( 5.*ec + ( 2.*ep - em ) )/6.
935 ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar)
936 da = ar - al ; ma = 0.5*( ar + al )
937 if ((ep-ec)*(ec-em) <= 0.)
then
939 elseif ( da*(ec-ma) > (da*da)/6. )
then
941 elseif ( da*(ec-ma) < - (da*da)/6. )
then
944 a6 = 6.*ec - 3. * (ar + al)
946 flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
949 flux_en(a) = dt * flux
956 subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle)
958 integer,
intent(in) :: NAngle
960 real,
dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
964 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
966 real,
intent(in) :: freq
967 real,
intent(in) :: dt
972 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
974 integer,
parameter :: stencil = 2
975 real,
dimension(SZIB_(G),SZJ_(G)) :: &
977 real,
dimension(SZI_(G),SZJB_(G)) :: &
979 real,
dimension(0:NAngle) :: &
981 real,
dimension(NAngle) :: &
982 Cgx_av, Cgy_av, dCgx, dCgy
984 real :: Angle_size, I_Angle_size, angle
988 integer :: is, ie, js, je, asd, aed, na
989 integer :: ish, ieh, jsh, jeh
992 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na =
size(en,3)
993 asd = 1-stencil ; aed = nangle+stencil
1007 jsh = js-1 ; jeh = je+1 ; ish = is-1 ; ieh = ie+1
1009 angle_size = (8.0*atan(1.0)) / real(nangle)
1010 i_angle_size = 1.0 / angle_size
1012 if (cs%corner_adv)
then
1018 do j=jsh-1,jeh ;
do i=ish-1,ieh
1019 f2 = g%CoriolisBu(i,j)**2
1020 speed(i,j) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * &
1021 sqrt(max(freq2 - f2, 0.0)) * ifreq
1025 lb%jsh = js ; lb%jeh = je ; lb%ish = is ; lb%ieh = ie
1033 angle = (real(a) - 0.5) * angle_size
1034 cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
1038 cgx_av(a) = (sin_angle(a) - sin_angle(a-1)) * i_angle_size
1039 cgy_av(a) = -(cos_angle(a) - cos_angle(a-1)) * i_angle_size
1040 dcgx(a) = sqrt(0.5 + 0.5*(sin_angle(a)*cos_angle(a) - &
1041 sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
1043 dcgy(a) = sqrt(0.5 - 0.5*(sin_angle(a)*cos_angle(a) - &
1044 sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
1048 do j=jsh,jeh ;
do i=ish-1,ieh
1049 f2 = 0.5 * (g%CoriolisBu(i,j)**2 + g%CoriolisBu(i,j-1)**2)
1050 speed_x(i,j) = 0.5*(cn(i,j) + cn(i+1,j)) * g%mask2dCu(i,j) * &
1051 sqrt(max(freq2 - f2, 0.0)) * ifreq
1053 do j=jsh-1,jeh ;
do i=ish,ieh
1054 f2 = 0.5 * (g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2)
1055 speed_y(i,j) = 0.5*(cn(i,j) + cn(i,j+1)) * g%mask2dCv(i,j) * &
1056 sqrt(max(freq2 - f2, 0.0)) * ifreq
1060 lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
1061 call propagate_x(en(:,:,:), speed_x, cgx_av(:), dcgx(:), dt, g, us, cs%nAngle, cs, lb)
1071 lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
1072 call propagate_y(en(:,:,:), speed_y, cgy_av(:), dcgy(:), dt, g, us, cs%nAngle, cs, lb)
1085 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
1088 real,
dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), &
1091 integer,
intent(in) :: energized_wedge
1092 integer,
intent(in) :: NAngle
1094 real,
intent(in) :: dt
1099 integer :: i, j, k, ish, ieh, jsh, jeh, m
1100 real :: TwoPi, Angle_size
1101 real :: energized_angle
1106 real :: I_Nsubwedges
1107 real :: cos_thetaDT, sin_thetaDT
1108 real :: xNE,xNW,xSW,xSE,yNE,yNW,ySW,ySE
1109 real :: CFL_xNE,CFL_xNW,CFL_xSW,CFL_xSE,CFL_yNE,CFL_yNW,CFL_ySW,CFL_ySE,CFL_max
1110 real :: xN,xS,xE,xW,yN,yS,yE,yW
1113 real :: slopeN,slopeW,slopeS,slopeE, bN,bW,bS,bE
1114 real :: aNE,aN,aNW,aW,aSW,aS,aSE,aE,aC
1117 real,
dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: x,y
1118 real,
dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: Idx,Idy
1119 real,
dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: dx,dy
1120 real,
dimension(2) :: E_new
1123 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1124 twopi = (8.0*atan(1.0))
1125 nsubrays = real(
size(e_new))
1126 i_nsubwedges = 1./(nsubrays - 1)
1128 angle_size = twopi / real(nangle)
1129 energized_angle = angle_size * real(energized_wedge - 1)
1132 do j=jsh-1,jeh ;
do i=ish-1,ieh
1135 x(i,j) = g%US%m_to_L*g%geoLonBu(i,j)
1136 y(i,j) = g%US%m_to_L*g%geoLatBu(i,j)
1137 idx(i,j) = g%IdxBu(i,j) ; dx(i,j) = g%dxBu(i,j)
1138 idy(i,j) = g%IdyBu(i,j) ; dy(i,j) = g%dyBu(i,j)
1141 do j=jsh,jeh ;
do i=ish,ieh
1142 do m=1,int(nsubrays)
1143 theta = energized_angle - 0.5*angle_size + real(m - 1)*angle_size*i_nsubwedges
1144 if (theta < 0.0)
then
1145 theta = theta + twopi
1146 elseif (theta > twopi)
then
1147 theta = theta - twopi
1149 cos_thetadt = cos(theta)*dt
1150 sin_thetadt = sin(theta)*dt
1153 xg = x(i,j); yg = y(i,j)
1154 xne = xg - speed(i,j)*cos_thetadt
1155 yne = yg - speed(i,j)*sin_thetadt
1156 cfl_xne = (xg-xne)*idx(i,j)
1157 cfl_yne = (yg-yne)*idy(i,j)
1159 xg = x(i-1,j); yg = y(i-1,j)
1160 xnw = xg - speed(i-1,j)*cos_thetadt
1161 ynw = yg - speed(i-1,j)*sin_thetadt
1162 cfl_xnw = (xg-xnw)*idx(i-1,j)
1163 cfl_ynw = (yg-ynw)*idy(i-1,j)
1165 xg = x(i-1,j-1); yg = y(i-1,j-1)
1166 xsw = xg - speed(i-1,j-1)*cos_thetadt
1167 ysw = yg - speed(i-1,j-1)*sin_thetadt
1168 cfl_xsw = (xg-xsw)*idx(i-1,j-1)
1169 cfl_ysw = (yg-ysw)*idy(i-1,j-1)
1171 xg = x(i,j-1); yg = y(i,j-1)
1172 xse = xg - speed(i,j-1)*cos_thetadt
1173 yse = yg - speed(i,j-1)*sin_thetadt
1174 cfl_xse = (xg-xse)*idx(i,j-1)
1175 cfl_yse = (yg-yse)*idy(i,j-1)
1177 cfl_max = max(abs(cfl_xne),abs(cfl_xnw),abs(cfl_xsw), &
1178 abs(cfl_xse),abs(cfl_yne),abs(cfl_ynw), &
1179 abs(cfl_ysw),abs(cfl_yse))
1180 if (cfl_max > 1.0)
then
1181 call mom_error(warning,
"propagate_corner_spread: CFL exceeds 1.", .true.)
1185 if (0.0 <= theta .and. theta < 0.25*twopi)
then
1188 elseif (0.25*twopi <= theta .and. theta < 0.5*twopi)
then
1191 elseif (0.5*twopi <= theta .and. theta < 0.75*twopi)
then
1194 elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi)
then
1202 slopen = (yne - ynw)/(xne - xnw)
1203 bn = -slopen*xne + yne
1206 if (xnw == xsw)
then
1209 slopew = (ynw - ysw)/(xnw - xsw)
1210 bw = -slopew*xnw + ynw
1211 xw = (yw - bw)/slopew
1214 slopes = (ysw - yse)/(xsw - xse)
1215 bs = -slopes*xsw + ysw
1218 if (xne == xse)
then
1221 slopee = (yse - yne)/(xse - xne)
1222 be = -slopee*xse + yse
1223 xe = (ye - be)/slopee
1227 ane = 0.0; an = 0.0; anw = 0.0;
1228 aw = 0.0; asw = 0.0; as = 0.0;
1229 ase = 0.0; ae = 0.0; ac = 0.0;
1230 if (0.0 <= theta .and. theta < 0.25*twopi)
then
1231 xcrn = x(i-1,j-1); ycrn = y(i-1,j-1)
1233 a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1234 a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1235 a3 = (yw - ynw)*(0.5*(xw + xnw))
1236 a4 = (ynw - yn)*(0.5*(xnw + xn))
1237 aw = a1 + a2 + a3 + a4
1239 a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1240 a2 = (ys - ysw)*(0.5*(xs + xsw))
1241 a3 = (ysw - yw)*(0.5*(xsw + xw))
1242 a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1243 asw = a1 + a2 + a3 + a4
1245 a1 = (ye - yse)*(0.5*(xe + xse))
1246 a2 = (yse - ys)*(0.5*(xse + xs))
1247 a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1248 a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1249 as = a1 + a2 + a3 + a4
1251 a1 = (yne - ye)*(0.5*(xne + xe))
1252 a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1253 a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1254 a4 = (yn - yne)*(0.5*(xn + xne))
1255 ac = a1 + a2 + a3 + a4
1256 elseif (0.25*twopi <= theta .and. theta < 0.5*twopi)
then
1257 xcrn = x(i,j-1); ycrn = y(i,j-1)
1259 a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1260 a2 = (ys - ysw)*(0.5*(xs + xsw))
1261 a3 = (ysw - yw)*(0.5*(xsw + xw))
1262 a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1263 as = a1 + a2 + a3 + a4
1265 a1 = (ye - yse)*(0.5*(xe + xse))
1266 a2 = (yse - ys)*(0.5*(xse + xs))
1267 a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1268 a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1269 ase = a1 + a2 + a3 + a4
1271 a1 = (yne - ye)*(0.5*(xne + xe))
1272 a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1273 a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1274 a4 = (yn - yne)*(0.5*(xn + xne))
1275 ae = a1 + a2 + a3 + a4
1277 a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1278 a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1279 a3 = (yw - ynw)*(0.5*(xw + xnw))
1280 a4 = (ynw - yn)*(0.5*(xnw + xn))
1281 ac = a1 + a2 + a3 + a4
1282 elseif (0.5*twopi <= theta .and. theta < 0.75*twopi)
then
1283 xcrn = x(i,j); ycrn = y(i,j)
1285 a1 = (ye - yse)*(0.5*(xe + xse))
1286 a2 = (yse - ys)*(0.5*(xse + xs))
1287 a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1288 a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1289 ae = a1 + a2 + a3 + a4
1291 a1 = (yne - ye)*(0.5*(xne + xe))
1292 a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1293 a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1294 a4 = (yn - yne)*(0.5*(xn + xne))
1295 ane = a1 + a2 + a3 + a4
1297 a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1298 a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1299 a3 = (yw - ynw)*(0.5*(xw + xnw))
1300 a4 = (ynw - yn)*(0.5*(xnw + xn))
1301 an = a1 + a2 + a3 + a4
1303 a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1304 a2 = (ys - ysw)*(0.5*(xs + xsw))
1305 a3 = (ysw - yw)*(0.5*(xsw + xw))
1306 a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1307 ac = a1 + a2 + a3 + a4
1308 elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi)
then
1309 xcrn = x(i-1,j); ycrn = y(i-1,j)
1311 a1 = (yne - ye)*(0.5*(xne + xe))
1312 a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1313 a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1314 a4 = (yn - yne)*(0.5*(xn + xne))
1315 an = a1 + a2 + a3 + a4
1317 a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1318 a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1319 a3 = (yw - ynw)*(0.5*(xw + xnw))
1320 a4 = (ynw - yn)*(0.5*(xnw + xn))
1321 anw = a1 + a2 + a3 + a4
1323 a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1324 a2 = (ys - ysw)*(0.5*(xs + xsw))
1325 a3 = (ysw - yw)*(0.5*(xsw + xw))
1326 a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1327 aw = a1 + a2 + a3 + a4
1329 a1 = (ye - yse)*(0.5*(xe + xse))
1330 a2 = (yse - ys)*(0.5*(xse + xs))
1331 a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1332 a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1333 ac = a1 + a2 + a3 + a4
1337 a_total = ane + an + anw + aw + asw + as + ase + ae + ac
1338 e_new(m) = (ane*en(i+1,j+1) + an*en(i,j+1) + anw*en(i-1,j+1) + &
1339 aw*en(i-1,j) + asw*en(i-1,j-1) + as*en(i,j-1) + &
1340 ase*en(i+1,j-1) + ae*en(i+1,j) + ac*en(i,j)) / (dx(i,j)*dy(i,j))
1343 en(i,j) = sum(e_new)/nsubrays
1348 subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB)
1350 integer,
intent(in) :: NAngle
1352 real,
dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
1355 real,
dimension(G%IsdB:G%IedB,G%jsd:G%jed), &
1356 intent(in) :: speed_x
1358 real,
dimension(Nangle),
intent(in) :: Cgx_av
1359 real,
dimension(Nangle),
intent(in) :: dCgx
1361 real,
intent(in) :: dt
1367 real,
dimension(SZI_(G),SZJ_(G)) :: &
1369 real,
dimension(SZIB_(G),SZJ_(G)) :: &
1371 real,
dimension(SZIB_(G)) :: &
1372 cg_p, cg_m, flux1, flux2
1374 real,
dimension(SZI_(G),SZJB_(G),Nangle) :: &
1376 integer :: i, j, k, ish, ieh, jsh, jeh, a
1378 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1381 if (cs%upwind_1st)
then
1382 do j=jsh,jeh ;
do i=ish-1,ieh+1
1383 enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
1392 cg_p(i) = speed_x(i,j) * (cgx_av(a))
1394 call zonal_flux_en(cg_p, en(:,j,a), enl(:,j), enr(:,j), flux1, &
1395 dt, g, us, j, ish, ieh, cs%vol_CFL)
1396 do i=ish-1,ieh ; flux_x(i,j) = flux1(i);
enddo
1399 do j=jsh,jeh ;
do i=ish,ieh
1400 fdt_m(i,j,a) = dt*flux_x(i-1,j)
1401 fdt_p(i,j,a) = -dt*flux_x(i,j)
1408 call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
1409 call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
1410 call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
1411 call teleport(fdt_p(:,:,:), nangle, cs, g, lb)
1414 do j=jsh,jeh ;
do i=ish,ieh
1417 en(i,j,:) = en(i,j,:) + g%IareaT(i,j)*(fdt_m(i,j,:) + fdt_p(i,j,:))
1423 subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB)
1425 integer,
intent(in) :: NAngle
1427 real,
dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
1430 real,
dimension(G%isd:G%ied,G%JsdB:G%JedB), &
1431 intent(in) :: speed_y
1433 real,
dimension(Nangle),
intent(in) :: Cgy_av
1434 real,
dimension(Nangle),
intent(in) :: dCgy
1436 real,
intent(in) :: dt
1442 real,
dimension(SZI_(G),SZJ_(G)) :: &
1444 real,
dimension(SZI_(G),SZJB_(G)) :: &
1446 real,
dimension(SZI_(G)) :: &
1447 cg_p, cg_m, flux1, flux2
1449 real,
dimension(SZI_(G),SZJB_(G),Nangle) :: &
1451 character(len=160) :: mesg
1452 integer :: i, j, k, ish, ieh, jsh, jeh, a
1454 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1457 if (cs%upwind_1st)
then
1458 do j=jsh-1,jeh+1 ;
do i=ish,ieh
1459 enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
1468 cg_p(i) = speed_y(i,j) * (cgy_av(a))
1470 call merid_flux_en(cg_p, en(:,:,a), enl(:,:), enr(:,:), flux1, &
1471 dt, g, us, j, ish, ieh, cs%vol_CFL)
1472 do i=ish,ieh ; flux_y(i,j) = flux1(i);
enddo
1475 do j=jsh,jeh ;
do i=ish,ieh
1476 fdt_m(i,j,a) = dt*flux_y(i,j-1)
1477 fdt_p(i,j,a) = -dt*flux_y(i,j)
1490 call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
1491 call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
1492 call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
1493 call teleport(fdt_p(:,:,:), nangle, cs, g, lb)
1496 do a=1,nangle ;
do j=jsh,jeh ;
do i=ish,ieh
1499 en(i,j,a) = en(i,j,a) + g%IareaT(i,j)*(fdt_m(i,j,a) + fdt_p(i,j,a))
1500 enddo ;
enddo ;
enddo
1505 subroutine zonal_flux_en(u, h, hL, hR, uh, dt, G, US, j, ish, ieh, vol_CFL)
1507 real,
dimension(SZIB_(G)),
intent(in) :: u
1508 real,
dimension(SZI_(G)),
intent(in) :: h
1510 real,
dimension(SZI_(G)),
intent(in) :: hL
1512 real,
dimension(SZI_(G)),
intent(in) :: hR
1514 real,
dimension(SZIB_(G)),
intent(inout) :: uh
1515 real,
intent(in) :: dt
1517 integer,
intent(in) :: j
1518 integer,
intent(in) :: ish
1519 integer,
intent(in) :: ieh
1520 logical,
intent(in) :: vol_CFL
1530 if (u(i) > 0.0)
then
1531 if (vol_cfl)
then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
1532 else ; cfl = u(i) * dt * g%IdxT(i,j) ;
endif
1533 curv_3 = (hl(i) + hr(i)) - 2.0*h(i)
1534 uh(i) = g%dy_Cu(i,j) * u(i) * &
1535 (hr(i) + cfl * (0.5*(hl(i) - hr(i)) + curv_3*(cfl - 1.5)))
1536 elseif (u(i) < 0.0)
then
1537 if (vol_cfl)
then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
1538 else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ;
endif
1539 curv_3 = (hl(i+1) + hr(i+1)) - 2.0*h(i+1)
1540 uh(i) = g%dy_Cu(i,j) * u(i) * &
1541 (hl(i+1) + cfl * (0.5*(hr(i+1)-hl(i+1)) + curv_3*(cfl - 1.5)))
1549 subroutine merid_flux_en(v, h, hL, hR, vh, dt, G, US, J, ish, ieh, vol_CFL)
1551 real,
dimension(SZI_(G)),
intent(in) :: v
1552 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h
1554 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: hL
1556 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: hR
1558 real,
dimension(SZI_(G)),
intent(inout) :: vh
1559 real,
intent(in) :: dt
1561 integer,
intent(in) :: J
1562 integer,
intent(in) :: ish
1563 integer,
intent(in) :: ieh
1564 logical,
intent(in) :: vol_CFL
1574 if (v(i) > 0.0)
then
1575 if (vol_cfl)
then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1576 else ; cfl = v(i) * dt * g%IdyT(i,j) ;
endif
1577 curv_3 = hl(i,j) + hr(i,j) - 2.0*h(i,j)
1578 vh(i) = g%dx_Cv(i,j) * v(i) * ( hr(i,j) + cfl * &
1579 (0.5*(hl(i,j) - hr(i,j)) + curv_3*(cfl - 1.5)) )
1580 elseif (v(i) < 0.0)
then
1581 if (vol_cfl)
then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1582 else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ;
endif
1583 curv_3 = hl(i,j+1) + hr(i,j+1) - 2.0*h(i,j+1)
1584 vh(i) = g%dx_Cv(i,j) * v(i) * ( hl(i,j+1) + cfl * &
1585 (0.5*(hr(i,j+1)-hl(i,j+1)) + curv_3*(cfl - 1.5)) )
1593 subroutine reflect(En, NAngle, CS, G, LB)
1595 integer,
intent(in) :: NAngle
1597 real,
dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1605 real,
dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
1607 real,
dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
1610 logical,
dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
1616 real,
dimension(1:NAngle) :: angle_i
1618 real,
dimension(1:Nangle) :: En_reflected
1619 integer :: i, j, a, a_r, na
1622 integer :: isc, iec, jsc, jec
1624 integer :: ish, ieh, jsh, jeh
1626 integer :: id_g, jd_g
1629 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1630 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1632 twopi = 8.0*atan(1.0)
1633 angle_size = twopi / (real(nangle))
1638 angle_i(a) = angle_size * real(a - 1)
1641 angle_c = cs%refl_angle
1642 part_refl = cs%refl_pref
1644 en_reflected(:) = 0.0
1653 if (angle_c(i,j) /= cs%nullangle)
then
1655 if (en(i,j,a) > 0.0)
then
1657 if (sin(angle_i(a) - angle_c(i,j)) >= 0.0)
then
1658 angle_wall = angle_c(i,j)
1660 elseif (ridge(i,j))
then
1661 angle_wall = angle_c(i,j) + 0.5*twopi
1662 if (angle_wall > twopi)
then
1663 angle_wall = angle_wall - twopi*floor(abs(angle_wall)/twopi)
1664 elseif (angle_wall < 0.0)
then
1665 angle_wall = angle_wall + twopi*ceiling(abs(angle_wall)/twopi)
1669 angle_wall = angle_c(i,j)
1672 if (sin(angle_i(a) - angle_wall) >= 0.0)
then
1673 angle_r = 2.0*angle_wall - angle_i(a)
1674 if (angle_r > twopi)
then
1675 angle_r = angle_r - twopi*floor(abs(angle_r)/twopi)
1676 elseif (angle_r < 0.0)
then
1677 angle_r = angle_r + twopi*ceiling(abs(angle_r)/twopi)
1679 a_r = nint(angle_r/angle_size) + 1
1680 do while (a_r > nangle) ; a_r = a_r - nangle ;
enddo
1682 en_reflected(a_r) = part_refl(i,j)*en(i,j,a)
1683 en(i,j,a) = (1.0-part_refl(i,j))*en(i,j,a)
1688 en(i,j,:) = en(i,j,:) + en_reflected(:)
1689 en_reflected(:) = 0.0
1707 subroutine teleport(En, NAngle, CS, G, LB)
1709 integer,
intent(in) :: NAngle
1711 real,
dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1719 real,
dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
1721 real,
dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
1724 logical,
dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell
1726 logical,
dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
1730 real,
dimension(1:NAngle) :: angle_i
1731 real,
dimension(1:NAngle) :: cos_angle, sin_angle
1733 character(len=160) :: mesg
1739 integer :: ish, ieh, jsh, jeh
1741 integer :: id_g, jd_g
1743 real :: cos_normal, sin_normal, angle_wall
1748 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1750 twopi = 8.0*atan(1.0)
1751 angle_size = twopi / (real(nangle))
1756 angle_i(a) = angle_size * real(a - 1)
1757 cos_angle(a) = cos(angle_i(a)) ; sin_angle(a) = sin(angle_i(a))
1760 angle_c = cs%refl_angle
1761 part_refl = cs%refl_pref
1762 pref_cell = cs%refl_pref_logical
1767 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
1768 if (pref_cell(i,j))
then
1770 if (en(i,j,a) > 0)
then
1772 if (sin(angle_i(a) - angle_c(i,j)) >= 0.0)
then
1773 angle_wall = angle_c(i,j)
1775 elseif (ridge(i,j))
then
1776 angle_wall = angle_c(i,j) + 0.5*twopi
1779 angle_wall = angle_c(i,j)
1782 if (sin(angle_i(a) - angle_wall) >= 0.0)
then
1784 cos_normal = cos(angle_wall + 0.25*twopi)
1785 sin_normal = sin(angle_wall + 0.25*twopi)
1787 ios = int(sign(1.,cos_normal))
1789 jos = int(sign(1.,sin_normal))
1791 if (.not. pref_cell(i+ios,j+jos))
then
1792 en(i,j,a) = en(i,j,a) - en_tele
1793 en(i+ios,j+jos,a) = en(i+ios,j+jos,a) + en_tele
1795 write(mesg,*)
'idg=',id_g,
'jd_g=',jd_g,
'a=',a
1796 call mom_error(fatal,
"teleport: no receptive ocean cell at "//trim(mesg), .true.)
1811 real,
dimension(:,:,:,:,:),
intent(inout) :: En
1814 real,
dimension(SZI_(G),SZJ_(G),2), &
1818 integer,
intent(in) :: NAngle
1821 real,
dimension(G%isd:G%ied,NAngle) :: En2d
1822 integer,
dimension(G%isd:G%ied) :: a_shift
1823 integer :: i_first, i_last, a_new
1824 integer :: a, i, j, isd, ied, jsd, jed, m, fr
1825 character(len=160) :: mesg
1826 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1829 i_first = ied+1 ; i_last = isd-1
1832 if (test(i,j,1) /= 1.0)
then
1833 if (i<i_first) i_first = i
1834 if (i>i_last) i_last = i
1836 if (test(i,j,1) == -1.0)
then ; a_shift(i) = nangle/2
1837 elseif (test(i,j,2) == 1.0)
then ; a_shift(i) = -nangle/4
1838 elseif (test(i,j,2) == -1.0)
then ; a_shift(i) = nangle/4
1840 write(mesg,
'("Unrecognized rotation test vector ",2ES9.2," at ",F7.2," E, ",&
1841 &F7.2," N; i,j=",2i4)') &
1842 test(i,j,1), test(i,j,2), g%GeoLonT(i,j), g%GeoLatT(i,j), i, j
1848 if (i_first <= i_last)
then
1850 do m=1,
size(en,5) ;
do fr=1,
size(en,4)
1851 do a=1,nangle ;
do i=i_first,i_last ;
if (a_shift(i) /= 0)
then
1852 a_new = a + a_shift(i)
1853 if (a_new < 1) a_new = a_new + nangle
1854 if (a_new > nangle) a_new = a_new - nangle
1855 en2d(i,a_new) = en(i,j,a,fr,m)
1856 endif ;
enddo ;
enddo
1857 do a=1,nangle ;
do i=i_first,i_last ;
if (a_shift(i) /= 0)
then
1858 en(i,j,a,fr,m) = en2d(i,a)
1859 endif ;
enddo ;
enddo
1868 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
1869 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_l
1870 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_r
1872 logical,
optional,
intent(in) :: simple_2nd
1876 real,
dimension(SZI_(G),SZJ_(G)) :: slp
1877 real,
parameter :: oneSixth = 1./6.
1878 real :: h_ip1, h_im1
1880 logical :: use_CW84, use_2nd
1881 character(len=256) :: mesg
1882 integer :: i, j, isl, iel, jsl, jel, stencil
1884 use_2nd = .false. ;
if (
present(simple_2nd)) use_2nd = simple_2nd
1885 isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1888 stencil = 2 ;
if (use_2nd) stencil = 1
1890 if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied))
then
1891 write(mesg,
'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
1892 & "x-halo that needs to be increased by ",i2,".")') &
1893 stencil + max(g%isd-isl,iel-g%ied)
1896 if ((jsl < g%jsd) .or. (jel > g%jed))
then
1897 write(mesg,
'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
1898 & "y-halo that needs to be increased by ",i2,".")') &
1899 max(g%jsd-jsl,jel-g%jed)
1904 do j=jsl,jel ;
do i=isl,iel
1905 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1906 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1907 h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1908 h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1911 do j=jsl,jel ;
do i=isl-1,iel+1
1912 if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0)
then
1916 slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1918 dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1919 dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1920 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1925 do j=jsl,jel ;
do i=isl,iel
1930 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1931 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1933 h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1934 h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1938 call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
1944 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
1945 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_l
1946 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_r
1948 logical,
optional,
intent(in) :: simple_2nd
1952 real,
dimension(SZI_(G),SZJ_(G)) :: slp
1953 real,
parameter :: oneSixth = 1./6.
1954 real :: h_jp1, h_jm1
1957 character(len=256) :: mesg
1958 integer :: i, j, isl, iel, jsl, jel, stencil
1960 use_2nd = .false. ;
if (
present(simple_2nd)) use_2nd = simple_2nd
1961 isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
1964 stencil = 2 ;
if (use_2nd) stencil = 1
1966 if ((isl < g%isd) .or. (iel > g%ied))
then
1967 write(mesg,
'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
1968 & "x-halo that needs to be increased by ",i2,".")') &
1969 max(g%isd-isl,iel-g%ied)
1972 if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed))
then
1973 write(mesg,
'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
1974 & "y-halo that needs to be increased by ",i2,".")') &
1975 stencil + max(g%jsd-jsl,jel-g%jed)
1980 do j=jsl,jel ;
do i=isl,iel
1981 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
1982 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
1983 h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
1984 h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
1987 do j=jsl-1,jel+1 ;
do i=isl,iel
1988 if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0)
then
1992 slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
1994 dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
1995 dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
1996 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
2001 do j=jsl,jel ;
do i=isl,iel
2004 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2005 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2007 h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2008 h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2012 call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
2019 subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
2021 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
2022 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: h_L
2023 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: h_R
2024 real,
intent(in) :: h_min
2026 integer,
intent(in) :: iis
2027 integer,
intent(in) :: iie
2028 integer,
intent(in) :: jis
2029 integer,
intent(in) :: jie
2031 real :: curv, dh, scale
2032 character(len=256) :: mesg
2035 do j=jis,jie ;
do i=iis,iie
2038 curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2039 if (curv > 0.0)
then
2040 dh = h_r(i,j) - h_l(i,j)
2041 if (abs(dh) < curv)
then
2042 if (h_in(i,j) <= h_min)
then
2043 h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2044 elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2))
then
2047 scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2048 h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2049 h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2100 type(time_type),
target,
intent(in) :: time
2106 type(
diag_ctrl),
target,
intent(in) :: diag
2112 real,
allocatable :: angles(:)
2113 real,
allocatable,
dimension(:,:) :: h2
2114 real :: kappa_itides, kappa_h2_factor
2117 real,
allocatable :: ridge_temp(:,:)
2120 logical :: use_int_tides, use_temperature
2122 integer :: num_angle, num_freq, num_mode, m, fr
2123 integer :: isd, ied, jsd, jed, a, id_ang, i, j
2126 #include "version_variable.h"
2127 character(len=40) :: mdl =
"MOM_internal_tides"
2128 character(len=16),
dimension(8) :: freq_name
2129 character(len=40) :: var_name
2130 character(len=160) :: var_descript
2131 character(len=200) :: filename
2132 character(len=200) :: refl_angle_file, land_mask_file
2133 character(len=200) :: refl_pref_file, refl_dbl_file
2134 character(len=200) :: dy_cu_file, dx_cv_file
2135 character(len=200) :: h2_file
2137 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2139 if (
associated(cs))
then
2140 call mom_error(warning,
"internal_tides_init called "//&
2141 "with an associated control structure.")
2147 use_int_tides = .false.
2148 call read_param(param_file,
"INTERNAL_TIDES", use_int_tides)
2149 cs%do_int_tides = use_int_tides
2150 if (.not.use_int_tides)
return
2152 use_temperature = .true.
2153 call read_param(param_file,
"ENABLE_THERMODYNAMICS", use_temperature)
2154 if (.not.use_temperature)
call mom_error(fatal, &
2155 "register_int_tide_restarts: internal_tides only works with "//&
2156 "ENABLE_THERMODYNAMICS defined.")
2159 num_freq = 1 ; num_angle = 24 ; num_mode = 1
2160 call read_param(param_file,
"INTERNAL_TIDE_FREQS", num_freq)
2161 call read_param(param_file,
"INTERNAL_TIDE_ANGLES", num_angle)
2162 call read_param(param_file,
"INTERNAL_TIDE_MODES", num_mode)
2163 if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0)))
return
2164 cs%nFreq = num_freq ; cs%nAngle = num_angle ; cs%nMode = num_mode
2167 allocate(cs%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode))
2168 cs%En(:,:,:,:,:) = 0.0
2171 allocate(cs%cp(isd:ied, jsd:jed, num_freq, num_mode))
2172 cs%cp(:,:,:,:) = 0.0
2175 allocate(cs%frequency(num_freq))
2176 call get_param(param_file, mdl,
"FIRST_MODE_PERIOD", period_1, units=
"s", scale=us%s_to_T)
2178 cs%frequency(fr) = (8.0*atan(1.0) * (real(fr)) / period_1)
2185 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, default=
".")
2186 cs%inputdir = slasher(cs%inputdir)
2190 call get_param(param_file, mdl,
"INTERNAL_TIDE_FREQS", num_freq, &
2191 "The number of distinct internal tide frequency bands "//&
2192 "that will be calculated.", default=1)
2193 call get_param(param_file, mdl,
"INTERNAL_TIDE_MODES", num_mode, &
2194 "The number of distinct internal tide modes "//&
2195 "that will be calculated.", default=1)
2196 call get_param(param_file, mdl,
"INTERNAL_TIDE_ANGLES", num_angle, &
2197 "The number of angular resolution bands for the internal "//&
2198 "tide calculations.", default=24)
2200 if (use_int_tides)
then
2201 if ((num_freq <= 0) .and. (num_mode <= 0) .and. (num_angle <= 0))
then
2202 call mom_error(warning,
"Internal tides were enabled, but the number "//&
2203 "of requested frequencies, modes and angles were not all positive.")
2207 if ((num_freq > 0) .and. (num_mode > 0) .and. (num_angle > 0))
then
2208 call mom_error(warning,
"Internal tides were not enabled, even though "//&
2209 "the number of requested frequencies, modes and angles were all "//&
2215 if (cs%NFreq /= num_freq)
call mom_error(fatal,
"Internal_tides_init: "//&
2216 "Inconsistent number of frequencies.")
2217 if (cs%NAngle /= num_angle)
call mom_error(fatal,
"Internal_tides_init: "//&
2218 "Inconsistent number of angles.")
2219 if (cs%NMode /= num_mode)
call mom_error(fatal,
"Internal_tides_init: "//&
2220 "Inconsistent number of modes.")
2221 if (4*(num_angle/4) /= num_angle)
call mom_error(fatal, &
2222 "Internal_tides_init: INTERNAL_TIDE_ANGLES must be a multiple of 4.")
2226 call get_param(param_file, mdl,
"INTERNAL_TIDE_DECAY_RATE", cs%decay_rate, &
2227 "The rate at which internal tide energy is lost to the "//&
2228 "interior ocean internal wave field.", &
2229 units=
"s-1", default=0.0, scale=us%T_to_s)
2230 call get_param(param_file, mdl,
"INTERNAL_TIDE_VOLUME_BASED_CFL", cs%vol_CFL, &
2231 "If true, use the ratio of the open face lengths to the "//&
2232 "tracer cell areas when estimating CFL numbers in the "//&
2233 "internal tide code.", default=.false.)
2234 call get_param(param_file, mdl,
"INTERNAL_TIDE_CORNER_ADVECT", cs%corner_adv, &
2235 "If true, internal tide ray-tracing advection uses a "//&
2236 "corner-advection scheme rather than PPM.", default=.false.)
2237 call get_param(param_file, mdl,
"INTERNAL_TIDE_SIMPLE_2ND_PPM", cs%simple_2nd, &
2238 "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2239 "(arithmetic mean) interpolation of the edge values. "//&
2240 "This may give better PV conservation properties. While "//&
2241 "it formally reduces the accuracy of the continuity "//&
2242 "solver itself in the strongly advective limit, it does "//&
2243 "not reduce the overall order of accuracy of the dynamic "//&
2244 "core.", default=.false.)
2245 call get_param(param_file, mdl,
"INTERNAL_TIDE_UPWIND_1ST", cs%upwind_1st, &
2246 "If true, the internal tide ray-tracing advection uses "//&
2247 "1st-order upwind advection. This scheme is highly "//&
2248 "continuity solver. This scheme is highly "//&
2249 "diffusive but may be useful for debugging.", default=.false.)
2250 call get_param(param_file, mdl,
"INTERNAL_TIDE_BACKGROUND_DRAG", &
2251 cs%apply_background_drag,
"If true, the internal tide "//&
2252 "ray-tracing advection uses a background drag term as a sink.",&
2254 call get_param(param_file, mdl,
"INTERNAL_TIDE_QUAD_DRAG", cs%apply_bottom_drag, &
2255 "If true, the internal tide ray-tracing advection uses "//&
2256 "a quadratic bottom drag term as a sink.", default=.false.)
2257 call get_param(param_file, mdl,
"INTERNAL_TIDE_WAVE_DRAG", cs%apply_wave_drag, &
2258 "If true, apply scattering due to small-scale roughness as a sink.", &
2260 call get_param(param_file, mdl,
"INTERNAL_TIDE_FROUDE_DRAG", cs%apply_Froude_drag, &
2261 "If true, apply wave breaking as a sink.", &
2263 call get_param(param_file, mdl,
"CDRAG", cs%cdrag, &
2264 "CDRAG is the drag coefficient relating the magnitude of "//&
2265 "the velocity field to the bottom stress.", units=
"nondim", &
2267 call get_param(param_file, mdl,
"INTERNAL_TIDE_ENERGIZED_ANGLE", cs%energized_angle, &
2268 "If positive, only one angular band of the internal tides "//&
2269 "gets all of the energy. (This is for debugging.)", default=-1)
2270 call get_param(param_file, mdl,
"USE_PPM_ANGULAR", cs%use_PPMang, &
2271 "If true, use PPM for advection of energy in angular space.", &
2273 call get_param(param_file, mdl,
"GAMMA_ITIDES", cs%q_itides, &
2274 "The fraction of the internal tidal energy that is "//&
2275 "dissipated locally with INT_TIDE_DISSIPATION. "//&
2276 "THIS NAME COULD BE BETTER.", &
2277 units=
"nondim", default=0.3333)
2278 call get_param(param_file, mdl,
"KAPPA_ITIDES", kappa_itides, &
2279 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
2280 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
2281 units=
"m-1", default=8.e-4*atan(1.0), scale=us%L_to_m)
2282 call get_param(param_file, mdl,
"KAPPA_H2_FACTOR", kappa_h2_factor, &
2283 "A scaling factor for the roughness amplitude with n"//&
2284 "INT_TIDE_DISSIPATION.", units=
"nondim", default=1.0)
2287 allocate(h2(isd:ied,jsd:jed)) ; h2(:,:) = 0.0
2288 allocate(cs%TKE_itidal_loss_fixed(isd:ied,jsd:jed))
2289 cs%TKE_itidal_loss_fixed = 0.0
2290 allocate(cs%TKE_leak_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2291 cs%TKE_leak_loss(:,:,:,:,:) = 0.0
2292 allocate(cs%TKE_quad_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2293 cs%TKE_quad_loss(:,:,:,:,:) = 0.0
2294 allocate(cs%TKE_itidal_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2295 cs%TKE_itidal_loss(:,:,:,:,:) = 0.0
2296 allocate(cs%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2297 cs%TKE_Froude_loss(:,:,:,:,:) = 0.0
2298 allocate(cs%tot_leak_loss(isd:ied,jsd:jed)) ; cs%tot_leak_loss(:,:) = 0.0
2299 allocate(cs%tot_quad_loss(isd:ied,jsd:jed) ) ; cs%tot_quad_loss(:,:) = 0.0
2300 allocate(cs%tot_itidal_loss(isd:ied,jsd:jed)) ; cs%tot_itidal_loss(:,:) = 0.0
2301 allocate(cs%tot_Froude_loss(isd:ied,jsd:jed)) ; cs%tot_Froude_loss(:,:) = 0.0
2304 call get_param(param_file, mdl,
"H2_FILE", h2_file, &
2305 "The path to the file containing the sub-grid-scale "//&
2306 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
2307 fail_if_missing=.true.)
2308 filename = trim(cs%inputdir) // trim(h2_file)
2309 call log_param(param_file, mdl,
"INPUTDIR/H2_FILE", filename)
2310 call mom_read_data(filename,
'h2', h2, g%domain, timelevel=1, scale=us%m_to_Z)
2311 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
2313 h2(i,j) = min(0.01*(g%bathyT(i,j))**2, h2(i,j))
2316 cs%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*gv%Rho0 * us%L_to_Z*kappa_itides * h2(i,j)
2322 call get_param(param_file, mdl,
"REFL_ANGLE_FILE", refl_angle_file, &
2323 "The path to the file containing the local angle of "//&
2324 "the coastline/ridge/shelf with respect to the equator.", &
2325 fail_if_missing=.false.)
2326 filename = trim(cs%inputdir) // trim(refl_angle_file)
2327 call log_param(param_file, mdl,
"INPUTDIR/REFL_ANGLE_FILE", filename)
2328 allocate(cs%refl_angle(isd:ied,jsd:jed)) ; cs%refl_angle(:,:) = cs%nullangle
2330 g%domain, timelevel=1)
2332 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
2333 if (is_nan(cs%refl_angle(i,j))) cs%refl_angle(i,j) = cs%nullangle
2335 call pass_var(cs%refl_angle,g%domain)
2338 call get_param(param_file, mdl,
"REFL_PREF_FILE", refl_pref_file, &
2339 "The path to the file containing the reflection coefficients.", &
2340 fail_if_missing=.false.)
2341 filename = trim(cs%inputdir) // trim(refl_pref_file)
2342 call log_param(param_file, mdl,
"INPUTDIR/REFL_PREF_FILE", filename)
2343 allocate(cs%refl_pref(isd:ied,jsd:jed)) ; cs%refl_pref(:,:) = 1.0
2344 call mom_read_data(filename,
'refl_pref', cs%refl_pref, g%domain, timelevel=1)
2346 call pass_var(cs%refl_pref,g%domain)
2349 allocate(cs%refl_pref_logical(isd:ied,jsd:jed)) ; cs%refl_pref_logical(:,:) = .false.
2353 if (cs%refl_angle(i,j) /= cs%nullangle .and. &
2354 cs%refl_pref(i,j) < 1.0 .and. cs%refl_pref(i,j) > 0.0)
then
2355 cs%refl_pref_logical(i,j) = .true.
2361 call get_param(param_file, mdl,
"REFL_DBL_FILE", refl_dbl_file, &
2362 "The path to the file containing the double-reflective ridge tags.", &
2363 fail_if_missing=.false.)
2364 filename = trim(cs%inputdir) // trim(refl_dbl_file)
2365 call log_param(param_file, mdl,
"INPUTDIR/REFL_DBL_FILE", filename)
2366 allocate(ridge_temp(isd:ied,jsd:jed)) ; ridge_temp(:,:) = 0.0
2367 call mom_read_data(filename,
'refl_dbl', ridge_temp, g%domain, timelevel=1)
2369 allocate(cs%refl_dbl(isd:ied,jsd:jed)) ; cs%refl_dbl(:,:) = .false.
2370 do i=isd,ied;
do j=jsd,jed
2371 if (ridge_temp(i,j) == 1) then; cs%refl_dbl(i,j) = .true.
2372 else ; cs%refl_dbl(i,j) = .false. ;
endif
2410 cs%id_refl_ang = register_diag_field(
'ocean_model',
'refl_angle', diag%axesT1, &
2411 time,
'Local angle of coastline/ridge/shelf with respect to equator',
'rad')
2412 cs%id_refl_pref = register_diag_field(
'ocean_model',
'refl_pref', diag%axesT1, &
2413 time,
'Partial reflection coefficients',
'')
2414 cs%id_dx_Cv = register_diag_field(
'ocean_model',
'dx_Cv', diag%axesT1, &
2415 time,
'North face unblocked width',
'm', conversion=us%L_to_m)
2416 cs%id_dy_Cu = register_diag_field(
'ocean_model',
'dy_Cu', diag%axesT1, &
2417 time,
'East face unblocked width',
'm', conversion=us%L_to_m)
2418 cs%id_land_mask = register_diag_field(
'ocean_model',
'land_mask', diag%axesT1, &
2419 time,
'Land mask',
'logical')
2421 if (cs%id_refl_ang > 0)
call post_data(cs%id_refl_ang, cs%refl_angle, cs%diag)
2422 if (cs%id_refl_pref > 0)
call post_data(cs%id_refl_pref, cs%refl_pref, cs%diag)
2423 if (cs%id_dx_Cv > 0)
call post_data(cs%id_dx_Cv, g%dx_Cv, cs%diag)
2424 if (cs%id_dy_Cu > 0)
call post_data(cs%id_dy_Cu, g%dy_Cu, cs%diag)
2425 if (cs%id_land_mask > 0)
call post_data(cs%id_land_mask, g%mask2dT, cs%diag)
2428 cs%id_tot_En = register_diag_field(
'ocean_model',
'ITide_tot_En', diag%axesT1, &
2429 time,
'Internal tide total energy density', &
2430 'J m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**2)
2432 cs%id_itide_drag = register_diag_field(
'ocean_model',
'ITide_drag', diag%axesT1, &
2433 time,
'Interior and bottom drag internal tide decay timescale',
's-1', conversion=us%s_to_T)
2435 cs%id_TKE_itidal_input = register_diag_field(
'ocean_model',
'TKE_itidal_input', diag%axesT1, &
2436 time,
'Conversion from barotropic to baroclinic tide, '//&
2437 'a fraction of which goes into rays', &
2438 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2440 cs%id_tot_leak_loss = register_diag_field(
'ocean_model',
'ITide_tot_leak_loss', diag%axesT1, &
2441 time,
'Internal tide energy loss to background drag', &
2442 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2443 cs%id_tot_quad_loss = register_diag_field(
'ocean_model',
'ITide_tot_quad_loss', diag%axesT1, &
2444 time,
'Internal tide energy loss to bottom drag', &
2445 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2446 cs%id_tot_itidal_loss = register_diag_field(
'ocean_model',
'ITide_tot_itidal_loss', diag%axesT1, &
2447 time,
'Internal tide energy loss to wave drag', &
2448 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2449 cs%id_tot_Froude_loss = register_diag_field(
'ocean_model',
'ITide_tot_Froude_loss', diag%axesT1, &
2450 time,
'Internal tide energy loss to wave breaking', &
2451 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2452 cs%id_tot_allprocesses_loss = register_diag_field(
'ocean_model',
'ITide_tot_allprocesses_loss', diag%axesT1, &
2453 time,
'Internal tide energy loss summed over all processes', &
2454 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2456 allocate(cs%id_En_mode(cs%nFreq,cs%nMode)) ; cs%id_En_mode(:,:) = -1
2457 allocate(cs%id_En_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_En_ang_mode(:,:) = -1
2458 allocate(cs%id_itidal_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_mode(:,:) = -1
2459 allocate(cs%id_allprocesses_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_allprocesses_loss_mode(:,:) = -1
2460 allocate(cs%id_itidal_loss_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_ang_mode(:,:) = -1
2461 allocate(cs%id_Ub_mode(cs%nFreq,cs%nMode)) ; cs%id_Ub_mode(:,:) = -1
2462 allocate(cs%id_cp_mode(cs%nFreq,cs%nMode)) ; cs%id_cp_mode(:,:) = -1
2464 allocate(angles(cs%NAngle)) ; angles(:) = 0.0
2465 angle_size = (8.0*atan(1.0)) / (real(num_angle))
2466 do a=1,num_angle ; angles(a) = (real(a) - 1) * angle_size ;
enddo
2468 id_ang = diag_axis_init(
"angle", angles,
"Radians",
"N",
"Angular Orienation of Fluxes")
2469 call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), axes_ang)
2471 do fr=1,cs%nFreq ;
write(freq_name(fr),
'("freq",i1)') fr ;
enddo
2472 do m=1,cs%nMode ;
do fr=1,cs%nFreq
2474 write(var_name,
'("Itide_En_freq",i1,"_mode",i1)') fr, m
2475 write(var_descript,
'("Internal tide energy density in frequency ",i1," mode ",i1)') fr, m
2476 cs%id_En_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2477 diag%axesT1, time, var_descript,
'J m-2', conversion=us%R_to_kg_m3*us%Z_to_m**2*us%s_to_T**3)
2478 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2481 write(var_name,
'("Itide_En_ang_freq",i1,"_mode",i1)') fr, m
2482 write(var_descript,
'("Internal tide angular energy density in frequency ",i1," mode ",i1)') fr, m
2483 cs%id_En_ang_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2484 axes_ang, time, var_descript,
'J m-2 band-1', conversion=us%R_to_kg_m3*us%Z_to_m**2*us%s_to_T**3)
2485 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2489 write(var_name,
'("Itide_wavedrag_loss_freq",i1,"_mode",i1)') fr, m
2490 write(var_descript,
'("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
2491 cs%id_itidal_loss_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2492 diag%axesT1, time, var_descript,
'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2493 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2495 write(var_name,
'("Itide_allprocesses_loss_freq",i1,"_mode",i1)') fr, m
2496 write(var_descript,
'("Internal tide energy loss due to all processes from frequency ",i1," mode ",i1)') fr, m
2497 cs%id_allprocesses_loss_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2498 diag%axesT1, time, var_descript,
'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2499 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2503 write(var_name,
'("Itide_wavedrag_loss_ang_freq",i1,"_mode",i1)') fr, m
2504 write(var_descript,
'("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
2505 cs%id_itidal_loss_ang_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2506 axes_ang, time, var_descript,
'W m-2 band-1', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2507 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2510 write(var_name,
'("Itide_Ub_freq",i1,"_mode",i1)') fr, m
2511 write(var_descript,
'("Near-bottom horizonal velocity for frequency ",i1," mode ",i1)') fr, m
2512 cs%id_Ub_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2513 diag%axesT1, time, var_descript,
'm s-1', conversion=us%L_T_to_m_s)
2514 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2517 write(var_name,
'("Itide_cp_freq",i1,"_mode",i1)') fr, m
2518 write(var_descript,
'("Horizonal phase velocity for frequency ",i1," mode ",i1)') fr, m
2519 cs%id_cp_mode(fr,m) = register_diag_field(
'ocean_model', var_name, &
2520 diag%axesT1, time, var_descript,
'm s-1', conversion=us%L_T_to_m_s)
2521 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
2535 if (
associated(cs))
then
2536 if (
associated(cs%En))
deallocate(cs%En)
2537 if (
allocated(cs%frequency))
deallocate(cs%frequency)
2538 if (
allocated(cs%id_En_mode))
deallocate(cs%id_En_mode)
2539 if (
allocated(cs%id_Ub_mode))
deallocate(cs%id_Ub_mode)
2540 if (
allocated(cs%id_cp_mode))
deallocate(cs%id_cp_mode)