From c86358b747c4894719b454637cc21406d57849a6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 20 Jan 2024 06:26:23 -0500 Subject: [PATCH] (*)Avoiding using RHO_0 in non-Boussinesq debugging When in non-Boussinesq mode, write out mass-based tracer inventories from calls to MOM_tracer_chkinv and mean thicknesses in mass per unit area from calls to MOM_state_stats, rather than the approximate volumes that depend on the Boussinesq reference density. Similarly the thicknesses and transports output by write_u_accel and write_v_accel are in mass per unit area or mass flux per unit face length in non-Boussinesq mode. This is done specifically by using GV%H_to_MKS in place of GV%H_to_m; these are identical in Boussinesq mode, but GV%H_to_MKS avoids rescaling by the Boussinesq reference density in non-Boussinesq mode. Several comments were also updated to reflect these changes. All solutions are identical, but there are changes in the units of several diagnostic output fields that the model can write out when debugging non-Boussinesq mode runs. --- src/core/MOM_checksum_packages.F90 | 13 +++++++------ src/diagnostics/MOM_PointAccel.F90 | 18 ++++++++++-------- src/tracer/MOM_tracer_registry.F90 | 20 ++++++++++++-------- 3 files changed, 29 insertions(+), 22 deletions(-) diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index 4a9df04c4d..dc60e0888f 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -268,10 +268,11 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: & tmp_A, & ! The area per cell [m2] (unscaled to permit reproducing sum). - tmp_V, & ! The column-integrated volume [m3] (unscaled to permit reproducing sum) - tmp_T, & ! The column-integrated temperature [degC m3] (unscaled to permit reproducing sum) - tmp_S ! The column-integrated salinity [ppt m3] (unscaled to permit reproducing sum) - real :: Vol, dV ! The total ocean volume and its change [m3] (unscaled to permit reproducing sum). + tmp_V, & ! The column-integrated volume [m3] or mass [kg] (unscaled to permit reproducing sum), + ! depending on whether the Boussinesq approximation is used + tmp_T, & ! The column-integrated temperature [degC m3] or [degC kg] (unscaled to permit reproducing sum) + tmp_S ! The column-integrated salinity [ppt m3] or [ppt kg] (unscaled to permit reproducing sum) + real :: Vol, dV ! The total ocean volume or mass and its change [m3] or [kg] (unscaled to permit reproducing sum). real :: Area ! The total ocean surface area [m2] (unscaled to permit reproducing sum). real :: h_minimum ! The minimum layer thicknesses [H ~> m or kg m-2] real :: T_scale ! The scaling conversion factor for temperatures [degC C-1 ~> 1] @@ -284,7 +285,7 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe ! assumption we will not turn this on with threads type(stats), save :: oldT, oldS logical, save :: firstCall = .true. - real, save :: oldVol ! The previous total ocean volume [m3] + real, save :: oldVol ! The previous total ocean volume [m3] or mass [kg] character(len=80) :: lMsg integer :: is, ie, js, je, nz, i, j, k @@ -308,7 +309,7 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe h_minimum = 1.E34*GV%m_to_H do k=1,nz ; do j=js,je ; do i=is,ie if (G%mask2dT(i,j)>0.) then - dV = US%L_to_m**2*G%areaT(i,j)*GV%H_to_m*h(i,j,k) + dV = US%L_to_m**2*G%areaT(i,j)*GV%H_to_MKS*h(i,j,k) tmp_V(i,j) = tmp_V(i,j) + dV if (do_TS .and. h(i,j,k)>0.) then T%minimum = min( T%minimum, T_scale*Temp(i,j,k) ) ; T%maximum = max( T%maximum, T_scale*Temp(i,j,k) ) diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index e9c1092ed7..ab3c104d0f 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -95,9 +95,10 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st real :: du ! A velocity change [L T-1 ~> m s-1] real :: Inorm(SZK_(GV)) ! The inverse of the normalized velocity change [L T-1 ~> m s-1] real :: e(SZK_(GV)+1) ! Simple estimates of interface heights based on the sum of thicknesses [m] - real :: h_scale ! A scaling factor for thicknesses [m H-1 ~> 1 or m3 kg-1] + real :: h_scale ! A scaling factor for thicknesses [m H-1 ~> 1] or [kg m-2 H-1 ~> 1] real :: vel_scale ! A scaling factor for velocities [m T s-1 L-1 ~> 1] - real :: uh_scale ! A scaling factor for transport per unit length [m2 T s-1 L-1 H-1 ~> 1 or m3 kg-1] + real :: uh_scale ! A scaling factor for transport per unit length [m2 T s-1 L-1 H-1 ~> 1] + ! or [kg T m-1 s-1 L-1 H-1 ~> 1] real :: temp_scale ! A scaling factor for temperatures [degC C-1 ~> 1] real :: saln_scale ! A scaling factor for salinities [ppt S-1 ~> 1] integer :: yr, mo, day, hr, minute, sec, yearday @@ -108,7 +109,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st integer :: file Angstrom = GV%Angstrom_H + GV%H_subroundoff - h_scale = GV%H_to_m ; vel_scale = US%L_T_to_m_s ; uh_scale = GV%H_to_m*US%L_T_to_m_s + h_scale = GV%H_to_mks ; vel_scale = US%L_T_to_m_s ; uh_scale = h_scale*vel_scale temp_scale = US%C_to_degC ; saln_scale = US%S_to_ppt ! if (.not.associated(CS)) return @@ -232,7 +233,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(I,j,k) ; enddo endif if (present(str)) then - write(file,'(/,"Stress: ",ES10.3)', advance='no') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) + write(file,'(/,"Stress: ",ES10.3)', advance='no') (uh_scale*GV%RZ_to_H) * (str*dt) endif if (associated(CS%u_accel_bt)) then @@ -435,9 +436,10 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st real :: dv ! A velocity change [L T-1 ~> m s-1] real :: Inorm(SZK_(GV)) ! The inverse of the normalized velocity change [L T-1 ~> m s-1] real :: e(SZK_(GV)+1) ! Simple estimates of interface heights based on the sum of thicknesses [m] - real :: h_scale ! A scaling factor for thicknesses [m H-1 ~> 1 or m3 kg-1] + real :: h_scale ! A scaling factor for thicknesses [m H-1 ~> 1] or [kg m-2 H-1 ~> 1] real :: vel_scale ! A scaling factor for velocities [m T s-1 L-1 ~> 1] - real :: uh_scale ! A scaling factor for transport per unit length [m2 T s-1 L-1 H-1 ~> 1 or m3 kg-1] + real :: uh_scale ! A scaling factor for transport per unit length [m2 T s-1 L-1 H-1 ~> 1] + ! or [kg T m-1 s-1 L-1 H-1 ~> 1] real :: temp_scale ! A scaling factor for temperatures [degC C-1 ~> 1] real :: saln_scale ! A scaling factor for salinities [ppt S-1 ~> 1] integer :: yr, mo, day, hr, minute, sec, yearday @@ -448,7 +450,7 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st integer :: file Angstrom = GV%Angstrom_H + GV%H_subroundoff - h_scale = GV%H_to_m ; vel_scale = US%L_T_to_m_s ; uh_scale = GV%H_to_m*US%L_T_to_m_s + h_scale = GV%H_to_mks ; vel_scale = US%L_T_to_m_s ; uh_scale = h_scale*vel_scale temp_scale = US%C_to_degC ; saln_scale = US%S_to_ppt ! if (.not.associated(CS)) return @@ -576,7 +578,7 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(i,J,k) ; enddo endif if (present(str)) then - write(file,'(/,"Stress: ",ES10.3)', advance='no') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) + write(file,'(/,"Stress: ",ES10.3)', advance='no') (uh_scale*GV%RZ_to_H) * (str*dt) endif if (associated(CS%v_accel_bt)) then diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index c01419f3f8..8f7f2a280c 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -786,14 +786,16 @@ subroutine tracer_array_chkinv(mesg, G, GV, h, Tr, ntr) integer, intent(in) :: ntr !< number of registered tracers ! Local variables - real :: vol_scale ! The dimensional scaling factor to convert volumes to m3 [m3 H-1 L-2 ~> 1 or m3 kg-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tr_inv ! Volumetric tracer inventory in each cell [conc m3] - real :: total_inv ! The total amount of tracer [conc m3] + real :: vol_scale ! The dimensional scaling factor to convert volumes to m3 [m3 H-1 L-2 ~> 1] or cell + ! masses to kg [kg H-1 L-2 ~> 1], depending on whether the Boussinesq approximation is used + real :: tr_inv(SZI_(G),SZJ_(G),SZK_(GV)) ! Volumetric or mass-based tracer inventory in + ! each cell [conc m3] or [conc kg] + real :: total_inv ! The total amount of tracer [conc m3] or [conc kg] integer :: is, ie, js, je, nz integer :: i, j, k, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - vol_scale = GV%H_to_m*G%US%L_to_m**2 + vol_scale = GV%H_to_MKS*G%US%L_to_m**2 do m=1,ntr do k=1,nz ; do j=js,je ; do i=is,ie tr_inv(i,j,k) = Tr(m)%conc_scale*Tr(m)%t(i,j,k) * (vol_scale * h(i,j,k) * G%areaT(i,j)*G%mask2dT(i,j)) @@ -814,16 +816,18 @@ subroutine tracer_Reg_chkinv(mesg, G, GV, h, Reg) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] ! Local variables - real :: vol_scale ! The dimensional scaling factor to convert volumes to m3 [m3 H-1 L-2 ~> 1 or m3 kg-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tr_inv ! Volumetric tracer inventory in each cell [conc m3] - real :: total_inv ! The total amount of tracer [conc m3] + real :: vol_scale ! The dimensional scaling factor to convert volumes to m3 [m3 H-1 L-2 ~> 1] or cell + ! masses to kg [kg H-1 L-2 ~> 1], depending on whether the Boussinesq approximation is used + real :: tr_inv(SZI_(G),SZJ_(G),SZK_(GV)) ! Volumetric or mass-based tracer inventory in + ! each cell [conc m3] or [conc kg] + real :: total_inv ! The total amount of tracer [conc m3] or [conc kg] integer :: is, ie, js, je, nz integer :: i, j, k, m if (.not.associated(Reg)) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - vol_scale = GV%H_to_m*G%US%L_to_m**2 + vol_scale = GV%H_to_MKS*G%US%L_to_m**2 do m=1,Reg%ntr do k=1,nz ; do j=js,je ; do i=is,ie tr_inv(i,j,k) = Reg%Tr(m)%conc_scale*Reg%Tr(m)%t(i,j,k) * (vol_scale * h(i,j,k) * G%areaT(i,j)*G%mask2dT(i,j))