From 6fd0efb277b09a8506016c1eaee3ae1e23f21233 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 30 Nov 2024 11:21:17 -0500 Subject: [PATCH] Fix rescaling of internal tide of debugging code Corrected the internal tide rescaling arguments so that some of the debugging variables have units that are consistent with their documented units. This involved changing the scale arguments to global_area_integral to tmp_scale arguments in 8 places so that the units of the output retain the scaling of the input variable. The multiplication by the reverse of the scaling factor was also added to 4 debugging output messages. Also replaced the scale argument to 35 chksum calls with the equivalent but preferred unscale argument, following the pattern elsewhere in the MOM6 code. All answers are bitwise identical, and only debugging code was modified. --- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- .../lateral/MOM_internal_tides.F90 | 100 ++++++++---------- .../vertical/MOM_set_diffusivity.F90 | 22 ++-- 3 files changed, 57 insertions(+), 67 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 4bac09b7cc..052fafbddc 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -269,7 +269,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (allocated(MEKE%mom_src)) & call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%mom_src_bh)) & - call hchksum(MEKE%mom_src_bh, 'MEKE mom_src_bh', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + call hchksum(MEKE%mom_src_bh, 'MEKE mom_src_bh', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%GME_snk)) & call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%GM_src)) & diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index c7101ac6b7..df6f9c8b6f 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -147,7 +147,7 @@ module MOM_internal_tides !! vertical profile squared, for each mode [H T-2 ~> m s-2 or kg m-2 s-2] real :: q_itides !< fraction of local dissipation [nondim] real :: mixing_effic !< mixing efficiency [nondim] - real :: En_sum !< global sum of energy for use in debugging, in MKS units [H Z2 T-2 L2 ~> m5 s-2 or J] + real :: En_sum !< global sum of energy for use in debugging, in MKS units [m5 s-2 or J] real :: En_underflow !< A minuscule amount of energy [H Z2 T-2 ~> m3 s-2 or J m-2] integer :: En_restart_power !< A power factor of 2 by which to multiply the energy in restart [nondim] type(time_type), pointer :: Time => NULL() !< A pointer to the model's clock. @@ -327,12 +327,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C real :: En_sumtmp ! Energies for debugging [H Z2 T-2 ~> m3 s-2 or J m-2] real :: En_initial, Delta_E_check ! Energies for debugging [H Z2 T-2 ~> m3 s-2 or J m-2] real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! Energy losses for debugging [H Z2 T-3 ~> m3 s-3 or W m-2] - real :: HZ2_T3_to_W_m2 ! unit conversion factor for TKE from internal to mks - ! [H Z2 T-3 ~> m3 s-3 or W m-2] real :: HZ2_T2_to_J_m2 ! unit conversion factor for Energy from internal to mks ! [H Z2 T-2 ~> m3 s-2 or J m-2] - real :: W_m2_to_HZ2_T3 ! unit conversion factor for TKE from mks to internal - ! [m3 s-3 or W m-2 ~> H Z2 T-3] real :: J_m2_to_HZ2_T2 ! unit conversion factor for Energy from mks to internal ! [m3 s-2 or J m-2 ~> H Z2 T-2] character(len=160) :: mesg ! The text of an error message @@ -346,9 +342,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nAngle = CS%NAngle - HZ2_T3_to_W_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**3) HZ2_T2_to_J_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**2) - W_m2_to_HZ2_T3 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**3) J_m2_to_HZ2_T2 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**2) cn_subRO = 1e-30*US%m_s_to_L_T @@ -402,7 +396,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C do m=1,CS%nMode ; do fr=1,CS%nFreq En_sumtmp = 0. do a=1,CS%nAngle - En_sumtmp = En_sumtmp + global_area_integral(CS%En(:,:,a,fr,m), G, scale=HZ2_T2_to_J_m2) + En_sumtmp = En_sumtmp + global_area_integral(CS%En(:,:,a,fr,m), G, tmp_scale=HZ2_T2_to_J_m2) enddo CS%En_ini_glo(fr,m) = En_sumtmp enddo ; enddo @@ -424,14 +418,14 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call pass_var(cn,G%Domain) if (CS%debug) then - call hchksum(cn(:,:,1), "CN mode 1", G%HI, haloshift=0, scale=US%L_to_m*US%s_to_T) + call hchksum(cn(:,:,1), "CN mode 1", G%HI, haloshift=0, unscale=US%L_to_m*US%s_to_T) call hchksum(CS%w_struct(:,:,:,1), "Wstruct mode 1", G%HI, haloshift=0) - call hchksum(CS%u_struct(:,:,:,1), "Ustruct mode 1", G%HI, haloshift=0, scale=US%m_to_Z) - call hchksum(CS%u_struct_bot(:,:,1), "Ustruct_bot mode 1", G%HI, haloshift=0, scale=US%m_to_Z) - call hchksum(CS%u_struct_max(:,:,1), "Ustruct_max mode 1", G%HI, haloshift=0, scale=US%m_to_Z) - call hchksum(CS%int_w2(:,:,1), "int_w2", G%HI, haloshift=0, scale=GV%H_to_MKS) - call hchksum(CS%int_U2(:,:,1), "int_U2", G%HI, haloshift=0, scale=GV%H_to_mks*US%m_to_Z**2) - call hchksum(CS%int_N2w2(:,:,1), "int_N2w2", G%HI, haloshift=0, scale=GV%H_to_mks*US%s_to_T**2) + call hchksum(CS%u_struct(:,:,:,1), "Ustruct mode 1", G%HI, haloshift=0, unscale=US%m_to_Z) + call hchksum(CS%u_struct_bot(:,:,1), "Ustruct_bot mode 1", G%HI, haloshift=0, unscale=US%m_to_Z) + call hchksum(CS%u_struct_max(:,:,1), "Ustruct_max mode 1", G%HI, haloshift=0, unscale=US%m_to_Z) + call hchksum(CS%int_w2(:,:,1), "int_w2", G%HI, haloshift=0, unscale=GV%H_to_MKS) + call hchksum(CS%int_U2(:,:,1), "int_U2", G%HI, haloshift=0, unscale=GV%H_to_mks*US%m_to_Z**2) + call hchksum(CS%int_N2w2(:,:,1), "int_N2w2", G%HI, haloshift=0, unscale=GV%H_to_mks*US%s_to_T**2) endif ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.********************** @@ -449,8 +443,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%debug) then call hchksum(TKE_itidal_input(:,:,1), "TKE_itidal_input", G%HI, haloshift=0, & - scale=GV%H_to_m*(US%Z_to_m**2)*(US%s_to_T)**3) - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides bf input", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + unscale=GV%H_to_m*(US%Z_to_m**2)*(US%s_to_T)**3) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides bf input", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) endif if (CS%energized_angle <= 0) then @@ -489,14 +483,14 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%init_forcing_only) CS%add_tke_forcing=.false. if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af input", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af input", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) ! save forcing for online budget do m=1,CS%nMode ; do fr=1,CS%nFreq En_sumtmp = 0. do a=1,CS%nAngle En_sumtmp = En_sumtmp + global_area_integral(dt*frac_per_sector*(1.0-CS%q_itides)* & CS%fraction_tidal_input(fr,m)*TKE_itidal_input(:,:,fr), & - G, scale=HZ2_T2_to_J_m2) + G, tmp_scale=HZ2_T2_to_J_m2) enddo CS%TKE_input_glo_dt(fr,m) = En_sumtmp enddo ; enddo @@ -508,7 +502,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call start_group_pass(pass_test, G%domain) if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after forcing') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after forcing', CS%En_sum @@ -535,7 +529,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after 1/2 refraction') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after 1/2 refraction', CS%En_sum @@ -546,7 +540,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging write(mesg,*) 'After first refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=',CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg)) !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") endif @@ -565,7 +559,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call correct_halo_rotation(CS%En, test, G, CS%nAngle, halo=En_halo_ij_stencil) if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo R", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo R", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after correct halo rotation') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after correct halo rotation', CS%En_sum @@ -596,7 +590,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af prop", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af prop", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after propagate') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after propagate', CS%En_sum @@ -608,7 +602,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset if (abs(CS%En(i,j,a,fr,m))>CS%En_check_tol) then ! only print if large write(mesg,*) 'After propagation: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=', CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg)) ! RD propagate produces very little negative energy (diff 2 large numbers), needs fix !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") @@ -638,7 +632,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr2", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr2", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after 2/2 refraction') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after 2/2 refraction', CS%En_sum @@ -649,7 +643,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging write(mesg,*) 'After second refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=', CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg)) !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") endif @@ -694,7 +688,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after leak", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after leak", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after background drag') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after background drag', CS%En_sum @@ -707,7 +701,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging write(mesg,*) 'After leak loss: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=',CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg), all_print=.true.) !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") endif @@ -718,7 +712,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C En_sumtmp = 0. do a=1,CS%nAngle En_sumtmp = En_sumtmp + global_area_integral(CS%TKE_leak_loss(:,:,a,fr,m)*dt, G, & - scale=HZ2_T2_to_J_m2) + tmp_scale=HZ2_T2_to_J_m2) enddo CS%TKE_leak_loss_glo_dt(fr,m) = En_sumtmp enddo ; enddo @@ -754,8 +748,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C enddo ; enddo ; enddo ; enddo endif - if (CS%debug) call hchksum(drag_scale(:,:,1,1), "dragscale", G%HI, haloshift=0, scale=US%s_to_T) - if (CS%debug) call hchksum(tot_vel_btTide2(:,:), "tot_vel_btTide2", G%HI, haloshift=0, scale=US%L_T_to_m_s**2) + if (CS%debug) call hchksum(drag_scale(:,:,1,1), "dragscale", G%HI, haloshift=0, unscale=US%s_to_T) + if (CS%debug) call hchksum(tot_vel_btTide2(:,:), "tot_vel_btTide2", G%HI, haloshift=0, unscale=US%L_T_to_m_s**2) do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale @@ -778,13 +772,13 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after quad", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after quad", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) ! save loss term for online budget do m=1,CS%nMode ; do fr=1,CS%nFreq En_sumtmp = 0. do a=1,CS%nAngle En_sumtmp = En_sumtmp + global_area_integral(CS%TKE_quad_loss(:,:,a,fr,m)*dt, G, & - scale=HZ2_T2_to_J_m2) + tmp_scale=HZ2_T2_to_J_m2) enddo CS%TKE_quad_loss_glo_dt(fr,m) = En_sumtmp enddo ; enddo @@ -794,7 +788,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging write(mesg,*) 'After bottom loss: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=',CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg), all_print=.true.) !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") endif @@ -865,7 +859,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after wave", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after wave", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: before Froude drag') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: before Froude drag', CS%En_sum @@ -875,7 +869,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C En_sumtmp = 0. do a=1,CS%nAngle En_sumtmp = En_sumtmp + global_area_integral(CS%TKE_itidal_loss(:,:,a,fr,m)*dt, G, & - scale=HZ2_T2_to_J_m2) + tmp_scale=HZ2_T2_to_J_m2) enddo CS%TKE_itidal_loss_glo_dt(fr,m) = En_sumtmp enddo ; enddo @@ -885,7 +879,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging write(mesg,*) 'After wave drag loss: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=',CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg), all_print=.true.) !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") endif @@ -939,7 +933,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after froude", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after froude", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after Froude drag') if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after Froude drag', CS%En_sum @@ -951,7 +945,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C En_sumtmp = 0. do a=1,CS%nAngle En_sumtmp = En_sumtmp + global_area_integral(CS%TKE_Froude_loss(:,:,a,fr,m)*dt, G, & - scale=HZ2_T2_to_J_m2) + tmp_scale=HZ2_T2_to_J_m2) enddo CS%TKE_Froude_loss_glo_dt(fr,m) = En_sumtmp enddo ; enddo @@ -961,7 +955,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset write(mesg,*) 'After Froude loss: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=',CS%En(i,j,a,fr,m) + 'En=', HZ2_T2_to_J_m2*CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg), all_print=.true.) !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") endif @@ -994,7 +988,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C endif if (CS%debug) then - call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after slope", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) + call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after slope", G%HI, haloshift=0, unscale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide') enddo ; enddo @@ -1003,7 +997,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C En_sumtmp = 0. do a=1,CS%nAngle En_sumtmp = En_sumtmp + global_area_integral(CS%TKE_residual_loss(:,:,a,fr,m)*dt, G, & - scale=HZ2_T2_to_J_m2) + tmp_scale=HZ2_T2_to_J_m2) enddo CS%TKE_residual_loss_glo_dt(fr,m) = En_sumtmp enddo ; enddo @@ -1015,7 +1009,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C do m=1,CS%nMode ; do fr=1,CS%nFreq En_sumtmp = 0. do a=1,CS%nAngle - En_sumtmp = En_sumtmp + global_area_integral(CS%En(:,:,a,fr,m), G, scale=HZ2_T2_to_J_m2) + En_sumtmp = En_sumtmp + global_area_integral(CS%En(:,:,a,fr,m), G, tmp_scale=HZ2_T2_to_J_m2) enddo CS%En_end_glo(fr,m) = En_sumtmp enddo ; enddo @@ -1025,7 +1019,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C CS%TKE_quad_loss_glo_dt(fr,m) - CS%TKE_itidal_loss_glo_dt(fr,m) - & CS%TKE_Froude_loss_glo_dt(fr,m) - CS%TKE_residual_loss_glo_dt(fr,m) - & CS%En_end_glo(fr,m) - if (is_root_pe()) write(stdout,'(A,F18.10)') "error in Energy budget", CS%error_mode(fr,m) + if (is_root_pe()) write(stdout,'(A,F18.10)') "error in Energy budget", HZ2_T2_to_J_m2*CS%error_mode(fr,m) enddo ; enddo endif @@ -1283,12 +1277,10 @@ subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixe real :: En_a, En_b ! energy before and after timestep [H Z2 T-2 ~> m3 s-2 or J m-2] real :: I_dt ! The inverse of the timestep [T-1 ~> s-1] real :: J_m2_to_HZ2_T2 ! unit conversion factor for Energy from mks to internal [m3 s-2 or J m-2 ~> H Z2 T-2] - real :: HZ2_T3_to_W_m2 ! unit conversion factor for Energy from internal to mks [H Z2 T-3 ~> m3 s-3 or W m-2] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec J_m2_to_HZ2_T2 = GV%m_to_H*(US%m_to_Z**2)*(US%T_to_s**2) - HZ2_T3_to_W_m2 = GV%H_to_m*(US%Z_to_m**2)*(US%s_to_T**3) I_dt = 1.0 / dt q_itides = CS%q_itides @@ -1301,9 +1293,9 @@ subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixe if (CS%debug) then call hchksum(TKE_loss_fixed, "TKE loss fixed", G%HI, haloshift=0, & - scale=US%RZ_to_kg_m2*(US%Z_to_m**3)*GV%m_to_H*(US%m_to_L**2)) - call hchksum(Nb(:,:), "Nbottom", G%HI, haloshift=0, scale=US%s_to_T) - call hchksum(Ub(:,:,1,1), "Ubottom", G%HI, haloshift=0, scale=US%L_to_m*US%s_to_T) + unscale=US%RZ_to_kg_m2*(US%Z_to_m**3)*GV%m_to_H*(US%m_to_L**2)) + call hchksum(Nb(:,:), "Nbottom", G%HI, haloshift=0, unscale=US%s_to_T) + call hchksum(Ub(:,:,1,1), "Ubottom", G%HI, haloshift=0, unscale=US%L_to_m*US%s_to_T) endif do j=js,je ; do i=is,ie ; do m=1,CS%nMode ; do fr=1,CS%nFreq @@ -2284,7 +2276,7 @@ subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS !energized_angle = Angle_size * real(energized_wedge - 1) + 2.0*Angle_size ! !energized_angle = Angle_size * real(energized_wedge - 1) + 0.5*Angle_size ! do J=jsh-1,jeh ; do I=ish-1,ieh - ! This will only work for a Cartesian grid for which G%geoLonBu is in the same units has dx. + !### This will only work for a Cartesian grid for which G%geoLonBu is in the same units has dx. ! This needs to be extensively revised to work for a general grid. x(I,J) = G%US%m_to_L*G%geoLonBu(I,J) y(I,J) = G%US%m_to_L*G%geoLatBu(I,J) @@ -2867,7 +2859,7 @@ subroutine reflect(En, NAngle, CS, G, LB) ! Check to make sure no energy gets onto land (only run for debugging) ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then - ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',i+G%idg_offset, 'jg_g=',j+G%jdg_offset + ! write (mesg,*) 'En=', HZ2_T2_to_J_m2*En(i,j,a), 'a=', a, 'ig_g=',i+G%idg_offset, 'jg_g=',j+G%jdg_offset ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.) ! endif ! enddo ; enddo ; enddo @@ -3397,7 +3389,6 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) real :: period ! A tidal period read from namelist [T ~> s] real :: HZ2_T2_to_J_m2 ! unit conversion factor for Energy from internal to mks [H Z2 T-2 ~> m3 s-2 or J m-2] real :: HZ2_T3_to_W_m2 ! unit conversion factor for TKE from internal to mks [H Z2 T-3 ~> m3 s-3 or W m-2] - real :: W_m2_to_HZ2_T3 ! unit conversion factor for TKE from mks to internal [m3 s-3 or W m-2 ~> H Z2 T-3] real :: J_m2_to_HZ2_T2 ! unit conversion factor for Energy from mks to internal [m3 s-2 or J m-2 ~> H Z2 T-2] integer :: num_angle, num_freq, num_mode, m, fr integer :: isd, ied, jsd, jed, a, id_ang, i, j, nz @@ -3422,7 +3413,6 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) HZ2_T2_to_J_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**2) HZ2_T3_to_W_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**3) - W_m2_to_HZ2_T3 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**3) J_m2_to_HZ2_T2 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**2) CS%initialized = .true. @@ -3545,7 +3535,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=0.2) call get_param(param_file, mdl, "MAX_TKE_TO_KD", CS%max_TKE_to_Kd, & "Limiter for TKE_to_Kd.", & - units="", default=1e9, scale=US%Z_to_m*US%s_to_T**2) + units="s2 m-1", default=1e9, scale=US%Z_to_m*US%s_to_T**2) call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", CS%decay_rate, & "The rate at which internal tide energy is lost to the "//& "interior ocean internal wave field.", & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index d15a364267..85cf9fdf7a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -716,17 +716,17 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i endif if (CS%debug) then - if (CS%id_prof_leak > 0) call hchksum(dd%prof_leak, "leakage_profile", G%HI, haloshift=0, scale=GV%m_to_H) - if (CS%id_prof_slope > 0) call hchksum(dd%prof_slope, "slope_profile", G%HI, haloshift=0, scale=GV%m_to_H) - if (CS%id_prof_Froude > 0) call hchksum(dd%prof_Froude, "Froude_profile", G%HI, haloshift=0, scale=GV%m_to_H) - if (CS%id_prof_quad > 0) call hchksum(dd%prof_quad, "quad_profile", G%HI, haloshift=0, scale=GV%m_to_H) - if (CS%id_prof_itidal > 0) call hchksum(dd%prof_itidal, "itidal_profile", G%HI, haloshift=0, scale=GV%m_to_H) - if (CS%id_TKE_to_Kd > 0) call hchksum(dd%TKE_to_Kd, "TKE_to_Kd", G%HI, haloshift=0, scale=US%m_to_Z*US%T_to_s**2) - if (CS%id_Kd_leak > 0) call hchksum(dd%Kd_leak, "Kd_leak", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) - if (CS%id_Kd_quad > 0) call hchksum(dd%Kd_quad, "Kd_quad", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) - if (CS%id_Kd_itidal > 0) call hchksum(dd%Kd_itidal, "Kd_itidal", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) - if (CS%id_Kd_Froude > 0) call hchksum(dd%Kd_Froude, "Kd_Froude", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) - if (CS%id_Kd_slope > 0) call hchksum(dd%Kd_slope, "Kd_slope", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) + if (CS%id_prof_leak > 0) call hchksum(dd%prof_leak, "leakage_profile", G%HI, haloshift=0, unscale=GV%m_to_H) + if (CS%id_prof_slope > 0) call hchksum(dd%prof_slope, "slope_profile", G%HI, haloshift=0, unscale=GV%m_to_H) + if (CS%id_prof_Froude > 0) call hchksum(dd%prof_Froude, "Froude_profile", G%HI, haloshift=0, unscale=GV%m_to_H) + if (CS%id_prof_quad > 0) call hchksum(dd%prof_quad, "quad_profile", G%HI, haloshift=0, unscale=GV%m_to_H) + if (CS%id_prof_itidal > 0) call hchksum(dd%prof_itidal, "itidal_profile", G%HI, haloshift=0, unscale=GV%m_to_H) + if (CS%id_TKE_to_Kd > 0) call hchksum(dd%TKE_to_Kd, "TKE_to_Kd", G%HI, haloshift=0, unscale=US%m_to_Z*US%T_to_s**2) + if (CS%id_Kd_leak > 0) call hchksum(dd%Kd_leak, "Kd_leak", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) + if (CS%id_Kd_quad > 0) call hchksum(dd%Kd_quad, "Kd_quad", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) + if (CS%id_Kd_itidal > 0) call hchksum(dd%Kd_itidal, "Kd_itidal", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) + if (CS%id_Kd_Froude > 0) call hchksum(dd%Kd_Froude, "Kd_Froude", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) + if (CS%id_Kd_slope > 0) call hchksum(dd%Kd_slope, "Kd_slope", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) endif ! post diagnostics