From 3acd777490645758a0844649dd3ad26a532adb1e Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 23 Dec 2020 16:01:44 -0500 Subject: [PATCH 01/31] Draft. Compiles but FPE in advection --- src/tracer/MOM_tracer_flow_control.F90 | 24 +- src/tracer/nw2_tracers.F90 | 319 +++++++++++++++++++++++++ 2 files changed, 338 insertions(+), 5 deletions(-) create mode 100644 src/tracer/nw2_tracers.F90 diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 4c7c27c7e6..c1885644cc 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -62,6 +62,8 @@ module MOM_tracer_flow_control use boundary_impulse_tracer, only : boundary_impulse_tracer_column_physics, boundary_impulse_tracer_surface_state use boundary_impulse_tracer, only : boundary_impulse_stock, boundary_impulse_tracer_end use boundary_impulse_tracer, only : boundary_impulse_tracer_CS +use nw2_tracers, only : nw2_tracers_CS, register_nw2_tracers, nw2_tracer_column_physics +use nw2_tracers, only : initialize_nw2_tracers, nw2_tracers_end implicit none ; private @@ -84,6 +86,7 @@ module MOM_tracer_flow_control logical :: use_pseudo_salt_tracer = .false. !< If true, use the psuedo_salt tracer package logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package logical :: use_dyed_obc_tracer = .false. !< If true, use the dyed OBC tracer package + logical :: use_nw2_tracers = .false. !< If true, use the ideal age tracer package !>@{ Pointers to the control strucures for the tracer packages type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL() type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL() @@ -98,6 +101,7 @@ module MOM_tracer_flow_control type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL() type(boundary_impulse_tracer_CS), pointer :: boundary_impulse_tracer_CSp => NULL() type(dyed_obc_tracer_CS), pointer :: dyed_obc_tracer_CSp => NULL() + type(nw2_tracers_CS), pointer :: nw2_tracers_CSp => NULL() !>@} end type tracer_flow_control_CS @@ -206,6 +210,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) call get_param(param_file, mdl, "USE_DYED_OBC_TRACER", CS%use_dyed_obc_tracer, & "If true, use the dyed_obc_tracer tracer package.", & default=.false.) + call get_param(param_file, mdl, "USE_NW2_TRACERS", CS%use_nw2_tracers, & + "If true, use the NeverWorld2 tracers.", & + default=.false.) ! Add other user-provided calls to register tracers for restarting here. Each ! tracer package registration call returns a logical false if it cannot be run @@ -249,7 +256,8 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) if (CS%use_dyed_obc_tracer) CS%use_dyed_obc_tracer = & register_dyed_obc_tracer(HI, GV, param_file, CS%dyed_obc_tracer_CSp, & tr_Reg, restart_CS) - + if (CS%use_nw2_tracers) CS%use_ideal_age = & + register_nw2_tracers(HI, GV, param_file, CS%nw2_tracers_CSp, tr_Reg, restart_CS) end subroutine call_tracer_register @@ -328,6 +336,8 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag sponge_CSp, tv) if (CS%use_dyed_obc_tracer) & call initialize_dyed_obc_tracer(restart, day, G, GV, h, diag, OBC, CS%dyed_obc_tracer_CSp) + if (CS%use_nw2_tracers) & + call initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS%nw2_tracers_CSp) end subroutine tracer_flow_control_init @@ -485,8 +495,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, G, GV, US, CS%dyed_obc_tracer_CSp, & evap_CFL_limit=evap_CFL_limit, & minimum_forcing_depth=minimum_forcing_depth) - - + if (CS%use_nw2_tracers) & + call nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, tv, CS%nw2_tracers_CSp, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) else ! Apply tracer surface fluxes using ea on the first layer if (CS%use_USER_tracer_example) & call tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & @@ -531,10 +544,10 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, if (CS%use_dyed_obc_tracer) & call dyed_obc_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, CS%dyed_obc_tracer_CSp) - + if (CS%use_nw2_tracers) call nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, tv, CS%nw2_tracers_CSp) endif - end subroutine call_tracer_column_fns !> This subroutine calls all registered tracer packages to enable them to @@ -774,6 +787,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp) if (CS%use_boundary_impulse_tracer) call boundary_impulse_tracer_end(CS%boundary_impulse_tracer_CSp) if (CS%use_dyed_obc_tracer) call dyed_obc_tracer_end(CS%dyed_obc_tracer_CSp) + if (CS%use_nw2_tracers) call nw2_tracers_end(CS%nw2_tracers_CSp) if (associated(CS)) deallocate(CS) end subroutine tracer_flow_control_end diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 new file mode 100644 index 0000000000..e71387f9c7 --- /dev/null +++ b/src/tracer/nw2_tracers.F90 @@ -0,0 +1,319 @@ +!> Ideal tracers designed to help diagnose a tracer diffusivity tensor in NeverWorld2 +module nw2_tracers + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_grid, only : ocean_grid_type +use MOM_hor_index, only : hor_index_type +use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_restart, only : query_initialized, MOM_restart_CS +use MOM_time_manager, only : time_type, time_type_to_real +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface, thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +public register_nw2_tracers +public initialize_nw2_tracers +public nw2_tracer_column_physics +public nw2_tracers_end + +integer, parameter :: NTR_MAX = 20 !< the maximum number of tracers in this module. + +!> The control structure for the nw2_tracers package +type, public :: nw2_tracers_CS ; private + integer :: ntr = NTR_MAX !< The number of tracers that are actually used. + type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. + type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry + real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this package, in g m-3? + real, dimension(NTR_MAX) :: restore_rate !< The exponential growth rate for restoration value [year-1]. + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart controls structure +end type nw2_tracers_CS + +contains + +!> Register the NW2 tracer fields to be used with MOM. +logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracer. + type(tracer_registry_type), pointer :: tr_Reg !< A pointer that is set to point to the control + !! structure for the tracer advection and + !! diffusion module + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure + +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "nw2_example" ! This module's name. + character(len=200) :: inputdir ! The directory where the input files are. + character(len=8) :: var_name ! The variable's name. + real, pointer :: tr_ptr(:,:,:) => NULL() + logical :: do_nw2 + integer :: isd, ied, jsd, jed, nz, m + type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + + if (associated(CS)) then + call MOM_error(FATAL, "register_nw2_tracer called with an "// & + "associated control structure.") + return + endif + allocate(CS) + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + + allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0 + + do m=1,CS%ntr + write(var_name(1:8),'(a6,i2.2)') 'tracer',m + tr_desc = var_desc(var_name, "1", "Ideal Tracer", caller=mdl) + ! This is needed to force the compiler not to do a copy in the registration + ! calls. Curses on the designers and implementers of Fortran90. + tr_ptr => CS%tr(:,:,:,m) + call query_vardesc(tr_desc, name=var_name, caller="register_nw2_tracers") + ! Register the tracer for horizontal advection, diffusion, and restarts. + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, & + registry_diags=.true., restart_CS=restart_CS, mandatory=.true.) + if (m<=10) then + CS%restore_rate(m) = 1.0 / 6.0 ! 6 years + else + CS%restore_rate(m) = 1.0 ! 1 year + endif + enddo + + CS%tr_Reg => tr_Reg + CS%restart_CSp => restart_CS + register_nw2_tracers = .true. +end function register_nw2_tracers + +!> Sets the NW2 traces to their initial values and sets up the tracer output +subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS) + logical, intent(in) :: restart !< .true. if the fields have already + !! been read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate + !! diagnostic output. + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracer. + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights + real :: rscl ! z* scaling factor + integer :: i, j, k, m + + if (.not.associated(CS)) return + + CS%Time => day + CS%diag => diag + + ! Calculate z* interface positions + if (GV%Boussinesq) then + ! First calculate interface positions in z-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,GV%ke+1) = - G%mask2dT(i,j) * G%bathyT(i,j) + enddo ; enddo + do k=GV%ke,1,-1 ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h(i,j,k) * GV%H_to_Z + enddo ; enddo ; enddo + ! Re-calculate for interface positions in z*-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (G%bathyT(i,j)>0.) then + rscl = G%bathyT(i,j) / ( eta(i,j,1) + G%bathyT(i,j) ) + do K=GV%ke, 1, -1 + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h(i,j,k) * GV%H_to_Z * rscl + enddo + endif + enddo ; enddo + else + call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode") + endif + + if (.not.restart) then + do m=1,CS%ntr + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! if (G%mask2dT(i,j) < 0.5) then + ! CS%tr(i,j,k,m) = 0. + ! else + CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k) + ! endif + enddo ; enddo ; enddo + enddo ! Tracer loop + endif ! restart + +end subroutine initialize_nw2_tracers + +!> Applies diapycnal diffusion, aging and regeneration at the surface to the NW2 tracers +subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, tv, CS, & + evap_CFL_limit, minimum_forcing_depth) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: ea !< an array to which the amount of fluid entrained + !! from the layer above during this call will be + !! added [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: eb !< an array to which the amount of fluid entrained + !! from the layer below during this call will be + !! added [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic + !! and tracer forcing fields. Unused fields have NULL ptrs. + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracer. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which + !! fluxes can be applied [H ~> m or kg m-2] +! This subroutine applies diapycnal diffusion and any other column +! tracer physics or chemistry to the tracers from this file. +! This is a simple example of a set of advected passive tracers. + +! The arguments to this subroutine are redundant in that +! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights + integer :: i, j, k, m + real :: dt_x_rate ! dt * restoring rate + real :: rscl ! z* scaling factor + real :: target_value ! tracer value + +! if (.not.associated(CS)) return + + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + do m=1,CS%ntr + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CS%tr(:,:,:,m), dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV) + enddo + else + do m=1,CS%ntr + call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV) + enddo + endif + + ! Calculate z* interface positions + if (GV%Boussinesq) then + ! First calculate interface positions in z-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,GV%ke+1) = - G%mask2dT(i,j) * G%bathyT(i,j) + enddo ; enddo + do k=GV%ke,1,-1 ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h_new(i,j,k) * GV%H_to_Z + enddo ; enddo ; enddo + ! Re-calculate for interface positions in z*-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (G%bathyT(i,j)>0.) then + rscl = G%bathyT(i,j) / ( eta(i,j,1) + G%bathyT(i,j) ) + do K=GV%ke, 1, -1 + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h_new(i,j,k) * GV%H_to_Z * rscl + enddo + endif + enddo ; enddo + else + call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode") + endif + + do m=1,CS%ntr + dt_x_rate = dt * ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) ) +!$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate) + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k) + CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + G%mask2dT(i,j) * dt_x_rate * ( target_value - CS%tr(i,j,k,m) ) + enddo ; enddo ; enddo + enddo + +end subroutine nw2_tracer_column_physics + +!> The target value of a NeverWorld2 tracer label m at non-dimensional +!! position x=lon/Lx, y=lat/Ly, z=eta/H +real function nw2_tracer_dist(m, G, GV, eta, i, j, k) + integer, intent(in) :: m !< Indicates the NW2 tracer + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),0:SZK_(G)), & + intent(in) :: eta !< Interface position [m] + integer, intent(in) :: i !< Cell index i + integer, intent(in) :: j !< Cell index j + integer, intent(in) :: k !< Layer index k + ! Local variables + real :: pi ! 3.1415... + real :: x, y, z ! non-dimensional positions + pi = 2.*acos(0.) + x = ( G%geolonT(i,j) - G%west_lon ) / G%len_lon ! 0 ... 1 + y = -G%geolatT(i,j) / G%south_lat ! -1 ... 1 + z = - 0.5 * ( eta(i,j,K-1) + eta(i,j,K) ) / GV%max_depth ! 0 ... 1 + select case ( mod(m-1,10) ) + case (0) ! y/L + nw2_tracer_dist = y + case (1) ! -z/L + nw2_tracer_dist = -z + case (2) ! cos(2 pi x/L) + nw2_tracer_dist = cos( 2.0 * pi * x ) + case (3) ! sin(2 pi x/L) + nw2_tracer_dist = sin( 2.0 * pi * x ) + case (4) ! sin(4 pi x/L) + nw2_tracer_dist = sin( 2.0 * pi * mod( 2.0*x, 1.0 ) ) + case (5) ! sin(pi y/L) + nw2_tracer_dist = sin( pi * y ) + case (6) ! cos(2 pi y/L) + nw2_tracer_dist = cos( 2.0 * pi * y ) + case (7) ! sin(2 pi y/L) + nw2_tracer_dist = sin( 2.0 * pi * y ) + case (8) ! cos(pi z/H) + nw2_tracer_dist = cos( pi * z ) + case (9) ! sin(pi z/H) + nw2_tracer_dist = sin( pi * z ) + case default + stop 'This should not happen. Died in nw2_tracer_dist()!' + end select + nw2_tracer_dist = nw2_tracer_dist * G%mask2dT(i,j) +end function nw2_tracer_dist + +!> Deallocate any memory associated with this tracer package +subroutine nw2_tracers_end(CS) + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracers. + + integer :: m + + if (associated(CS)) then + if (associated(CS%tr)) deallocate(CS%tr) + deallocate(CS) + endif +end subroutine nw2_tracers_end + +!> \namespace nw2_tracers +!! +!! TBD + +end module nw2_tracers From ca517073610271080c5be107d7b36c63b38295fd Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Sun, 3 Jan 2021 17:19:23 -0500 Subject: [PATCH 02/31] Fix FPE problems with NW2 tracers - I had a "*" in place of a "/" for calculating a non-dimensional quantity that led to OOB numbers. --- src/tracer/nw2_tracers.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 index e71387f9c7..6b2bdf2c0d 100644 --- a/src/tracer/nw2_tracers.F90 +++ b/src/tracer/nw2_tracers.F90 @@ -244,7 +244,7 @@ subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US endif do m=1,CS%ntr - dt_x_rate = dt * ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) ) + dt_x_rate = dt / ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) ) !$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate) do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k) From 77c67fcfb416e10b49eb14b2a2efe5be77032271 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Sun, 3 Jan 2021 17:33:44 -0500 Subject: [PATCH 03/31] Allow NW2 tracesr to be added mid-run - I'd previously made the NW2 tracers mandatory in the restart field which prohibits adding tracers on a restart. This commit allows a restart without tracers to be used. --- src/tracer/nw2_tracers.F90 | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 index 6b2bdf2c0d..b3d751d084 100644 --- a/src/tracer/nw2_tracers.F90 +++ b/src/tracer/nw2_tracers.F90 @@ -87,7 +87,7 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS call query_vardesc(tr_desc, name=var_name, caller="register_nw2_tracers") ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, & - registry_diags=.true., restart_CS=restart_CS, mandatory=.true.) + registry_diags=.true., restart_CS=restart_CS, mandatory=.false.) if (m<=10) then CS%restore_rate(m) = 1.0 / 6.0 ! 6 years else @@ -148,14 +148,13 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS) call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode") endif - if (.not.restart) then + ! Initialize only if this is not a restart or we are using a restart + ! in which the tracers were not present + if ((.not.restart) .or. & + (.not. query_initialized(CS%tr(:,:,:,m),'nw2_tracer',CS%restart_CSp))) then do m=1,CS%ntr do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec - ! if (G%mask2dT(i,j) < 0.5) then - ! CS%tr(i,j,k,m) = 0. - ! else CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k) - ! endif enddo ; enddo ; enddo enddo ! Tracer loop endif ! restart From b333739df4e3c95f0681078ae1db9a3e63da338c Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Sun, 3 Jan 2021 17:35:16 -0500 Subject: [PATCH 04/31] Cleanup NW2 registration - I noticed an unnecessary function call which is now removed. --- src/tracer/nw2_tracers.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 index b3d751d084..3199a5b271 100644 --- a/src/tracer/nw2_tracers.F90 +++ b/src/tracer/nw2_tracers.F90 @@ -9,7 +9,7 @@ module nw2_tracers use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_time_manager, only : time_type, time_type_to_real use MOM_tracer_registry, only : register_tracer, tracer_registry_type @@ -84,7 +84,6 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS ! This is needed to force the compiler not to do a copy in the registration ! calls. Curses on the designers and implementers of Fortran90. tr_ptr => CS%tr(:,:,:,m) - call query_vardesc(tr_desc, name=var_name, caller="register_nw2_tracers") ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, & registry_diags=.true., restart_CS=restart_CS, mandatory=.false.) From adb3554add3c787093e8aeb0f7327f2d0feaede9 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Sun, 3 Jan 2021 18:00:20 -0500 Subject: [PATCH 05/31] Fixed NW2 tracer restarts - An inverted do/if loop was using uninitialized variables and was overwriting restart data. --- src/tracer/nw2_tracers.F90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 index 3199a5b271..5385771445 100644 --- a/src/tracer/nw2_tracers.F90 +++ b/src/tracer/nw2_tracers.F90 @@ -118,6 +118,7 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS) ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights real :: rscl ! z* scaling factor + character(len=8) :: var_name ! The variable's name. integer :: i, j, k, m if (.not.associated(CS)) return @@ -147,16 +148,17 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS) call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode") endif - ! Initialize only if this is not a restart or we are using a restart - ! in which the tracers were not present - if ((.not.restart) .or. & - (.not. query_initialized(CS%tr(:,:,:,m),'nw2_tracer',CS%restart_CSp))) then - do m=1,CS%ntr + do m=1,CS%ntr + ! Initialize only if this is not a restart or we are using a restart + ! in which the tracers were not present + write(var_name(1:8),'(a6,i2.2)') 'tracer',m + if ((.not.restart) .or. & + (.not. query_initialized(CS%tr(:,:,:,m),var_name,CS%restart_CSp))) then do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k) enddo ; enddo ; enddo - enddo ! Tracer loop - endif ! restart + endif ! restart + enddo ! Tracer loop end subroutine initialize_nw2_tracers From d20ecf7b6f382f6e62a7aee940ca7bf79d5c6263 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 22 Feb 2021 12:56:55 -0500 Subject: [PATCH 06/31] Remove use of parameter for NTR_MAX in NW2 - The number of tracers was hard coded via an integer parameter which is unnecessarily restrictive. --- src/tracer/nw2_tracers.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 index 5385771445..5417373b24 100644 --- a/src/tracer/nw2_tracers.F90 +++ b/src/tracer/nw2_tracers.F90 @@ -27,15 +27,13 @@ module nw2_tracers public nw2_tracer_column_physics public nw2_tracers_end -integer, parameter :: NTR_MAX = 20 !< the maximum number of tracers in this module. - !> The control structure for the nw2_tracers package type, public :: nw2_tracers_CS ; private - integer :: ntr = NTR_MAX !< The number of tracers that are actually used. + integer :: ntr = 0 !< The number of tracers that are actually used. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this package, in g m-3? - real, dimension(NTR_MAX) :: restore_rate !< The exponential growth rate for restoration value [year-1]. + real, allocatable , dimension(:) :: restore_rate !< The exponential growth rate for restoration value [year-1]. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart controls structure @@ -57,7 +55,7 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS ! This include declares and sets the variable "version". #include "version_variable.h" - character(len=40) :: mdl = "nw2_example" ! This module's name. + character(len=40) :: mdl = "nw2_tracers" ! This module's name. character(len=200) :: inputdir ! The directory where the input files are. character(len=8) :: var_name ! The variable's name. real, pointer :: tr_ptr(:,:,:) => NULL() @@ -76,7 +74,9 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") + CS%ntr=20 allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0 + allocate(CS%restore_rate(CS%ntr)) do m=1,CS%ntr write(var_name(1:8),'(a6,i2.2)') 'tracer',m From 8f6ec58bc0246e923a5a6430a2524a5bb693fbdc Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 22 Feb 2021 15:08:14 -0500 Subject: [PATCH 07/31] Replace NW2 2x10 tracers with 3x3 - The first interation of NW2 tracers implemented 20 tracesr with 10 different spatial structures and 2 different restoration time scales. This has been replaced with 9 tracers with 3 spatial, corresponding to sin(2x), y, z, and 3 time scales. - Number of tracer groups and time scales are now run-time configurable. --- src/tracer/nw2_tracers.F90 | 51 +++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 28 deletions(-) diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 index 5417373b24..e75c5c5d38 100644 --- a/src/tracer/nw2_tracers.F90 +++ b/src/tracer/nw2_tracers.F90 @@ -60,7 +60,9 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS character(len=8) :: var_name ! The variable's name. real, pointer :: tr_ptr(:,:,:) => NULL() logical :: do_nw2 - integer :: isd, ied, jsd, jed, nz, m + integer :: isd, ied, jsd, jed, nz, m, ig + integer :: n_groups ! Number of groups of three tracers (i.e. # tracers/3) + real, allocatable, dimension(:) :: timescale_in_days type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke @@ -73,8 +75,18 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") - - CS%ntr=20 + call get_param(param_file, mdl, "NW2_TRACER_GROUPS", n_groups, & + "The number of tracer groups where a group is of three tracers "//& + "initialized and restored to sin(x), y and z, respectively. Each "//& + "group is restored with an independent restoration rate.", & + default=3) + allocate(timescale_in_days(n_groups)) + timescale_in_days = (/365., 730., 1460./) + call get_param(param_file, mdl, "NW2_TRACER_RESTORE_TIMESCALE", timescale_in_days, & + "A list of timescales, one for each tracer group.", & + units="days") + + CS%ntr = 3 * n_groups allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0 allocate(CS%restore_rate(CS%ntr)) @@ -87,11 +99,8 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, & registry_diags=.true., restart_CS=restart_CS, mandatory=.false.) - if (m<=10) then - CS%restore_rate(m) = 1.0 / 6.0 ! 6 years - else - CS%restore_rate(m) = 1.0 ! 1 year - endif + ig = int( (m+2)/3 ) ! maps (1,2,3)->1, (4,5,6)->2, ... + CS%restore_rate(m) = 1.0 / ( timescale_in_days(ig) * 86400.0 ) enddo CS%tr_Reg => tr_Reg @@ -244,7 +253,7 @@ subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US endif do m=1,CS%ntr - dt_x_rate = dt / ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) ) + dt_x_rate = ( dt * CS%restore_rate(m) ) * US%T_to_s !$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate) do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k) @@ -272,27 +281,13 @@ real function nw2_tracer_dist(m, G, GV, eta, i, j, k) x = ( G%geolonT(i,j) - G%west_lon ) / G%len_lon ! 0 ... 1 y = -G%geolatT(i,j) / G%south_lat ! -1 ... 1 z = - 0.5 * ( eta(i,j,K-1) + eta(i,j,K) ) / GV%max_depth ! 0 ... 1 - select case ( mod(m-1,10) ) - case (0) ! y/L + select case ( mod(m-1,3) ) + case (0) ! sin(2 pi x/L) + nw2_tracer_dist = sin( 2.0 * pi * x ) + case (1) ! y/L nw2_tracer_dist = y - case (1) ! -z/L + case (2) ! -z/L nw2_tracer_dist = -z - case (2) ! cos(2 pi x/L) - nw2_tracer_dist = cos( 2.0 * pi * x ) - case (3) ! sin(2 pi x/L) - nw2_tracer_dist = sin( 2.0 * pi * x ) - case (4) ! sin(4 pi x/L) - nw2_tracer_dist = sin( 2.0 * pi * mod( 2.0*x, 1.0 ) ) - case (5) ! sin(pi y/L) - nw2_tracer_dist = sin( pi * y ) - case (6) ! cos(2 pi y/L) - nw2_tracer_dist = cos( 2.0 * pi * y ) - case (7) ! sin(2 pi y/L) - nw2_tracer_dist = sin( 2.0 * pi * y ) - case (8) ! cos(pi z/H) - nw2_tracer_dist = cos( pi * z ) - case (9) ! sin(pi z/H) - nw2_tracer_dist = sin( pi * z ) case default stop 'This should not happen. Died in nw2_tracer_dist()!' end select From ffbd82a809374e6df5b7b8b4a7d085b574c7b07d Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 29 Apr 2021 10:33:52 -0400 Subject: [PATCH 08/31] Performance monitoring; GitHub integration This patch introduces two new testing targets to the verification suite based on a small configuration based on the `benchmark` regression test. The profile test is saved a `p0` in `.testing`. Future tests can be included if appropriate. The new targets: * `make profile`: Run the model and record the FMS timings. * `make perf`: Run the model through the `perf` tool and record timings for the resolvable functions (as symbols). In both cases, the timings are converted to JSON output files and the top results are reported to stdout, and readable in GitHub actions output. It can also be run locally. Support Python scripts have been included to do this work. This will require a functional Python environment. Some system and configuration data is logged alongside the timings, but this is still rather incomplete and needs some further planning. Times are compared to the target build (usually dev/gfdl). ANSI terminal highlighting (i.e. color) is to used to highlight excessive differences. Current issues: - Model configuration - GitHub timings are still rather unreliable, and should currently only be treated as crude estimates. This should be considered a work in progress. - The GitHub profiling rule still builds the standard configuration, evem though it is unused. - Additional tools are required to push the timings to some database, either a local sqlite3 or an external one. --- .github/actions/testing-setup/action.yml | 1 + .github/workflows/perfmon.yml | 36 ++ .testing/Makefile | 76 +++- .testing/p0/MOM_input | 505 +++++++++++++++++++++++ .testing/p0/MOM_override | 0 .testing/p0/diag_table | 91 ++++ .testing/p0/input.nml | 22 + .testing/tools/compare_clocks.py | 88 ++++ .testing/tools/compare_perf.py | 110 +++++ .testing/tools/parse_fms_clocks.py | 120 ++++++ .testing/tools/parse_perf.py | 71 ++++ 11 files changed, 1119 insertions(+), 1 deletion(-) create mode 100644 .github/workflows/perfmon.yml create mode 100644 .testing/p0/MOM_input create mode 100644 .testing/p0/MOM_override create mode 100644 .testing/p0/diag_table create mode 100644 .testing/p0/input.nml create mode 100755 .testing/tools/compare_clocks.py create mode 100755 .testing/tools/compare_perf.py create mode 100755 .testing/tools/parse_fms_clocks.py create mode 100755 .testing/tools/parse_perf.py diff --git a/.github/actions/testing-setup/action.yml b/.github/actions/testing-setup/action.yml index c47270af0d..1ab96aa3df 100644 --- a/.github/actions/testing-setup/action.yml +++ b/.github/actions/testing-setup/action.yml @@ -33,6 +33,7 @@ runs: echo "::group::Install linux packages" sudo apt-get update sudo apt-get install netcdf-bin libnetcdf-dev libnetcdff-dev mpich libmpich-dev + sudo apt-get install linux-tools-common echo "::endgroup::" - name: Compile FMS library diff --git a/.github/workflows/perfmon.yml b/.github/workflows/perfmon.yml new file mode 100644 index 0000000000..00e645c4fd --- /dev/null +++ b/.github/workflows/perfmon.yml @@ -0,0 +1,36 @@ +name: Performance Monitor + +on: [pull_request] + +jobs: + build-test-perfmon: + + runs-on: ubuntu-latest + defaults: + run: + working-directory: .testing + + steps: + - uses: actions/checkout@v2 + with: + submodules: recursive + + - uses: ./.github/actions/testing-setup + + - name: Compile optimized models + run: >- + make -j build.prof + MOM_TARGET_SLUG=$GITHUB_REPOSITORY + MOM_TARGET_LOCAL_BRANCH=$GITHUB_BASE_REF + DO_REGRESSION_TESTS=true + + - name: Generate profile data + run: >- + pip install f90nml && + make profile + DO_REGRESSION_TESTS=true + + - name: Generate perf data + run: | + sudo sysctl -w kernel.perf_event_paranoid=2 + make perf DO_REGRESSION_TESTS=true diff --git a/.testing/Makefile b/.testing/Makefile index 45d05cd23f..89ea64e1e2 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -34,6 +34,7 @@ # Build configuration: # FCFLAGS_DEBUG Testing ("debug") compiler flags # FCFLAGS_REPRO Production ("repro") compiler flags +# FCFLAGS_OPT Aggressive optimization compiler flags # FCFLAGS_INIT Variable initialization flags # FCFLAGS_COVERAGE Code coverage flags # @@ -72,6 +73,7 @@ export MPIFC # NOTE: FMS will be built using FCFLAGS_DEBUG FCFLAGS_DEBUG ?= -g -O0 FCFLAGS_REPRO ?= -g -O2 +FCFLAGS_OPT ?= -g -O3 -mavx -fno-omit-frame-pointer FCFLAGS_INIT ?= FCFLAGS_COVERAGE ?= # Additional notes: @@ -96,6 +98,7 @@ DO_REPRO_TESTS ?= # Time measurement (configurable by the CI) TIME ?= time + #--- # Dependencies DEPS = deps @@ -122,6 +125,11 @@ ifeq ($(DO_REPRO_TESTS), true) TESTS += repros endif +# Profiling +ifeq ($(DO_PROFILE), false) + BUILDS += opt opt_target +endif + # The following variables are configured by Travis: # DO_REGRESSION_TESTS: true if $(TRAVIS_PULL_REQUEST) is a PR number # MOM_TARGET_SLUG: TRAVIS_REPO_SLUG @@ -195,6 +203,7 @@ endif .PHONY: all build.regressions all: $(foreach b,$(BUILDS),build/$(b)/MOM6) $(VENV_PATH) build.regressions: $(foreach b,symmetric target,build/$(b)/MOM6) +build.prof: $(foreach b,opt opt_target,build/$(b)/MOM6) # Executable BUILD_TARGETS = MOM6 Makefile path_names @@ -217,6 +226,7 @@ PATH_FMS = PATH="${PATH}:../../$(DEPS)/bin" SYMMETRIC_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(COVERAGE) $(FCFLAGS_FMS)" ASYMMETRIC_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" REPRO_FCFLAGS := FCFLAGS="$(FCFLAGS_REPRO) $(FCFLAGS_FMS)" +OPT_FCFLAGS := FCFLAGS="$(FCFLAGS_OPT) $(FCFLAGS_FMS)" OPENMP_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" TARGET_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" @@ -230,6 +240,8 @@ build/asymmetric/Makefile: MOM_ENV=$(PATH_FMS) $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLA build/repro/Makefile: MOM_ENV=$(PATH_FMS) $(REPRO_FCFLAGS) $(MOM_LDFLAGS) build/openmp/Makefile: MOM_ENV=$(PATH_FMS) $(OPENMP_FCFLAGS) $(MOM_LDFLAGS) build/target/Makefile: MOM_ENV=$(PATH_FMS) $(TARGET_FCFLAGS) $(MOM_LDFLAGS) +build/opt/Makefile: MOM_ENV=$(PATH_FMS) $(OPT_FCFLAGS) $(MOM_LDFLAGS) +build/opt_target/Makefile: MOM_ENV=$(PATH_FMS) $(OPT_FCFLAGS) $(MOM_LDFLAGS) build/coupled/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) build/nuopc/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) build/mct/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) @@ -240,12 +252,15 @@ build/asymmetric/Makefile: MOM_ACFLAGS=--enable-asymmetric build/repro/Makefile: MOM_ACFLAGS= build/openmp/Makefile: MOM_ACFLAGS=--enable-openmp build/target/Makefile: MOM_ACFLAGS= +build/opt/Makefile: MOM_ACFLAGS= +build/opt_target/Makefile: MOM_ACFLAGS= build/coupled/Makefile: MOM_ACFLAGS=--with-driver=coupled_driver build/nuopc/Makefile: MOM_ACFLAGS=--with-driver=nuopc_driver build/mct/Makefile: MOM_ACFLAGS=--with-driver=mct_driver # Fetch regression target source code build/target/Makefile: | $(TARGET_CODEBASE) +build/opt_target/Makefile: | $(TARGET_CODEBASE) # Define source code dependencies @@ -276,7 +291,8 @@ build/%/Makefile: ../ac/configure ../ac/Makefile.in $(DEPS)/lib/libFMS.a $(MKMF) # Fetch the regression target codebase -build/target/Makefile: $(TARGET_CODEBASE)/ac/configure $(DEPS)/lib/libFMS.a $(MKMF) $(LIST_PATHS) +build/target/Makefile build/opt_target/Makefile: \ + $(TARGET_CODEBASE)/ac/configure $(DEPS)/lib/libFMS.a $(MKMF) $(LIST_PATHS) mkdir -p $(@D) cd $(@D) \ && $(MOM_ENV) ../../$(TARGET_CODEBASE)/ac/configure $(MOM_ACFLAGS) \ @@ -620,6 +636,64 @@ test.summary: fi +#--- +# Profiling +# XXX: This is experimental work to track, log, and report changes in runtime +PCONFIGS = p0 + +.PHONY: profile +profile: $(foreach p,$(PCONFIGS), prof.$(p)) + +.PHONY: prof.p0 +prof.p0: work/p0/opt/clocks.json work/p0/opt_target/clocks.json + python tools/compare_clocks.py $^ + +work/p0/%/clocks.json: work/p0/%/std.out + python tools/parse_fms_clocks.py -d $(@D) $^ > $@ + +work/p0/opt/std.out: build/opt/MOM6 +work/p0/opt_target/std.out: build/opt_target/MOM6 + +work/p0/%/std.out: + mkdir -p $(@D) + cp -RL p0/* $(@D) + mkdir -p $(@D)/RESTART + echo -e "" > $(@D)/MOM_override + cd $(@D) \ + && $(MPIRUN) -n 1 ../../../$< 2> std.err > std.out + +#--- +# Same but with perf + +# TODO: This expects the -e flag, can I handle it in the command? +PERF_EVENTS ?= + +.PHONY: perf +perf: $(foreach p,$(PCONFIGS), perf.$(p)) + +.PHONY: prof.p0 +perf.p0: work/p0/opt/profile.json work/p0/opt_target/profile.json + python tools/compare_perf.py $^ + +work/p0/%/profile.json: work/p0/%/perf.data + python tools/parse_perf.py -f $< > $@ + +work/p0/opt/perf.data: build/opt/MOM6 +work/p0/opt_target/perf.data: build/opt_target/MOM6 + +work/p0/%/perf.data: + mkdir -p $(@D) + cp -RL p0/* $(@D) + mkdir -p $(@D)/RESTART + echo -e "" > $(@D)/MOM_override + cd $(@D) \ + && perf record \ + -F 3999 \ + ${PERF_EVENTS} \ + ../../../$< 2> std.perf.err > std.perf.out \ + || cat std.perf.err + + #---- # NOTE: These tests assert that we are in the .testing directory. diff --git a/.testing/p0/MOM_input b/.testing/p0/MOM_input new file mode 100644 index 0000000000..8f751d7bf1 --- /dev/null +++ b/.testing/p0/MOM_input @@ -0,0 +1,505 @@ +! This input file provides the adjustable run-time parameters for version 6 of the Modular Ocean Model (MOM6). +! Where appropriate, parameters use usually given in MKS units. + +! This particular file is for the example in benchmark. + +! This MOM_input file typically contains only the non-default values that are needed to reproduce this example. +! A full list of parameters for this example can be found in the corresponding MOM_parameter_doc.all file +! which is generated by the model at run-time. + +! === module MOM_domains === +NIGLOBAL = 32 ! + ! The total number of thickness grid points in the x-direction in the physical + ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. +NJGLOBAL = 32 ! + ! The total number of thickness grid points in the y-direction in the physical + ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. +LAYOUT = 1, 1 + ! The processor layout that was actually used. + +! === module MOM === +THICKNESSDIFFUSE = True ! [Boolean] default = False + ! If true, interface heights are diffused with a coefficient of KHTH. +THICKNESSDIFFUSE_FIRST = True ! [Boolean] default = False + ! If true, do thickness diffusion before dynamics. This is only used if + ! THICKNESSDIFFUSE is true. +DT = 900.0 ! [s] + ! The (baroclinic) dynamics time step. The time-step that is actually used will + ! be an integer fraction of the forcing time-step (DT_FORCING in ocean-only mode + ! or the coupling timestep in coupled mode.) +DT_THERM = 3600.0 ! [s] default = 900.0 + ! The thermodynamic and tracer advection time step. Ideally DT_THERM should be + ! an integer multiple of DT and less than the forcing or coupling time-step, + ! unless THERMO_SPANS_COUPLING is true, in which case DT_THERM can be an integer + ! multiple of the coupling timestep. By default DT_THERM is set to DT. +DTBT_RESET_PERIOD = 0.0 ! [s] default = 3600.0 + ! The period between recalculations of DTBT (if DTBT <= 0). If DTBT_RESET_PERIOD + ! is negative, DTBT is set based only on information available at + ! initialization. If 0, DTBT will be set every dynamics time step. The default + ! is set by DT_THERM. This is only used if SPLIT is true. +FRAZIL = True ! [Boolean] default = False + ! If true, water freezes if it gets too cold, and the accumulated heat deficit + ! is returned in the surface state. FRAZIL is only used if + ! ENABLE_THERMODYNAMICS is true. +C_P = 3925.0 ! [J kg-1 K-1] default = 3991.86795711963 + ! The heat capacity of sea water, approximated as a constant. This is only used + ! if ENABLE_THERMODYNAMICS is true. The default value is from the TEOS-10 + ! definition of conservative temperature. +SAVE_INITIAL_CONDS = True ! [Boolean] default = False + ! If true, write the initial conditions to a file given by IC_OUTPUT_FILE. +IC_OUTPUT_FILE = "GOLD_IC" ! default = "MOM_IC" + ! The file into which to write the initial conditions. + +! === module MOM_fixed_initialization === +INPUTDIR = "INPUT" ! default = "." + ! The directory in which input files are found. + +! === module MOM_grid_init === +GRID_CONFIG = "cartesian" ! + ! A character string that determines the method for defining the horizontal + ! grid. Current options are: + ! mosaic - read the grid from a mosaic (supergrid) + ! file set by GRID_FILE. + ! cartesian - use a (flat) Cartesian grid. + ! spherical - use a simple spherical grid. + ! mercator - use a Mercator spherical grid. +SOUTHLAT = -41.0 ! [degrees] + ! The southern latitude of the domain. +LENLAT = 41.0 ! [degrees] + ! The latitudinal length of the domain. +LENLON = 90.0 ! [degrees] + ! The longitudinal length of the domain. +ISOTROPIC = True ! [Boolean] default = False + ! If true, an isotropic grid on a sphere (also known as a Mercator grid) is + ! used. With an isotropic grid, the meridional extent of the domain (LENLAT), + ! the zonal extent (LENLON), and the number of grid points in each direction are + ! _not_ independent. In MOM the meridional extent is determined to fit the zonal + ! extent and the number of grid points, while grid is perfectly isotropic. +TOPO_CONFIG = "benchmark" ! + ! This specifies how bathymetry is specified: + ! file - read bathymetric information from the file + ! specified by (TOPO_FILE). + ! flat - flat bottom set to MAXIMUM_DEPTH. + ! bowl - an analytically specified bowl-shaped basin + ! ranging between MAXIMUM_DEPTH and MINIMUM_DEPTH. + ! spoon - a similar shape to 'bowl', but with an vertical + ! wall at the southern face. + ! halfpipe - a zonally uniform channel with a half-sine + ! profile in the meridional direction. + ! bbuilder - build topography from list of functions. + ! benchmark - use the benchmark test case topography. + ! Neverworld - use the Neverworld test case topography. + ! DOME - use a slope and channel configuration for the + ! DOME sill-overflow test case. + ! ISOMIP - use a slope and channel configuration for the + ! ISOMIP test case. + ! DOME2D - use a shelf and slope configuration for the + ! DOME2D gravity current/overflow test case. + ! Kelvin - flat but with rotated land mask. + ! seamount - Gaussian bump for spontaneous motion test case. + ! dumbbell - Sloshing channel with reservoirs on both ends. + ! shelfwave - exponential slope for shelfwave test case. + ! Phillips - ACC-like idealized topography used in the Phillips config. + ! dense - Denmark Strait-like dense water formation and overflow. + ! USER - call a user modified routine. + +! === module benchmark_initialize_topography === +MINIMUM_DEPTH = 1.0 ! [m] default = 0.0 + ! The minimum depth of the ocean. +MAXIMUM_DEPTH = 5500.0 ! [m] + ! The maximum depth of the ocean. + +! === module MOM_verticalGrid === +! Parameters providing information about the vertical grid. +NK = 75 ! [nondim] + ! The number of model layers. + +! === module MOM_EOS === + +! === module MOM_tracer_flow_control === +USE_IDEAL_AGE_TRACER = True ! [Boolean] default = False + ! If true, use the ideal_age_example tracer package. + +! === module ideal_age_example === + +! === module MOM_coord_initialization === +COORD_CONFIG = "ts_range" ! default = "none" + ! This specifies how layers are to be defined: + ! ALE or none - used to avoid defining layers in ALE mode + ! file - read coordinate information from the file + ! specified by (COORD_FILE). + ! BFB - Custom coords for buoyancy-forced basin case + ! based on SST_S, T_BOT and DRHO_DT. + ! linear - linear based on interfaces not layers + ! layer_ref - linear based on layer densities + ! ts_ref - use reference temperature and salinity + ! ts_range - use range of temperature and salinity + ! (T_REF and S_REF) to determine surface density + ! and GINT calculate internal densities. + ! gprime - use reference density (RHO_0) for surface + ! density and GINT calculate internal densities. + ! ts_profile - use temperature and salinity profiles + ! (read from COORD_FILE) to set layer densities. + ! USER - call a user modified routine. +TS_RANGE_T_LIGHT = 25.0 ! [degC] default = 10.0 + ! The initial temperature of the lightest layer when COORD_CONFIG is set to + ! ts_range. +TS_RANGE_T_DENSE = 3.0 ! [degC] default = 10.0 + ! The initial temperature of the densest layer when COORD_CONFIG is set to + ! ts_range. +TS_RANGE_RESOLN_RATIO = 5.0 ! [nondim] default = 1.0 + ! The ratio of density space resolution in the densest part of the range to that + ! in the lightest part of the range when COORD_CONFIG is set to ts_range. Values + ! greater than 1 increase the resolution of the denser water. + +! === module MOM_state_initialization === +THICKNESS_CONFIG = "benchmark" ! default = "uniform" + ! A string that determines how the initial layer thicknesses are specified for a + ! new run: + ! file - read interface heights from the file specified + ! thickness_file - read thicknesses from the file specified + ! by (THICKNESS_FILE). + ! coord - determined by ALE coordinate. + ! uniform - uniform thickness layers evenly distributed + ! between the surface and MAXIMUM_DEPTH. + ! list - read a list of positive interface depths. + ! DOME - use a slope and channel configuration for the + ! DOME sill-overflow test case. + ! ISOMIP - use a configuration for the + ! ISOMIP test case. + ! benchmark - use the benchmark test case thicknesses. + ! Neverworld - use the Neverworld test case thicknesses. + ! search - search a density profile for the interface + ! densities. This is not yet implemented. + ! circle_obcs - the circle_obcs test case is used. + ! DOME2D - 2D version of DOME initialization. + ! adjustment2d - 2D lock exchange thickness ICs. + ! sloshing - sloshing gravity thickness ICs. + ! seamount - no motion test with seamount ICs. + ! dumbbell - sloshing channel ICs. + ! soliton - Equatorial Rossby soliton. + ! rossby_front - a mixed layer front in thermal wind balance. + ! USER - call a user modified routine. + +! === module benchmark_initialize_thickness === +TS_CONFIG = "benchmark" ! + ! A string that determines how the initial tempertures and salinities are + ! specified for a new run: + ! file - read velocities from the file specified + ! by (TS_FILE). + ! fit - find the temperatures that are consistent with + ! the layer densities and salinity S_REF. + ! TS_profile - use temperature and salinity profiles + ! (read from TS_FILE) to set layer densities. + ! benchmark - use the benchmark test case T & S. + ! linear - linear in logical layer space. + ! DOME2D - 2D DOME initialization. + ! ISOMIP - ISOMIP initialization. + ! adjustment2d - 2d lock exchange T/S ICs. + ! sloshing - sloshing mode T/S ICs. + ! seamount - no motion test with seamount ICs. + ! dumbbell - sloshing channel ICs. + ! rossby_front - a mixed layer front in thermal wind balance. + ! SCM_CVMix_tests - used in the SCM CVMix tests. + ! USER - call a user modified routine. + +! === module MOM_diag_mediator === + +! === module MOM_lateral_mixing_coeffs === +USE_VARIABLE_MIXING = True ! [Boolean] default = False + ! If true, the variable mixing code will be called. This allows diagnostics to + ! be created even if the scheme is not used. If KHTR_SLOPE_CFF>0 or + ! KhTh_Slope_Cff>0, this is set to true regardless of what is in the parameter + ! file. +USE_VISBECK = True ! [Boolean] default = False + ! If true, use the Visbeck et al. (1997) formulation for + ! thickness diffusivity. +RESOLN_SCALED_KH = True ! [Boolean] default = False + ! If true, the Laplacian lateral viscosity is scaled away when the first + ! baroclinic deformation radius is well resolved. +RESOLN_SCALED_KHTH = True ! [Boolean] default = False + ! If true, the interface depth diffusivity is scaled away when the first + ! baroclinic deformation radius is well resolved. +RESOLN_SCALED_KHTR = True ! [Boolean] default = False + ! If true, the epipycnal tracer diffusivity is scaled away when the first + ! baroclinic deformation radius is well resolved. +KHTH_SLOPE_CFF = 0.1 ! [nondim] default = 0.0 + ! The nondimensional coefficient in the Visbeck formula for the interface depth + ! diffusivity +KHTR_SLOPE_CFF = 0.1 ! [nondim] default = 0.0 + ! The nondimensional coefficient in the Visbeck formula for the epipycnal tracer + ! diffusivity +VARMIX_KTOP = 6 ! [nondim] default = 2 + ! The layer number at which to start vertical integration of S*N for purposes of + ! finding the Eady growth rate. +VISBECK_L_SCALE = 3.0E+04 ! [m] default = 0.0 + ! The fixed length scale in the Visbeck formula. + +! === module MOM_set_visc === +PRANDTL_TURB = 0.0 ! [nondim] default = 1.0 + ! The turbulent Prandtl number applied to shear instability. +DYNAMIC_VISCOUS_ML = True ! [Boolean] default = False + ! If true, use a bulk Richardson number criterion to determine the mixed layer + ! thickness for viscosity. +ML_OMEGA_FRAC = 1.0 ! [nondim] default = 0.0 + ! When setting the decay scale for turbulence, use this fraction of the absolute + ! rotation rate blended with the local value of f, as sqrt((1-of)*f^2 + + ! of*4*omega^2). +HBBL = 10.0 ! [m] + ! The thickness of a bottom boundary layer with a viscosity of KVBBL if + ! BOTTOMDRAGLAW is not defined, or the thickness over which near-bottom + ! velocities are averaged for the drag law if BOTTOMDRAGLAW is defined but + ! LINEAR_DRAG is not. +DRAG_BG_VEL = 0.1 ! [m s-1] default = 0.0 + ! DRAG_BG_VEL is either the assumed bottom velocity (with LINEAR_DRAG) or an + ! unresolved velocity that is combined with the resolved velocity to estimate + ! the velocity magnitude. DRAG_BG_VEL is only used when BOTTOMDRAGLAW is + ! defined. +BBL_THICK_MIN = 0.1 ! [m] default = 0.0 + ! The minimum bottom boundary layer thickness that can be used with + ! BOTTOMDRAGLAW. This might be Kv/(cdrag*drag_bg_vel) to give Kv as the minimum + ! near-bottom viscosity. +KV = 1.0E-04 ! [m2 s-1] + ! The background kinematic viscosity in the interior. The molecular value, ~1e-6 + ! m2 s-1, may be used. + +! === module MOM_thickness_diffuse === +KHTH = 1.0 ! [m2 s-1] default = 0.0 + ! The background horizontal thickness diffusivity. +KHTH_MAX = 900.0 ! [m2 s-1] default = 0.0 + ! The maximum horizontal thickness diffusivity. + +! === module MOM_dynamics_split_RK2 === + +! === module MOM_continuity === + +! === module MOM_continuity_PPM === +ETA_TOLERANCE = 1.0E-06 ! [m] default = 1.1E-09 + ! The tolerance for the differences between the barotropic and baroclinic + ! estimates of the sea surface height due to the fluxes through each face. The + ! total tolerance for SSH is 4 times this value. The default is + ! 0.5*NK*ANGSTROM, and this should not be set less than about + ! 10^-15*MAXIMUM_DEPTH. +VELOCITY_TOLERANCE = 0.001 ! [m s-1] default = 3.0E+08 + ! The tolerance for barotropic velocity discrepancies between the barotropic + ! solution and the sum of the layer thicknesses. + +! === module MOM_CoriolisAdv === +BOUND_CORIOLIS = True ! [Boolean] default = False + ! If true, the Coriolis terms at u-points are bounded by the four estimates of + ! (f+rv)v from the four neighboring v-points, and similarly at v-points. This + ! option would have no effect on the SADOURNY Coriolis scheme if it were + ! possible to use centered difference thickness fluxes. + +! === module MOM_PressureForce === + +! === module MOM_PressureForce_AFV === + +! === module MOM_hor_visc === +AH_VEL_SCALE = 0.05 ! [m s-1] default = 0.0 + ! The velocity scale which is multiplied by the cube of the grid spacing to + ! calculate the biharmonic viscosity. The final viscosity is the largest of this + ! scaled viscosity, the Smagorinsky and Leith viscosities, and AH. +SMAGORINSKY_AH = True ! [Boolean] default = False + ! If true, use a biharmonic Smagorinsky nonlinear eddy viscosity. +SMAG_BI_CONST = 0.06 ! [nondim] default = 0.0 + ! The nondimensional biharmonic Smagorinsky constant, typically 0.015 - 0.06. + +! === module MOM_vert_friction === +MAXVEL = 10.0 ! [m s-1] default = 3.0E+08 + ! The maximum velocity allowed before the velocity components are truncated. + +! === module MOM_barotropic === +BOUND_BT_CORRECTION = True ! [Boolean] default = False + ! If true, the corrective pseudo mass-fluxes into the barotropic solver are + ! limited to values that require less than maxCFL_BT_cont to be accommodated. +NONLINEAR_BT_CONTINUITY = True ! [Boolean] default = False + ! If true, use nonlinear transports in the barotropic continuity equation. This + ! does not apply if USE_BT_CONT_TYPE is true. +BT_PROJECT_VELOCITY = True ! [Boolean] default = False + ! If true, step the barotropic velocity first and project out the velocity + ! tendency by 1+BEBT when calculating the transport. The default (false) is to + ! use a predictor continuity step to find the pressure field, and then to do a + ! corrector continuity step using a weighted average of the old and new + ! velocities, with weights of (1-BEBT) and BEBT. +BEBT = 0.2 ! [nondim] default = 0.1 + ! BEBT determines whether the barotropic time stepping uses the forward-backward + ! time-stepping scheme or a backward Euler scheme. BEBT is valid in the range + ! from 0 (for a forward-backward treatment of nonrotating gravity waves) to 1 + ! (for a backward Euler treatment). In practice, BEBT must be greater than about + ! 0.05. +DTBT = -0.95 ! [s or nondim] default = -0.98 + ! The barotropic time step, in s. DTBT is only used with the split explicit time + ! stepping. To set the time step automatically based the maximum stable value + ! use 0, or a negative value gives the fraction of the stable value. Setting + ! DTBT to 0 is the same as setting it to -0.98. The value of DTBT that will + ! actually be used is an integer fraction of DT, rounding down. + +! === module MOM_mixed_layer_restrat === +MIXEDLAYER_RESTRAT = True ! [Boolean] default = False + ! If true, a density-gradient dependent re-stratifying flow is imposed in the + ! mixed layer. Can be used in ALE mode without restriction but in layer mode can + ! only be used if BULKMIXEDLAYER is true. +FOX_KEMPER_ML_RESTRAT_COEF = 5.0 ! [nondim] default = 0.0 + ! A nondimensional coefficient that is proportional to the ratio of the + ! deformation radius to the dominant lengthscale of the submesoscale mixed layer + ! instabilities, times the minimum of the ratio of the mesoscale eddy kinetic + ! energy to the large-scale geostrophic kinetic energy or 1 plus the square of + ! the grid spacing over the deformation radius, as detailed by Fox-Kemper et al. + ! (2010) + +! === module MOM_diagnostics === + +! === module MOM_diabatic_driver === +! The following parameters are used for diabatic processes. + +! === module MOM_entrain_diffusive === +MAX_ENT_IT = 20 ! default = 5 + ! The maximum number of iterations that may be used to calculate the interior + ! diapycnal entrainment. +TOLERANCE_ENT = 1.0E-05 ! [m] default = 1.341640786499874E-05 + ! The tolerance with which to solve for entrainment values. + +! === module MOM_set_diffusivity === + +! === module MOM_bkgnd_mixing === +! Adding static vertical background mixing coefficients +KD = 2.0E-05 ! [m2 s-1] default = 0.0 + ! The background diapycnal diffusivity of density in the interior. Zero or the + ! molecular value, ~1e-7 m2 s-1, may be used. + +! === module MOM_kappa_shear === +! Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008 +USE_JACKSON_PARAM = True ! [Boolean] default = False + ! If true, use the Jackson-Hallberg-Legg (JPO 2008) shear mixing + ! parameterization. +MAX_RINO_IT = 25 ! [nondim] default = 50 + ! The maximum number of iterations that may be used to estimate the Richardson + ! number driven mixing. + +! === module MOM_diabatic_aux === +! The following parameters are used for auxiliary diabatic processes. +RECLAIM_FRAZIL = False ! [Boolean] default = True + ! If true, try to use any frazil heat deficit to cool any overlying layers down + ! to the freezing point, thereby avoiding the creation of thin ice when the SST + ! is above the freezing point. + +! === module MOM_mixed_layer === +MSTAR = 0.3 ! [nondim] default = 1.2 + ! The ratio of the friction velocity cubed to the TKE input to the mixed layer. +BULK_RI_ML = 0.05 ! [nondim] + ! The efficiency with which mean kinetic energy released by mechanically forced + ! entrainment of the mixed layer is converted to turbulent kinetic energy. +ABSORB_ALL_SW = True ! [Boolean] default = False + ! If true, all shortwave radiation is absorbed by the ocean, instead of passing + ! through to the bottom mud. +TKE_DECAY = 10.0 ! [nondim] default = 2.5 + ! TKE_DECAY relates the vertical rate of decay of the TKE available for + ! mechanical entrainment to the natural Ekman depth. +HMIX_MIN = 2.0 ! [m] default = 0.0 + ! The minimum mixed layer depth if the mixed layer depth is determined + ! dynamically. +LIMIT_BUFFER_DETRAIN = True ! [Boolean] default = False + ! If true, limit the detrainment from the buffer layers to not be too different + ! from the neighbors. +DEPTH_LIMIT_FLUXES = 0.1 ! [m] default = 0.2 + ! The surface fluxes are scaled away when the total ocean depth is less than + ! DEPTH_LIMIT_FLUXES. +CORRECT_ABSORPTION_DEPTH = True ! [Boolean] default = False + ! If true, the average depth at which penetrating shortwave radiation is + ! absorbed is adjusted to match the average heating depth of an exponential + ! profile by moving some of the heating upward in the water column. + +! === module MOM_opacity === +PEN_SW_SCALE = 15.0 ! [m] default = 0.0 + ! The vertical absorption e-folding depth of the penetrating shortwave + ! radiation. +PEN_SW_FRAC = 0.42 ! [nondim] default = 0.0 + ! The fraction of the shortwave radiation that penetrates below the surface. + +! === module MOM_tracer_advect === + +! === module MOM_tracer_hor_diff === +KHTR = 1.0 ! [m2 s-1] default = 0.0 + ! The background along-isopycnal tracer diffusivity. +KHTR_MAX = 900.0 ! [m2 s-1] default = 0.0 + ! The maximum along-isopycnal tracer diffusivity. +DIFFUSE_ML_TO_INTERIOR = True ! [Boolean] default = False + ! If true, enable epipycnal mixing between the surface boundary layer and the + ! interior. +ML_KHTR_SCALE = 0.0 ! [nondim] default = 1.0 + ! With Diffuse_ML_interior, the ratio of the truly horizontal diffusivity in the + ! mixed layer to the epipycnal diffusivity. The valid range is 0 to 1. + +! === module MOM_sum_output === +MAXTRUNC = 5000 ! [truncations save_interval-1] default = 0 + ! The run will be stopped, and the day set to a very large value if the velocity + ! is truncated more than MAXTRUNC times between energy saves. Set MAXTRUNC to 0 + ! to stop if there is any truncation of velocities. +ENERGYSAVEDAYS = 0.25 ! [days] default = 1.0 + ! The interval in units of TIMEUNIT between saves of the energies of the run and + ! other globally summed diagnostics. + +! === module MOM_surface_forcing === +BUOY_CONFIG = "linear" ! default = "zero" + ! The character string that indicates how buoyancy forcing is specified. Valid + ! options include (file), (zero), (linear), (USER), (BFB) and (NONE). +WIND_CONFIG = "gyres" ! default = "zero" + ! The character string that indicates how wind forcing is specified. Valid + ! options include (file), (2gyre), (1gyre), (gyres), (zero), and (USER). +TAUX_SIN_AMP = 0.1 ! [Pa] default = 0.0 + ! With the gyres wind_config, the sine amplitude in the zonal wind stress + ! profile: B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L). +TAUX_N_PIS = 1.0 ! [nondim] default = 0.0 + ! With the gyres wind_config, the number of gyres in the zonal wind stress + ! profile: n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L). +RESTOREBUOY = True ! [Boolean] default = False + ! If true, the buoyancy fluxes drive the model back toward some specified + ! surface state with a rate given by FLUXCONST. +FLUXCONST = 0.5 ! [m day-1] default = 0.0 + ! The constant that relates the restoring surface fluxes to the relative surface + ! anomalies (akin to a piston velocity). Note the non-MKS units. +SST_NORTH = 27.0 ! [deg C] default = 0.0 + ! With buoy_config linear, the sea surface temperature at the northern end of + ! the domain toward which to to restore. +SST_SOUTH = 3.0 ! [deg C] default = 0.0 + ! With buoy_config linear, the sea surface temperature at the southern end of + ! the domain toward which to to restore. +GUST_CONST = 0.02 ! [Pa] default = 0.0 + ! The background gustiness in the winds. + +! === module MOM_main (MOM_driver) === +DT_FORCING = 3600.0 ! [s] default = 900.0 + ! The time step for changing forcing, coupling with other components, or + ! potentially writing certain diagnostics. The default value is given by DT. +DAYMAX = 3.0 ! [days] + ! The final time of the whole simulation, in units of TIMEUNIT seconds. This + ! also sets the potential end time of the present run segment if the end time is + ! not set via ocean_solo_nml in input.nml. +RESTART_CONTROL = 3 ! default = 1 + ! An integer whose bits encode which restart files are written. Add 2 (bit 1) + ! for a time-stamped file, and odd (bit 0) for a non-time-stamped file. A + ! non-time-stamped restart file is saved at the end of the run segment for any + ! non-negative value. +RESTINT = 365.0 ! [days] default = 0.0 + ! The interval between saves of the restart file in units of TIMEUNIT. Use 0 + ! (the default) to not save incremental restart files at all. + +! === module MOM_write_cputime === +MAXCPU = 2.88E+04 ! [wall-clock seconds] default = -1.0 + ! The maximum amount of cpu time per processor for which MOM should run before + ! saving a restart file and quitting with a return value that indicates that a + ! further run is required to complete the simulation. If automatic restarts are + ! not desired, use a negative value for MAXCPU. MAXCPU has units of wall-clock + ! seconds, so the actual CPU time used is larger by a factor of the number of + ! processors used. + +! Debugging parameters set to non-default values +U_TRUNC_FILE = "U_velocity_truncations" ! default = "" + ! The absolute path to a file into which the accelerations leading to zonal + ! velocity truncations are written. Undefine this for efficiency if this + ! diagnostic is not needed. +V_TRUNC_FILE = "V_velocity_truncations" ! default = "" + ! The absolute path to a file into which the accelerations leading to meridional + ! velocity truncations are written. Undefine this for efficiency if this + ! diagnostic is not needed. diff --git a/.testing/p0/MOM_override b/.testing/p0/MOM_override new file mode 100644 index 0000000000..e69de29bb2 diff --git a/.testing/p0/diag_table b/.testing/p0/diag_table new file mode 100644 index 0000000000..68c71dd2c4 --- /dev/null +++ b/.testing/p0/diag_table @@ -0,0 +1,91 @@ +"MOM benchmark Experiment" +1 1 1 0 0 0 +"prog", 1,"days",1,"days","time", +#"ave_prog", 5,"days",1,"days","Time",365,"days" +#"cont", 5,"days",1,"days","Time",365,"days" + +#This is the field section of the diag_table. + +# Prognostic Ocean fields: +#========================= + +"ocean_model","u","u","prog","all",.false.,"none",2 +"ocean_model","v","v","prog","all",.false.,"none",2 +"ocean_model","h","h","prog","all",.false.,"none",1 +"ocean_model","e","e","prog","all",.false.,"none",2 +"ocean_model","temp","temp","prog","all",.false.,"none",2 +#"ocean_model","salt","salt","prog","all",.false.,"none",2 + +#"ocean_model","u","u","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","v","v","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","h","h","ave_prog_%4yr_%3dy","all",.true.,"none",1 +#"ocean_model","e","e","ave_prog_%4yr_%3dy","all",.true.,"none",2 + +# Auxilary Tracers: +#================== +#"ocean_model","vintage","vintage","prog_%4yr_%3dy","all",.false.,"none",2 +#"ocean_model","age","age","prog_%4yr_%3dy","all",.false.,"none",2 + +# Continuity Equation Terms: +#=========================== +#"ocean_model","dhdt","dhdt","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","wd","wd","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","uh","uh","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","vh","vh","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","h_rho","h_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","uh_rho","uh_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","vh_rho","vh_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","uhGM_rho","uhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","vhGM_rho","vhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 + +# +# Tracer Fluxes: +#================== +#"ocean_model","T_adx", "T_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","T_ady", "T_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","T_diffx","T_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","T_diffy","T_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_adx", "S_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_ady", "S_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_diffx","S_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_diffy","S_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 + +# testing +# ======= +#"ocean_model","Kv_u","Kv_u","prog","all",.false.,"none",2 +#"ocean_model","Kv_v","Kv_v","prog","all",.false.,"none",2 + +#============================================================================================= +# +#===- This file can be used with diag_manager/v2.0a (or higher) ==== +# +# +# FORMATS FOR FILE ENTRIES (not all input values are used) +# ------------------------ +# +#"file_name", output_freq, "output_units", format, "time_units", "time_long_name", ... +# (opt) new_file_frequecy, (opt) "new_file_freq_units", "new_file_start_date" +# +# +#output_freq: > 0 output frequency in "output_units" +# = 0 output frequency every time step +# =-1 output frequency at end of run +# +#output_units = units used for output frequency +# (years, months, days, minutes, hours, seconds) +# +#time_units = units used to label the time axis +# (days, minutes, hours, seconds) +# +# +# FORMAT FOR FIELD ENTRIES (not all input values are used) +# ------------------------ +# +#"module_name", "field_name", "output_name", "file_name" "time_sampling", time_avg, "other_opts", packing +# +#time_avg = .true. or .false. +# +#packing = 1 double precision +# = 2 float +# = 4 packed 16-bit integers +# = 8 packed 1-byte (not tested?) diff --git a/.testing/p0/input.nml b/.testing/p0/input.nml new file mode 100644 index 0000000000..41555b8822 --- /dev/null +++ b/.testing/p0/input.nml @@ -0,0 +1,22 @@ + &MOM_input_nml + output_directory = './', + input_filename = 'n' + restart_input_dir = 'INPUT/', + restart_output_dir = 'RESTART/', + parameter_filename = 'MOM_input', + 'MOM_override' / + + &diag_manager_nml + / + + &fms_nml + clock_grain='ROUTINE' + clock_flags='SYNC' + !domains_stack_size = 955296 + domains_stack_size = 14256000 + stack_size =0 / + +!&ocean_solo_nml +! hours = 1 +! !days = 1 +!/ diff --git a/.testing/tools/compare_clocks.py b/.testing/tools/compare_clocks.py new file mode 100755 index 0000000000..77198fda6a --- /dev/null +++ b/.testing/tools/compare_clocks.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python +import argparse +import json + +# Ignore timers below this threshold (in seconds) +DEFAULT_THRESHOLD = 0.05 + +# Thresholds for reporting +DT_WARN = 0.10 # Slowdown warning +DT_FAIL = 0.25 # Slowdown abort + +ANSI_RED = '\033[31m' +ANSI_GREEN = '\033[32m' +ANSI_YELLOW = '\033[33m' +ANSI_RESET = '\033[0m' + + +def main(): + desc = ( + 'Compare two FMS clock output files and report any differences within ' + 'a defined threshold.' + ) + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('expt') + parser.add_argument('ref') + parser.add_argument('--threshold') + parser.add_argument('--verbose', action='store_true') + args = parser.parse_args() + + threshold = float(args.threshold) if args.threshold else DEFAULT_THRESHOLD + verbose = args.verbose + + clock_cmp = {} + + print('{:35s}{:8s} {:8s}'.format('', 'Profile', 'Reference')) + print() + + with open(args.expt) as log_expt, open(args.ref) as log_ref: + clocks_expt = json.load(log_expt)['clocks'] + clocks_ref = json.load(log_ref)['clocks'] + + # Gather timers which appear in both clocks + clock_tags = [clk for clk in clocks_expt if clk in clocks_ref] + + for clk in clock_tags: + clock_cmp[clk] = {} + + # For now, we only comparge tavg, the rank-averaged timing + rec = 'tavg' + + t_expt = clocks_expt[clk][rec] + t_ref = clocks_ref[clk][rec] + + # Compare the relative runtimes + if all(t > threshold for t in (t_expt, t_ref)): + dclk = (t_expt - t_ref) / t_ref + else: + dclk = 0. + clock_cmp[clk][rec] = dclk + + # Skip trivially low clocks + if all(t < threshold for t in (t_expt, t_ref)) and not verbose: + continue + + # Report the time differences + ansi_color = ANSI_RESET + + if abs(t_expt - t_ref) > threshold: + if dclk > DT_FAIL: + ansi_color = ANSI_RED + elif dclk > DT_WARN: + ansi_color = ANSI_YELLOW + elif dclk < -DT_WARN: + ansi_color = ANSI_GREEN + + print('{}{}: {:7.3f}s, {:7.3f}s ({:.1f}%){}'.format( + ansi_color, + ' ' * (32 - len(clk)) + clk, + t_expt, + t_ref, + 100. * dclk, + ANSI_RESET, + )) + + +if __name__ == '__main__': + main() diff --git a/.testing/tools/compare_perf.py b/.testing/tools/compare_perf.py new file mode 100755 index 0000000000..e4a651c709 --- /dev/null +++ b/.testing/tools/compare_perf.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python +import argparse +import json + +# Ignore timers below this threshold (in seconds) +DEFAULT_THRESHOLD = 0.05 + +# Thresholds for reporting +DT_WARN = 0.10 # Slowdown warning +DT_FAIL = 0.25 # Slowdown abort + +ANSI_RED = '\033[31m' +ANSI_GREEN = '\033[32m' +ANSI_YELLOW = '\033[33m' +ANSI_RESET = '\033[0m' + + +def main(): + desc = ( + 'Compare two FMS clock output files and report any differences within ' + 'a defined threshold.' + ) + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('expt') + parser.add_argument('ref') + parser.add_argument('--threshold') + parser.add_argument('--verbose', action='store_true') + args = parser.parse_args() + + threshold = float(args.threshold) if args.threshold else DEFAULT_THRESHOLD + verbose = args.verbose + + clock_cmp = {} + + print('{:35s}{:8s} {:8s}'.format('', 'Profile', 'Reference')) + print() + + with open(args.expt) as profile_expt, open(args.ref) as profile_ref: + perf_expt = json.load(profile_expt) + perf_ref = json.load(profile_ref) + + events = [ev for ev in perf_expt if ev in perf_ref] + + for event in events: + # For now, only report the times + if event not in ('task-clock', 'cpu-clock'): + continue + + count_expt = perf_expt[event]['count'] + count_ref = perf_ref[event]['count'] + + symbols_expt = perf_expt[event]['symbol'] + symbols_ref = perf_ref[event]['symbol'] + + symbols = [ + s for s in symbols_expt + if s in symbols_ref + and not s.startswith('0x') + ] + + for symbol in symbols: + t_expt = float(symbols_expt[symbol]) / 1e9 + t_ref = float(symbols_ref[symbol]) / 1e9 + + # Compare the relative runtimes + if all(t > threshold for t in (t_expt, t_ref)): + dclk = (t_expt - t_ref) / t_ref + else: + dclk = 0. + + # Skip trivially low clocks + if all(t < threshold for t in (t_expt, t_ref)) and not verbose: + continue + + # Report the time differences + ansi_color = ANSI_RESET + + if abs(t_expt - t_ref) > threshold: + if dclk > DT_FAIL: + ansi_color = ANSI_RED + elif dclk > DT_WARN: + ansi_color = ANSI_YELLOW + elif dclk < -DT_WARN: + ansi_color = ANSI_GREEN + + # Remove module name + sname = symbol.split('_MOD_', 1)[-1] + + # Strip version from glibc calls + sname = sname.split('@')[0] + + # Remove GCC optimization renaming + sname = sname.replace('.constprop.0', '') + + if len(sname) > 32: + sname = sname[:29] + '...' + + print('{}{}: {:7.3f}s, {:7.3f}s ({:.1f}%){}'.format( + ansi_color, + ' ' * (32 - len(sname)) + sname, + t_expt, + t_ref, + 100. * dclk, + ANSI_RESET, + )) + + +if __name__ == '__main__': + main() diff --git a/.testing/tools/parse_fms_clocks.py b/.testing/tools/parse_fms_clocks.py new file mode 100755 index 0000000000..b57fc481ab --- /dev/null +++ b/.testing/tools/parse_fms_clocks.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python +import argparse +import collections +import json +import os +import sys + +import f90nml + +record_type = collections.defaultdict(lambda: float) +for rec in ('grain', 'pemin', 'pemax',): + record_type[rec] = int + + +def main(): + desc = 'Parse MOM6 model stdout and return clock data in JSON format.' + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('--format', '-f', action='store_true') + parser.add_argument('--dir', '-d') + parser.add_argument('log') + args = parser.parse_args() + + config = {} + + if args.dir: + # Gather model configuration + input_nml = os.path.join(args.dir, 'input.nml') + nml = f90nml.read(input_nml) + config['input.nml'] = nml.todict() + + parameter_filenames = [ + ('params', 'MOM_parameter_doc.all'), + ('layout', 'MOM_parameter_doc.layout'), + ('debug', 'MOM_parameter_doc.debugging'), + ] + for key, fname in parameter_filenames: + config[key] = {} + with open(os.path.join(args.dir, fname)) as param_file: + params = parse_mom6_param(param_file) + config[key].update(params) + + # Get log path + if os.path.isfile(args.log): + log_path = args.log + elif os.path.isfile(os.path.join(args.dir, args.log)): + log_path = os.path.join(args.dir, args.log) + else: + sys.exit('stdout log not found.') + + # Parse timings + with open(log_path) as log: + clocks = parse_clocks(log) + + config['clocks'] = clocks + + if args.format: + print(json.dumps(config, indent=4)) + else: + print(json.dumps(config)) + + +def parse_mom6_param(param_file): + params = {} + for line in param_file: + param_stmt = line.split('!')[0].strip() + if param_stmt: + key, val = [s.strip() for s in param_stmt.split('=')] + + # TODO: Convert to equivalent Python types + if val in ('True', 'False'): + params[key] = bool(val) + else: + params[key] = val + + return params + + +def parse_clocks(log): + clock_start_msg = 'Tabulating mpp_clock statistics across' + clock_end_msg = 'MPP_STACK high water mark=' + + fields = [] + for line in log: + if line.startswith(clock_start_msg): + npes = line.lstrip(clock_start_msg).split()[0] + + # Get records + fields = [] + line = next(log) + + # Skip blank lines + while line.isspace(): + line = next(log) + + fields = line.split() + + # Exit this loop, begin clock parsing + break + + clocks = {} + for line in log: + # Treat MPP_STACK usage as end of clock reports + if line.lstrip().startswith(clock_end_msg): + break + + record = line.split()[-len(fields):] + + clk = line.split(record[0])[0].strip() + clocks[clk] = {} + + for fld, rec in zip(fields, record): + rtype = record_type[fld] + clocks[clk][fld] = rtype(rec) + + return clocks + + +if __name__ == '__main__': + main() diff --git a/.testing/tools/parse_perf.py b/.testing/tools/parse_perf.py new file mode 100755 index 0000000000..b86b1cc106 --- /dev/null +++ b/.testing/tools/parse_perf.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python +import argparse +import collections +import json +import os +import shlex +import subprocess +import sys + + +def main(): + desc = 'Parse perf.data and return in JSON format.' + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('--format', '-f', action='store_true') + parser.add_argument('data') + args = parser.parse_args() + + profile = parse_perf_report(args.data) + + if args.format: + print(json.dumps(profile, indent=4)) + else: + print(json.dumps(profile)) + + +def parse_perf_report(perf_data_path): + profile = {} + + cmd = shlex.split( + 'perf report -s symbol,period -i {}'.format(perf_data_path) + ) + with subprocess.Popen(cmd, stdout=subprocess.PIPE, text=True) as proc: + event_name = None + for line in proc.stdout: + # Skip blank lines: + if not line or line.isspace(): + continue + + # Set the current event + if line.startswith('# Samples: '): + event_name = line.split()[-1].strip("'") + + # Remove perf modifiers for now + event_name = event_name.rsplit(':', 1)[0] + + profile[event_name] = {} + profile[event_name]['symbol'] = {} + + # Get total count + elif line.startswith('# Event count '): + event_count = int(line.split()[-1]) + profile[event_name]['count'] = event_count + + # skip all other 'comment' lines + elif line.startswith('#'): + continue + + # get per-symbol count + else: + tokens = line.split() + symbol = tokens[2] + period = int(tokens[3]) + + profile[event_name]['symbol'][symbol] = period + + return profile + + +if __name__ == '__main__': + main() From 7a650261c2765d3c60ccd931a1ed609914842079 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Sat, 24 Jul 2021 13:46:41 -0600 Subject: [PATCH 09/31] Implement visc remnant versions of PFu, CAu, u_accel_bt, diffu --- src/core/MOM_barotropic.F90 | 12 ++- src/core/MOM_dynamics_split_RK2.F90 | 85 +++++++++++++++++++ src/core/MOM_variables.F90 | 12 ++- .../lateral/MOM_hor_visc.F90 | 36 ++++++++ 4 files changed, 143 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index a8262608f8..42f91b2d8e 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2685,7 +2685,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k) enddo ; enddo ; enddo endif - + if (present(ADp) .and. (associated(ADp%visc_rem_u))) then + do k=1,nz ; do j=js,je ; do I=is-1,ie + ADp%visc_rem_u(I,j,k) = visc_rem_u(I,j,k) + enddo ; enddo ; enddo + endif + if (present(ADp) .and. (associated(ADp%visc_rem_u))) then + do k=1,nz ; do J=js-1,je ; do i=is,ie + ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k) + enddo ; enddo ; enddo + endif + if (G%nonblocking_updates) then if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain) call complete_group_pass(CS%pass_ubta_uhbta, G%Domain) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 07a302b2b0..c55dbce684 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -163,10 +163,12 @@ module MOM_dynamics_split_RK2 integer :: id_h_PFu = -1, id_h_PFv = -1 integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1 integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1 + integer :: id_PFu_visc_rem = -1, id_PFv_visc_rem = -1 ! integer :: id_hf_CAu = -1, id_hf_CAv = -1 integer :: id_h_CAu = -1, id_h_CAv = -1 integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1 integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1 + integer :: id_CAu_visc_rem = -1, id_CAv_visc_rem = -1 ! Split scheme only. integer :: id_uav = -1, id_vav = -1 @@ -175,6 +177,7 @@ module MOM_dynamics_split_RK2 integer :: id_h_u_BT_accel = -1, id_h_v_BT_accel = -1 integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1 integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1 + integer :: id_u_BT_accel_visc_rem = -1, id_v_BT_accel_visc_rem = -1 !>@} type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the @@ -353,6 +356,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [H L T-2 ~> m2 s-2]. real, dimension(SZI_(G),SZJB_(G)) :: & intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2]. + + ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv], + real, allocatable, dimension(:,:,:) :: & + PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2]. + CAu_visc_rem, CAv_visc_rem, & ! Coriolis force accel. x visc_rem_[uv] [L T-2 ~> m s-2]. + u_BT_accel_visc_rem, v_BT_accel_visc_rem ! barotropic correction accel. x visc_rem_[uv] [L T-2 ~> m s-2]. real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. @@ -1102,6 +1111,61 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s deallocate(h_v_BT_accel) endif + if (CS%id_PFu_visc_rem > 0) then + allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + PFu_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag) + deallocate(PFu_visc_rem) + endif + if (CS%id_PFv_visc_rem > 0) then + allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + PFv_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag) + deallocate(PFv_visc_rem) + endif + if (CS%id_CAu_visc_rem > 0) then + allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + CAu_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag) + deallocate(CAu_visc_rem) + endif + if (CS%id_CAv_visc_rem > 0) then + allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + CAv_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag) + deallocate(CAv_visc_rem) + endif + if (CS%id_u_BT_accel_visc_rem > 0) then + allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + u_BT_accel_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag) + deallocate(u_BT_accel_visc_rem) + endif + if (CS%id_v_BT_accel_visc_rem > 0) then + allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + v_BT_accel_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag) + deallocate(v_BT_accel_visc_rem) + endif + if (CS%debug) then call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) @@ -1605,6 +1669,27 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 'Depth-integral of Barotropic Anomaly Meridional Acceleration', & 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) + + CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, & + 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, & + 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + + CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, & + 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, & + 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + + CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, & + 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, & + 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f966ab2ad2..c81d1f6f8d 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -163,16 +163,24 @@ module MOM_variables real, pointer, dimension(:,:,:) :: & diffu => NULL(), & !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] diffv => NULL(), & !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] + diffu_visc_rem => NULL(), & !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] + diffv_visc_rem => NULL(), & !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] CAu => NULL(), & !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2] CAv => NULL(), & !< Meridional Coriolis and momentum advection accelerations [L T-2 ~> m s-2] + CAu_visc_rem => NULL(), & !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] + CAv_visc_rem => NULL(), & !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] PFu => NULL(), & !< Zonal acceleration due to pressure forces [L T-2 ~> m s-2] PFv => NULL(), & !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2] + PFu_visc_rem => NULL(), & !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] + PFv_visc_rem => NULL(), & !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2] dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2] du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal mixing [L T-2 ~> m s-2] dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal mixing [L T-2 ~> m s-2] u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2] - v_accel_bt => NULL() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2] + v_accel_bt => NULL(), & !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2] + u_accel_bt_visc_rem => NULL(), &!< Pointer to the zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] + v_accel_bt_visc_rem => NULL() !< Pointer to the meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] real, pointer, dimension(:,:,:) :: du_other => NULL() !< Zonal velocity changes due to any other processes that are !! not due to any explicit accelerations [L T-1 ~> m s-1]. @@ -191,6 +199,8 @@ module MOM_variables real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points + real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points + real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points end type accel_diag_ptrs !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances. diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 2324970254..bafce4b209 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -190,6 +190,7 @@ module MOM_hor_visc integer :: id_h_diffu = -1, id_h_diffv = -1 integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1 integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1 + integer :: id_diffu_visc_rem = -1, id_diffv_visc_rem = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 @@ -280,6 +281,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: h_diffv ! h x diffv [H L T-2 ~> m2 s-2] + real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m2 s-2] + real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m2 s-2] real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] @@ -1703,6 +1706,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call post_data(CS%id_h_diffv, h_diffv, CS%diag) deallocate(h_diffv) endif + + if (present(ADp) .and. (CS%id_diffu_visc_rem > 0)) then + allocate(diffu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + diffu_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + diffu_visc_rem(I,j,k) = diffu(I,j,k) * ADp%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_diffu_visc_rem, diffu_visc_rem, CS%diag) + deallocate(diffu_visc_rem) + endif + if (present(ADp) .and. (CS%id_diffv_visc_rem > 0)) then + allocate(diffv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + diffv_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + diffv_visc_rem(i,J,k) = diffv(i,J,k) * ADp%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_diffv_visc_rem, diffv_visc_rem, CS%diag) + deallocate(diffv_visc_rem) + endif end subroutine horizontal_viscosity @@ -2447,6 +2469,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) endif + + CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, & + 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke) + endif + + CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, & + 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) + endif if (CS%biharmonic) then CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, & From 2601b47285ebcc9ed6fc3742a960f43c92d38809 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Sat, 24 Jul 2021 16:03:01 -0600 Subject: [PATCH 10/31] Add d[uv]_dt_visc_rem diagnostics --- src/core/MOM_variables.F90 | 4 ++ src/diagnostics/MOM_diagnostics.F90 | 7 +++ .../vertical/MOM_vert_friction.F90 | 53 +++++++++++++++++++ 3 files changed, 64 insertions(+) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index c81d1f6f8d..7272a5994f 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -161,6 +161,8 @@ module MOM_variables ! Each of the following fields has nz layers. real, pointer, dimension(:,:,:) :: & + du_dt => NULL(), & !< Zonal acceleration [L T-2 ~> m s-2] + dv_dt => NULL(), & !< Meridional acceleration [L T-2 ~> m s-2] diffu => NULL(), & !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] diffv => NULL(), & !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] diffu_visc_rem => NULL(), & !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] @@ -175,6 +177,8 @@ module MOM_variables PFv_visc_rem => NULL(), & !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2] dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2] + du_dt_visc_rem => NULL(), &!< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] + dv_dt_visc_rem => NULL(), &!< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal mixing [L T-2 ~> m s-2] dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal mixing [L T-2 ~> m s-2] u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2] diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 1d59969011..00ed7682d7 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -290,6 +290,13 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state) if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state) + + do k=1,nz ; do j=js,je ; do I=is-1,ie + ADp%du_dt(I,j,k) = CS%du_dt(I,j,k) + enddo ; enddo ; enddo + do k=1,nz ; do J=js-1,je ; do i=is,ie + ADp%dv_dt(i,J,k) = CS%dv_dt(i,J,k) + enddo ; enddo ; enddo !! Diagnostics for terms multiplied by fractional thicknesses diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 081179eb41..d958188a00 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -127,6 +127,7 @@ module MOM_vert_friction ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1 integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1 integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1 + integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1 !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure @@ -216,6 +217,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2] + + real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u) [L T-2 ~> m2 s-2] + real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v) [L T-2 ~> m2 s-2] logical :: do_i(SZIB_(G)) logical :: DoStokesMixing @@ -335,6 +339,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt enddo ; enddo ; endif + + !if (associated(ADp%du_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq + ! ADp%du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + & + ! (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k) + !enddo ; enddo ; endif if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq visc%taux_shelf(I,j) = -GV%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf? @@ -416,6 +425,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt enddo ; enddo ; endif + + !if (associated(ADp%dv_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq + ! ADp%dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + & + ! (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k) + !enddo ; enddo ; endif if (associated(visc%tauy_shelf)) then ; do i=is,ie visc%tauy_shelf(i,J) = -GV%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf? @@ -521,6 +535,27 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag) deallocate(h_dv_dt_visc) endif + + if (CS%id_du_dt_visc_rem > 0) then + allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + du_dt_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + & + (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag) + deallocate(du_dt_visc_rem) + endif + if (CS%id_dv_dt_visc_rem > 0) then + allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + dv_dt_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + & + (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag) + deallocate(dv_dt_visc_rem) + endif end subroutine vertvisc @@ -1876,6 +1911,24 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz) endif + + CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, & + 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_du_dt_visc_rem > 0) then + call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, & + 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_dv_dt_visc_rem > 0) then + call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) + endif if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) & call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp) From 20b8bfc9e6cad5fd773e0d1bffe64a6d13e7bf84 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Sun, 25 Jul 2021 16:14:11 -0600 Subject: [PATCH 11/31] Attempt to implement KE budget modified by visc rem --- src/core/MOM_variables.F90 | 24 ++-- src/diagnostics/MOM_diagnostics.F90 | 178 +++++++++++++++++++++++++++- 2 files changed, 185 insertions(+), 17 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 7272a5994f..2f9c3cf643 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -165,26 +165,16 @@ module MOM_variables dv_dt => NULL(), & !< Meridional acceleration [L T-2 ~> m s-2] diffu => NULL(), & !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] diffv => NULL(), & !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] - diffu_visc_rem => NULL(), & !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] - diffv_visc_rem => NULL(), & !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] CAu => NULL(), & !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2] CAv => NULL(), & !< Meridional Coriolis and momentum advection accelerations [L T-2 ~> m s-2] - CAu_visc_rem => NULL(), & !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] - CAv_visc_rem => NULL(), & !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] PFu => NULL(), & !< Zonal acceleration due to pressure forces [L T-2 ~> m s-2] PFv => NULL(), & !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2] - PFu_visc_rem => NULL(), & !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] - PFv_visc_rem => NULL(), & !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2] dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2] - du_dt_visc_rem => NULL(), &!< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] - dv_dt_visc_rem => NULL(), &!< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal mixing [L T-2 ~> m s-2] dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal mixing [L T-2 ~> m s-2] u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2] - v_accel_bt => NULL(), & !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2] - u_accel_bt_visc_rem => NULL(), &!< Pointer to the zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] - v_accel_bt_visc_rem => NULL() !< Pointer to the meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] + v_accel_bt => NULL() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2] real, pointer, dimension(:,:,:) :: du_other => NULL() !< Zonal velocity changes due to any other processes that are !! not due to any explicit accelerations [L T-1 ~> m s-1]. @@ -205,6 +195,18 @@ module MOM_variables real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points + + real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] + end type accel_diag_ptrs !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances. diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 00ed7682d7..28b5764079 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -108,7 +108,16 @@ module MOM_diagnostics KE_adv => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3] KE_visc => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] - KE_dia => NULL() !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] + KE_dia => NULL(), & !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] + + ! The following arrays hold diagnostics in the modified layer-integrated energy budget. + ! Modification is through using the visc_rem_[uv]-filtered momentum equation + PE_to_KE_visc_rem => NULL(), & !< potential energy to KE term [m3 s-3] + KE_BT_visc_rem => NULL(), & !< barotropic contribution to KE term [m3 s-3] + KE_CorAdv_visc_rem => NULL(), & !< KE source from the combined Coriolis and + !! advection terms [H L2 T-3 ~> m3 s-3]. + KE_visc_rem => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] + KE_horvisc_rem => NULL() !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] !>@{ Diagnostic IDs integer :: id_u = -1, id_v = -1, id_h = -1 @@ -124,6 +133,10 @@ module MOM_diagnostics integer :: id_KE_Coradv = -1 integer :: id_KE_adv = -1, id_KE_visc = -1 integer :: id_KE_horvisc = -1, id_KE_dia = -1 + integer :: id_PE_to_KE_visc_rem = -1, id_KE_BT_visc_rem = -1 + integer :: id_KE_Coradv_visc_rem = -1 + integer :: id_KE_visc_rem = -1 + integer :: id_KE_horvisc_rem = -1 integer :: id_uh_Rlay = -1, id_vh_Rlay = -1 integer :: id_uhGM_Rlay = -1, id_vhGM_Rlay = -1 integer :: id_h_Rlay = -1, id_Rd1 = -1 @@ -1055,7 +1068,10 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (.not.G%symmetric) then if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. & associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. & - associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) ) then + associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) .OR. & + associated(CS%PE_to_KE_visc_rem) .OR. associated(CS%KE_BT_visc_rem) .OR. & + associated(CS%KE_CorAdv_visc_rem) .OR. associated(CS%KE_visc_rem) .OR. & + associated(CS%KE_horvisc_rem)) then call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East) endif endif @@ -1223,6 +1239,100 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS enddo if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag) endif + + if (associated(CS%PE_to_KE_visc_rem)) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu_visc_rem(I,j,k) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv_visc_rem(i,J,k) + enddo ; enddo + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do j=js,je ; do i=is,ie + CS%PE_to_KE_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + enddo ; enddo + enddo + if (CS%id_PE_to_KE_visc_rem > 0) call post_data(CS%id_PE_to_KE_visc_rem, CS%PE_to_KE_visc_rem, CS%diag) + endif + + if (associated(CS%KE_BT_visc_rem)) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt_visc_rem(I,j,k) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt_visc_rem(i,J,k) + enddo ; enddo + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do j=js,je ; do i=is,ie + CS%KE_BT_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + enddo ; enddo + enddo + if (CS%id_KE_BT_visc_rem > 0) call post_data(CS%id_KE_BT_visc_rem, CS%KE_BT_visc_rem, CS%diag) + endif + + if (associated(CS%KE_CorAdv_visc_rem)) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%CAu_visc_rem(I,j,k) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%CAv_visc_rem(i,J,k) + enddo ; enddo + do j=js,je ; do i=is,ie + KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) & + * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) + enddo ; enddo + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do j=js,je ; do i=is,ie + CS%KE_CorAdv_visc_rem(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + enddo ; enddo + enddo + if (CS%id_KE_CorAdv_visc_rem > 0) call post_data(CS%id_KE_Coradv_visc_rem, CS%KE_Coradv_visc_rem, CS%diag) + endif + + if (associated(CS%KE_visc_rem)) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_rem(I,j,k) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc_rem(i,J,k) + enddo ; enddo + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do j=js,je ; do i=is,ie + CS%KE_visc(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + enddo ; enddo + enddo + if (CS%id_KE_visc_rem > 0) call post_data(CS%id_KE_visc_rem, CS%KE_visc_rem, CS%diag) + endif + + if (associated(CS%KE_horvisc_rem)) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu_visc_rem(I,j,k) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv_visc_rem(i,J,k) + enddo ; enddo + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do j=js,je ; do i=is,ie + CS%KE_horvisc_rem(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + enddo ; enddo + enddo + if (CS%id_KE_horvisc_rem > 0) call post_data(CS%id_KE_horvisc_rem, CS%KE_horvisc_rem, CS%diag) + endif end subroutine calculate_energy_diagnostics @@ -1883,18 +1993,31 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'Potential to Kinetic Energy Conversion of Layer', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz) + + CS%id_PE_to_KE_visc_rem = register_diag_field('ocean_model', 'PE_to_KE_visc_rem', diag%axesTL, Time, & + 'Potential to Kinetic Energy Conversion multiplied by viscous remnant fraction', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_PE_to_KE_visc_rem>0) call safe_alloc_ptr(CS%PE_to_KE_visc_rem,isd,ied,jsd,jed,nz) if (split) then CS%id_KE_BT = register_diag_field('ocean_model', 'KE_BT', diag%axesTL, Time, & 'Barotropic contribution to Kinetic Energy', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_KE_BT>0) call safe_alloc_ptr(CS%KE_BT,isd,ied,jsd,jed,nz) + if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz) + CS%id_KE_BT_visc_rem = register_diag_field('ocean_model', 'KE_BT_visc_rem', diag%axesTL, Time, & + 'Barotropic contribution to Kinetic Energy multiplied by viscous remnant fraction', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz) endif CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, & 'Kinetic Energy Source from Coriolis and Advection', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_Coradv>0) call safe_alloc_ptr(CS%KE_Coradv,isd,ied,jsd,jed,nz) + CS%id_KE_Coradv_visc_rem = register_diag_field('ocean_model', 'KE_Coradv_visc_rem', diag%axesTL, Time, & + 'Kinetic Energy Source from Coriolis and Advection multiplied by viscous remnant fraction', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_KE_Coradv_visc_rem>0) call safe_alloc_ptr(CS%KE_Coradv_visc_rem,isd,ied,jsd,jed,nz) CS%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, Time, & 'Kinetic Energy Source from Advection', & @@ -1905,11 +2028,19 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'Kinetic Energy Source from Vertical Viscosity and Stresses', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz) + CS%id_KE_visc_rem = register_diag_field('ocean_model', 'KE_visc_rem', diag%axesTL, Time, & + 'Kinetic Energy Source from Vertical Viscosity and Stresses multiplied by viscous remnant fraction', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_KE_visc_rem>0) call safe_alloc_ptr(CS%KE_visc_rem,isd,ied,jsd,jed,nz) CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, & 'Kinetic Energy Source from Horizontal Viscosity', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_horvisc>0) call safe_alloc_ptr(CS%KE_horvisc,isd,ied,jsd,jed,nz) + CS%id_KE_horvisc_rem = register_diag_field('ocean_model', 'KE_horvisc_rem', diag%axesTL, Time, & + 'Kinetic Energy Source from Horizontal Viscosity multiplied by viscous remnant fraction', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_KE_horvisc_rem>0) call safe_alloc_ptr(CS%KE_horvisc_rem,isd,ied,jsd,jed,nz) if (.not. adiabatic) then CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, & @@ -2307,8 +2438,11 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. & associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. & associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. & - associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) & + associated(CS%KE_horvisc) .or. associated(CS%KE_dia) .or. & + associated(CS%PE_to_KE_visc_rem) .or. & associated(CS%KE_BT_visc_rem) .or. & + associated(CS%KE_CorAdv_visc_rem) .or. associated(CS%KE_horvisc_rem)) then call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) + endif if (associated(CS%dKE_dt)) then if (.not.associated(CS%du_dt)) then @@ -2329,11 +2463,30 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) call safe_alloc_ptr(ADp%gradKEu,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%gradKEv,isd,ied,JsdB,JedB,nz) endif - if (associated(CS%KE_visc)) then call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) endif + if (associated(CS%KE_CorAdv_visc_rem)) then + call safe_alloc_ptr(ADp%CAu_visc_rem,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%CAv_visc_rem,isd,ied,JsdB,JedB,nz) + endif + if (associated(CS%PE_to_KE_visc_rem)) then + call safe_alloc_ptr(ADp%PFu_visc_rem,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%PFv_visc_rem,isd,ied,JsdB,JedB,nz) + endif + if (associated(CS%KE_BT_visc_rem)) then + call safe_alloc_ptr(ADp%u_accel_BT_visc_rem,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%v_accel_BT_visc_rem,isd,ied,JsdB,JedB,nz) + endif + if (associated(CS%KE_visc_rem)) then + call safe_alloc_ptr(ADp%du_dt_visc_rem,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%dv_dt_visc_rem,isd,ied,JsdB,JedB,nz) + endif + if (associated(CS%KE_horvisc_rem)) then + call safe_alloc_ptr(ADp%diffu_visc_rem,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%diffv_visc_rem,isd,ied,JsdB,JedB,nz) + endif if (associated(CS%KE_dia)) then call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) @@ -2359,11 +2512,16 @@ subroutine MOM_diagnostics_end(CS, ADp) if (associated(CS%KE)) deallocate(CS%KE) if (associated(CS%dKE_dt)) deallocate(CS%dKE_dt) if (associated(CS%PE_to_KE)) deallocate(CS%PE_to_KE) + if (associated(CS%PE_to_KE_visc_rem)) deallocate(CS%PE_to_KE_visc_rem) if (associated(CS%KE_BT)) deallocate(CS%KE_BT) + if (associated(CS%KE_BT_visc_rem)) deallocate(CS%KE_BT_visc_rem) if (associated(CS%KE_Coradv)) deallocate(CS%KE_Coradv) + if (associated(CS%KE_Coradv_visc_rem)) deallocate(CS%KE_Coradv_visc_rem) if (associated(CS%KE_adv)) deallocate(CS%KE_adv) if (associated(CS%KE_visc)) deallocate(CS%KE_visc) + if (associated(CS%KE_visc_rem)) deallocate(CS%KE_visc_rem) if (associated(CS%KE_horvisc)) deallocate(CS%KE_horvisc) + if (associated(CS%KE_horvisc_rem)) deallocate(CS%KE_horvisc_rem) if (associated(CS%KE_dia)) deallocate(CS%KE_dia) if (associated(CS%dv_dt)) deallocate(CS%dv_dt) if (associated(CS%dh_dt)) deallocate(CS%dh_dt) @@ -2375,13 +2533,21 @@ subroutine MOM_diagnostics_end(CS, ADp) if (associated(CS%vhGM_Rlay)) deallocate(CS%vhGM_Rlay) if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu) - if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu) + if (associated(ADp%gradKEv)) deallocate(ADp%gradKEv) if (associated(ADp%du_dt_visc)) deallocate(ADp%du_dt_visc) if (associated(ADp%dv_dt_visc)) deallocate(ADp%dv_dt_visc) if (associated(ADp%du_dt_dia)) deallocate(ADp%du_dt_dia) if (associated(ADp%dv_dt_dia)) deallocate(ADp%dv_dt_dia) if (associated(ADp%du_other)) deallocate(ADp%du_other) if (associated(ADp%dv_other)) deallocate(ADp%dv_other) + if (associated(ADp%CAu_visc_rem)) deallocate(ADp%CAu_visc_rem) + if (associated(ADp%CAv_visc_rem)) deallocate(ADp%CAv_visc_rem) + if (associated(ADp%PFu_visc_rem)) deallocate(ADp%PFu_visc_rem) + if (associated(ADp%PFv_visc_rem)) deallocate(ADp%PFv_visc_rem) + if (associated(ADp%u_accel_bt_visc_rem)) deallocate(ADp%u_accel_bt_visc_rem) + if (associated(ADp%v_accel_bt_visc_rem)) deallocate(ADp%v_accel_bt_visc_rem) + if (associated(ADp%du_dt_visc_rem)) deallocate(ADp%du_dt_visc_rem) + if (associated(ADp%dv_dt_visc_rem)) deallocate(ADp%dv_dt_visc_rem) if (associated(ADp%diag_hfrac_u)) deallocate(ADp%diag_hfrac_u) if (associated(ADp%diag_hfrac_v)) deallocate(ADp%diag_hfrac_v) From 4208525a1400641f941ac98953520cf8a4875372 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Sat, 31 Jul 2021 12:56:31 -0600 Subject: [PATCH 12/31] Fix syntax error --- src/diagnostics/MOM_diagnostics.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 28b5764079..549cf5d0e2 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -2438,9 +2438,12 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. & associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. & associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. & - associated(CS%KE_horvisc) .or. associated(CS%KE_dia) .or. & - associated(CS%PE_to_KE_visc_rem) .or. & associated(CS%KE_BT_visc_rem) .or. & - associated(CS%KE_CorAdv_visc_rem) .or. associated(CS%KE_horvisc_rem)) then + associated(CS%KE_horvisc) .or. & + associated(CS%KE_dia) .or. & + associated(CS%PE_to_KE_visc_rem) .or. & + associated(CS%KE_BT_visc_rem) .or. & + associated(CS%KE_CorAdv_visc_rem) .or. & + associated(CS%KE_horvisc_rem)) then call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) endif From f4d48393a7cd160bbf086e30f20e900af6205073 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Sat, 31 Jul 2021 13:46:28 -0600 Subject: [PATCH 13/31] Fix some style errors --- src/core/MOM_barotropic.F90 | 2 +- src/core/MOM_dynamics_split_RK2.F90 | 8 ++--- src/core/MOM_variables.F90 | 31 ++++++++++++------- src/diagnostics/MOM_diagnostics.F90 | 8 ++--- .../lateral/MOM_hor_visc.F90 | 4 +-- .../vertical/MOM_vert_friction.F90 | 19 +++--------- 6 files changed, 35 insertions(+), 37 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 42f91b2d8e..5b6fa03bb8 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2695,7 +2695,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k) enddo ; enddo ; enddo endif - + if (G%nonblocking_updates) then if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain) call complete_group_pass(CS%pass_ubta_uhbta, G%Domain) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index c55dbce684..52cd711ed9 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -356,7 +356,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [H L T-2 ~> m2 s-2]. real, dimension(SZI_(G),SZJB_(G)) :: & intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2]. - + ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv], real, allocatable, dimension(:,:,:) :: & PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2]. @@ -1669,21 +1669,21 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 'Depth-integral of Barotropic Anomaly Meridional Acceleration', & 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) - + CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, & 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) - + CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, & 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, & 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) - + CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, & 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 2f9c3cf643..b33eb6f92f 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -195,17 +195,26 @@ module MOM_variables real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points - - real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity + ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity + ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations + ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations + ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied + ! by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces + ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied + ! by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied + ! by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied + ! by viscous remnant fraction [L T-2 ~> m s-2] + real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied + ! by viscous remnant fraction [L T-2 ~> m s-2] end type accel_diag_ptrs diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 549cf5d0e2..0c5412bd72 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -109,7 +109,6 @@ module MOM_diagnostics KE_visc => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] KE_dia => NULL(), & !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] - ! The following arrays hold diagnostics in the modified layer-integrated energy budget. ! Modification is through using the visc_rem_[uv]-filtered momentum equation PE_to_KE_visc_rem => NULL(), & !< potential energy to KE term [m3 s-3] @@ -303,7 +302,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state) if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state) - do k=1,nz ; do j=js,je ; do I=is-1,ie ADp%du_dt(I,j,k) = CS%du_dt(I,j,k) enddo ; enddo ; enddo @@ -1239,7 +1237,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS enddo if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag) endif - + if (associated(CS%PE_to_KE_visc_rem)) then do k=1,nz do j=js,je ; do I=Isq,Ieq @@ -2439,8 +2437,8 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. & associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. & associated(CS%KE_horvisc) .or. & - associated(CS%KE_dia) .or. & - associated(CS%PE_to_KE_visc_rem) .or. & + associated(CS%KE_dia) .or. & + associated(CS%PE_to_KE_visc_rem) .or. & associated(CS%KE_BT_visc_rem) .or. & associated(CS%KE_CorAdv_visc_rem) .or. & associated(CS%KE_horvisc_rem)) then diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index bafce4b209..ad32f81f71 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1706,7 +1706,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call post_data(CS%id_h_diffv, h_diffv, CS%diag) deallocate(h_diffv) endif - + if (present(ADp) .and. (CS%id_diffu_visc_rem > 0)) then allocate(diffu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) diffu_visc_rem(:,:,:) = 0.0 @@ -2469,7 +2469,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) endif - + CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, & 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index d958188a00..4549a5f190 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -217,9 +217,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2] - - real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u) [L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v) [L T-2 ~> m2 s-2] + + real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u) + ! [L T-2 ~> m2 s-2] + real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v) + ! [L T-2 ~> m2 s-2] logical :: do_i(SZIB_(G)) logical :: DoStokesMixing @@ -339,11 +341,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt enddo ; enddo ; endif - - !if (associated(ADp%du_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq - ! ADp%du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + & - ! (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k) - !enddo ; enddo ; endif if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq visc%taux_shelf(I,j) = -GV%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf? @@ -425,11 +422,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt enddo ; enddo ; endif - - !if (associated(ADp%dv_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq - ! ADp%dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + & - ! (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k) - !enddo ; enddo ; endif if (associated(visc%tauy_shelf)) then ; do i=is,ie visc%tauy_shelf(i,J) = -GV%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf? @@ -1920,7 +1912,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) endif - CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, & 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) From b1357bc2933e42fee71f64d5224fe6fb3a29d098 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Sat, 31 Jul 2021 14:00:47 -0600 Subject: [PATCH 14/31] Fix more style errors --- src/core/MOM_variables.F90 | 18 +++++++++--------- src/diagnostics/MOM_diagnostics.F90 | 2 +- .../vertical/MOM_vert_friction.F90 | 4 ++-- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index b33eb6f92f..8a1d0c8c60 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -195,25 +195,25 @@ module MOM_variables real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points - real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity + real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity + real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations + real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations + real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied + real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces + real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied + real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied + real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied + real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied ! by viscous remnant fraction [L T-2 ~> m s-2] end type accel_diag_ptrs diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 0c5412bd72..a81fbca4da 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1991,7 +1991,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'Potential to Kinetic Energy Conversion of Layer', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz) - + CS%id_PE_to_KE_visc_rem = register_diag_field('ocean_model', 'PE_to_KE_visc_rem', diag%axesTL, Time, & 'Potential to Kinetic Energy Conversion multiplied by viscous remnant fraction', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 4549a5f190..83e5a391fa 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -527,7 +527,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag) deallocate(h_dv_dt_visc) endif - + if (CS%id_du_dt_visc_rem > 0) then allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) du_dt_visc_rem(:,:,:) = 0.0 @@ -1903,7 +1903,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz) endif - + CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, & 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) From be09de03ebd114b43bfca5f2c1ae2687f849366c Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Tue, 3 Aug 2021 18:58:58 -0600 Subject: [PATCH 15/31] Remove added KE budget diagnostics --- src/core/MOM_variables.F90 | 20 ---- src/diagnostics/MOM_diagnostics.F90 | 175 +--------------------------- 2 files changed, 4 insertions(+), 191 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 8a1d0c8c60..6f6cf12214 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -195,26 +195,6 @@ module MOM_variables real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points - real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity - ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity - ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations - ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations - ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied - ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces - ! multiplied by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied - ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied - ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied - ! by viscous remnant fraction [L T-2 ~> m s-2] - real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied - ! by viscous remnant fraction [L T-2 ~> m s-2] end type accel_diag_ptrs diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index a81fbca4da..88a805fd04 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -108,15 +108,7 @@ module MOM_diagnostics KE_adv => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3] KE_visc => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] - KE_dia => NULL(), & !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] - ! The following arrays hold diagnostics in the modified layer-integrated energy budget. - ! Modification is through using the visc_rem_[uv]-filtered momentum equation - PE_to_KE_visc_rem => NULL(), & !< potential energy to KE term [m3 s-3] - KE_BT_visc_rem => NULL(), & !< barotropic contribution to KE term [m3 s-3] - KE_CorAdv_visc_rem => NULL(), & !< KE source from the combined Coriolis and - !! advection terms [H L2 T-3 ~> m3 s-3]. - KE_visc_rem => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] - KE_horvisc_rem => NULL() !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] + KE_dia => NULL(), !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] !>@{ Diagnostic IDs integer :: id_u = -1, id_v = -1, id_h = -1 @@ -132,10 +124,6 @@ module MOM_diagnostics integer :: id_KE_Coradv = -1 integer :: id_KE_adv = -1, id_KE_visc = -1 integer :: id_KE_horvisc = -1, id_KE_dia = -1 - integer :: id_PE_to_KE_visc_rem = -1, id_KE_BT_visc_rem = -1 - integer :: id_KE_Coradv_visc_rem = -1 - integer :: id_KE_visc_rem = -1 - integer :: id_KE_horvisc_rem = -1 integer :: id_uh_Rlay = -1, id_vh_Rlay = -1 integer :: id_uhGM_Rlay = -1, id_vhGM_Rlay = -1 integer :: id_h_Rlay = -1, id_Rd1 = -1 @@ -1066,10 +1054,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (.not.G%symmetric) then if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. & associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. & - associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) .OR. & - associated(CS%PE_to_KE_visc_rem) .OR. associated(CS%KE_BT_visc_rem) .OR. & - associated(CS%KE_CorAdv_visc_rem) .OR. associated(CS%KE_visc_rem) .OR. & - associated(CS%KE_horvisc_rem)) then + associated(CS%KE_horvisc) .OR. associated(CS%KE_dia)) then call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East) endif endif @@ -1238,100 +1223,6 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag) endif - if (associated(CS%PE_to_KE_visc_rem)) then - do k=1,nz - do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu_visc_rem(I,j,k) - enddo ; enddo - do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv_visc_rem(i,J,k) - enddo ; enddo - if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) - do j=js,je ; do i=is,ie - CS%PE_to_KE_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) - enddo ; enddo - enddo - if (CS%id_PE_to_KE_visc_rem > 0) call post_data(CS%id_PE_to_KE_visc_rem, CS%PE_to_KE_visc_rem, CS%diag) - endif - - if (associated(CS%KE_BT_visc_rem)) then - do k=1,nz - do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt_visc_rem(I,j,k) - enddo ; enddo - do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt_visc_rem(i,J,k) - enddo ; enddo - if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) - do j=js,je ; do i=is,ie - CS%KE_BT_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) - enddo ; enddo - enddo - if (CS%id_KE_BT_visc_rem > 0) call post_data(CS%id_KE_BT_visc_rem, CS%KE_BT_visc_rem, CS%diag) - endif - - if (associated(CS%KE_CorAdv_visc_rem)) then - do k=1,nz - do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%CAu_visc_rem(I,j,k) - enddo ; enddo - do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%CAv_visc_rem(i,J,k) - enddo ; enddo - do j=js,je ; do i=is,ie - KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) & - * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) - enddo ; enddo - if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) - do j=js,je ; do i=is,ie - CS%KE_CorAdv_visc_rem(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) - enddo ; enddo - enddo - if (CS%id_KE_CorAdv_visc_rem > 0) call post_data(CS%id_KE_Coradv_visc_rem, CS%KE_Coradv_visc_rem, CS%diag) - endif - - if (associated(CS%KE_visc_rem)) then - do k=1,nz - do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_rem(I,j,k) - enddo ; enddo - do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc_rem(i,J,k) - enddo ; enddo - if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) - do j=js,je ; do i=is,ie - CS%KE_visc(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) - enddo ; enddo - enddo - if (CS%id_KE_visc_rem > 0) call post_data(CS%id_KE_visc_rem, CS%KE_visc_rem, CS%diag) - endif - - if (associated(CS%KE_horvisc_rem)) then - do k=1,nz - do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu_visc_rem(I,j,k) - enddo ; enddo - do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv_visc_rem(i,J,k) - enddo ; enddo - if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) - do j=js,je ; do i=is,ie - CS%KE_horvisc_rem(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) - enddo ; enddo - enddo - if (CS%id_KE_horvisc_rem > 0) call post_data(CS%id_KE_horvisc_rem, CS%KE_horvisc_rem, CS%diag) - endif - end subroutine calculate_energy_diagnostics !> This subroutine registers fields to calculate a diagnostic time derivative. @@ -1992,30 +1883,17 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz) - CS%id_PE_to_KE_visc_rem = register_diag_field('ocean_model', 'PE_to_KE_visc_rem', diag%axesTL, Time, & - 'Potential to Kinetic Energy Conversion multiplied by viscous remnant fraction', & - 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_PE_to_KE_visc_rem>0) call safe_alloc_ptr(CS%PE_to_KE_visc_rem,isd,ied,jsd,jed,nz) - if (split) then CS%id_KE_BT = register_diag_field('ocean_model', 'KE_BT', diag%axesTL, Time, & 'Barotropic contribution to Kinetic Energy', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz) - CS%id_KE_BT_visc_rem = register_diag_field('ocean_model', 'KE_BT_visc_rem', diag%axesTL, Time, & - 'Barotropic contribution to Kinetic Energy multiplied by viscous remnant fraction', & - 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz) + if (CS%id_KE_BT>0) call safe_alloc_ptr(CS%KE_BT,isd,ied,jsd,jed,nz) endif CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, & 'Kinetic Energy Source from Coriolis and Advection', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_Coradv>0) call safe_alloc_ptr(CS%KE_Coradv,isd,ied,jsd,jed,nz) - CS%id_KE_Coradv_visc_rem = register_diag_field('ocean_model', 'KE_Coradv_visc_rem', diag%axesTL, Time, & - 'Kinetic Energy Source from Coriolis and Advection multiplied by viscous remnant fraction', & - 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_KE_Coradv_visc_rem>0) call safe_alloc_ptr(CS%KE_Coradv_visc_rem,isd,ied,jsd,jed,nz) CS%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, Time, & 'Kinetic Energy Source from Advection', & @@ -2026,19 +1904,11 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'Kinetic Energy Source from Vertical Viscosity and Stresses', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz) - CS%id_KE_visc_rem = register_diag_field('ocean_model', 'KE_visc_rem', diag%axesTL, Time, & - 'Kinetic Energy Source from Vertical Viscosity and Stresses multiplied by viscous remnant fraction', & - 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_KE_visc_rem>0) call safe_alloc_ptr(CS%KE_visc_rem,isd,ied,jsd,jed,nz) CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, & 'Kinetic Energy Source from Horizontal Viscosity', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_horvisc>0) call safe_alloc_ptr(CS%KE_horvisc,isd,ied,jsd,jed,nz) - CS%id_KE_horvisc_rem = register_diag_field('ocean_model', 'KE_horvisc_rem', diag%axesTL, Time, & - 'Kinetic Energy Source from Horizontal Viscosity multiplied by viscous remnant fraction', & - 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) - if (CS%id_KE_horvisc_rem>0) call safe_alloc_ptr(CS%KE_horvisc_rem,isd,ied,jsd,jed,nz) if (.not. adiabatic) then CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, & @@ -2437,11 +2307,7 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. & associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. & associated(CS%KE_horvisc) .or. & - associated(CS%KE_dia) .or. & - associated(CS%PE_to_KE_visc_rem) .or. & - associated(CS%KE_BT_visc_rem) .or. & - associated(CS%KE_CorAdv_visc_rem) .or. & - associated(CS%KE_horvisc_rem)) then + associated(CS%KE_dia)) then call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) endif @@ -2468,26 +2334,6 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) endif - if (associated(CS%KE_CorAdv_visc_rem)) then - call safe_alloc_ptr(ADp%CAu_visc_rem,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%CAv_visc_rem,isd,ied,JsdB,JedB,nz) - endif - if (associated(CS%PE_to_KE_visc_rem)) then - call safe_alloc_ptr(ADp%PFu_visc_rem,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%PFv_visc_rem,isd,ied,JsdB,JedB,nz) - endif - if (associated(CS%KE_BT_visc_rem)) then - call safe_alloc_ptr(ADp%u_accel_BT_visc_rem,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%v_accel_BT_visc_rem,isd,ied,JsdB,JedB,nz) - endif - if (associated(CS%KE_visc_rem)) then - call safe_alloc_ptr(ADp%du_dt_visc_rem,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%dv_dt_visc_rem,isd,ied,JsdB,JedB,nz) - endif - if (associated(CS%KE_horvisc_rem)) then - call safe_alloc_ptr(ADp%diffu_visc_rem,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%diffv_visc_rem,isd,ied,JsdB,JedB,nz) - endif if (associated(CS%KE_dia)) then call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) @@ -2513,16 +2359,11 @@ subroutine MOM_diagnostics_end(CS, ADp) if (associated(CS%KE)) deallocate(CS%KE) if (associated(CS%dKE_dt)) deallocate(CS%dKE_dt) if (associated(CS%PE_to_KE)) deallocate(CS%PE_to_KE) - if (associated(CS%PE_to_KE_visc_rem)) deallocate(CS%PE_to_KE_visc_rem) if (associated(CS%KE_BT)) deallocate(CS%KE_BT) - if (associated(CS%KE_BT_visc_rem)) deallocate(CS%KE_BT_visc_rem) if (associated(CS%KE_Coradv)) deallocate(CS%KE_Coradv) - if (associated(CS%KE_Coradv_visc_rem)) deallocate(CS%KE_Coradv_visc_rem) if (associated(CS%KE_adv)) deallocate(CS%KE_adv) if (associated(CS%KE_visc)) deallocate(CS%KE_visc) - if (associated(CS%KE_visc_rem)) deallocate(CS%KE_visc_rem) if (associated(CS%KE_horvisc)) deallocate(CS%KE_horvisc) - if (associated(CS%KE_horvisc_rem)) deallocate(CS%KE_horvisc_rem) if (associated(CS%KE_dia)) deallocate(CS%KE_dia) if (associated(CS%dv_dt)) deallocate(CS%dv_dt) if (associated(CS%dh_dt)) deallocate(CS%dh_dt) @@ -2541,14 +2382,6 @@ subroutine MOM_diagnostics_end(CS, ADp) if (associated(ADp%dv_dt_dia)) deallocate(ADp%dv_dt_dia) if (associated(ADp%du_other)) deallocate(ADp%du_other) if (associated(ADp%dv_other)) deallocate(ADp%dv_other) - if (associated(ADp%CAu_visc_rem)) deallocate(ADp%CAu_visc_rem) - if (associated(ADp%CAv_visc_rem)) deallocate(ADp%CAv_visc_rem) - if (associated(ADp%PFu_visc_rem)) deallocate(ADp%PFu_visc_rem) - if (associated(ADp%PFv_visc_rem)) deallocate(ADp%PFv_visc_rem) - if (associated(ADp%u_accel_bt_visc_rem)) deallocate(ADp%u_accel_bt_visc_rem) - if (associated(ADp%v_accel_bt_visc_rem)) deallocate(ADp%v_accel_bt_visc_rem) - if (associated(ADp%du_dt_visc_rem)) deallocate(ADp%du_dt_visc_rem) - if (associated(ADp%dv_dt_visc_rem)) deallocate(ADp%dv_dt_visc_rem) if (associated(ADp%diag_hfrac_u)) deallocate(ADp%diag_hfrac_u) if (associated(ADp%diag_hfrac_v)) deallocate(ADp%diag_hfrac_v) From fbdac20c99dfadd1f1dbba8824755c37ec8372a0 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Tue, 3 Aug 2021 20:30:10 -0600 Subject: [PATCH 16/31] Fix small bug --- src/diagnostics/MOM_diagnostics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 88a805fd04..af2b357c3c 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -108,7 +108,7 @@ module MOM_diagnostics KE_adv => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3] KE_visc => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] - KE_dia => NULL(), !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] + KE_dia => NULL() !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] !>@{ Diagnostic IDs integer :: id_u = -1, id_v = -1, id_h = -1 From 86741cd1d3b049303db325d61089263eda90979b Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Tue, 3 Aug 2021 21:51:46 -0600 Subject: [PATCH 17/31] Small style fix --- src/diagnostics/MOM_diagnostics.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index af2b357c3c..2c567ccceb 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1054,7 +1054,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (.not.G%symmetric) then if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. & associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. & - associated(CS%KE_horvisc) .OR. associated(CS%KE_dia)) then + associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) ) then call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East) endif endif @@ -2306,10 +2306,8 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. & associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. & associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. & - associated(CS%KE_horvisc) .or. & - associated(CS%KE_dia)) then + associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) & call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) - endif if (associated(CS%dKE_dt)) then if (.not.associated(CS%du_dt)) then From a40a90f3f0a18b9f2b5bd71672ef835a8f478de1 Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Wed, 4 Aug 2021 15:43:43 -0600 Subject: [PATCH 18/31] Fix allocation problems --- src/core/MOM_dynamics_split_RK2.F90 | 16 +++++++++++----- src/core/MOM_variables.F90 | 2 ++ src/diagnostics/MOM_diagnostics.F90 | 2 +- .../vertical/MOM_vert_friction.F90 | 18 +++++++++--------- 4 files changed, 23 insertions(+), 15 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 52cd711ed9..924e9614d4 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1115,7 +1115,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) PFu_visc_rem(:,:,:) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq - PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%visc_rem_u(I,j,k) + PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%ADp%visc_rem_u(I,j,k) enddo ; enddo ; enddo call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag) deallocate(PFu_visc_rem) @@ -1124,7 +1124,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) PFv_visc_rem(:,:,:) = 0.0 do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%visc_rem_v(i,J,k) + PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%ADp%visc_rem_v(i,J,k) enddo ; enddo ; enddo call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag) deallocate(PFv_visc_rem) @@ -1133,7 +1133,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) CAu_visc_rem(:,:,:) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq - CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%visc_rem_u(I,j,k) + CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%ADp%visc_rem_u(I,j,k) enddo ; enddo ; enddo call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag) deallocate(CAu_visc_rem) @@ -1142,7 +1142,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) CAv_visc_rem(:,:,:) = 0.0 do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%visc_rem_v(i,J,k) + CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%ADp%visc_rem_v(i,J,k) enddo ; enddo ; enddo call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag) deallocate(CAv_visc_rem) @@ -1151,7 +1151,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) u_BT_accel_visc_rem(:,:,:) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq - u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%visc_rem_u(I,j,k) + u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%visc_rem_u(I,j,k) enddo ; enddo ; enddo call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag) deallocate(u_BT_accel_visc_rem) @@ -1673,23 +1673,29 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, & 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, & 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, & 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, & 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, & 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 6f6cf12214..387ef98492 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -195,6 +195,8 @@ module MOM_variables real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points + + integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1 end type accel_diag_ptrs diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 2c567ccceb..e715edcaf0 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -215,7 +215,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & !! [H L2 T-1 ~> m3 s-1 or kg s-1]. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables. - type(accel_diag_ptrs), intent(in) :: ADp !< structure with pointers to + type(accel_diag_ptrs), intent(inout) :: ADp !< structure with pointers to !! accelerations in momentum equation. type(cont_diag_ptrs), intent(in) :: CDp !< structure with pointers to !! terms in continuity equation. diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 83e5a391fa..0182e14f90 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -127,7 +127,7 @@ module MOM_vert_friction ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1 integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1 integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1 - integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1 + ! integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1 ! moved to MOM_variables !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure @@ -528,24 +528,24 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & deallocate(h_dv_dt_visc) endif - if (CS%id_du_dt_visc_rem > 0) then + if (ADp%id_du_dt_visc_rem > 0) then allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) du_dt_visc_rem(:,:,:) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + & (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k) enddo ; enddo ; enddo - call post_data(CS%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag) + call post_data(ADp%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag) deallocate(du_dt_visc_rem) endif - if (CS%id_dv_dt_visc_rem > 0) then + if (ADp%id_dv_dt_visc_rem > 0) then allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) dv_dt_visc_rem(:,:,:) = 0.0 do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + & (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k) enddo ; enddo ; enddo - call post_data(CS%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag) + call post_data(ADp%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag) deallocate(dv_dt_visc_rem) endif @@ -1904,18 +1904,18 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz) endif - CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, & + ADp%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, & 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) - if (CS%id_du_dt_visc_rem > 0) then + if (ADp%id_du_dt_visc_rem > 0) then call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) endif - CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, & + ADp%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, & 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) - if (CS%id_dv_dt_visc_rem > 0) then + if (ADp%id_dv_dt_visc_rem > 0) then call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) From 156d02c122207cb273d286db3dfdac8242cd5fea Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 4 Aug 2021 16:26:24 -0600 Subject: [PATCH 19/31] Make use_drag_rate independent of MEKE_damping Remove MEKE_damping from use_drag_rate so we can use linear dissipation without bottom drag. This is needed for the GEOMETRIC scheme. --- src/parameterizations/lateral/MOM_MEKE.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 850f94cff2..50fd35bae9 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -175,8 +175,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (.not.associated(MEKE)) call MOM_error(FATAL, & "MOM_MEKE: MEKE must be initialized before it is used.") - if ((CS%MEKE_damping > 0.0) .or. (CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) & - .or. CS%visc_drag) then + if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then use_drag_rate = .true. else use_drag_rate = .false. From f31e5e64b3992bea0ff6f4b54abfccb7226f92cc Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Wed, 4 Aug 2021 16:39:51 -0600 Subject: [PATCH 20/31] Remove du_dt_visc_rem diagnostic --- src/core/MOM_variables.F90 | 6 +-- src/diagnostics/MOM_diagnostics.F90 | 6 --- .../vertical/MOM_vert_friction.F90 | 44 ------------------- 3 files changed, 1 insertion(+), 55 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 387ef98492..f7f35ed2d1 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -161,8 +161,6 @@ module MOM_variables ! Each of the following fields has nz layers. real, pointer, dimension(:,:,:) :: & - du_dt => NULL(), & !< Zonal acceleration [L T-2 ~> m s-2] - dv_dt => NULL(), & !< Meridional acceleration [L T-2 ~> m s-2] diffu => NULL(), & !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] diffv => NULL(), & !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2] CAu => NULL(), & !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2] @@ -174,7 +172,7 @@ module MOM_variables du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal mixing [L T-2 ~> m s-2] dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal mixing [L T-2 ~> m s-2] u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2] - v_accel_bt => NULL() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2] + v_accel_bt => NULL() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2] real, pointer, dimension(:,:,:) :: du_other => NULL() !< Zonal velocity changes due to any other processes that are !! not due to any explicit accelerations [L T-1 ~> m s-1]. @@ -195,8 +193,6 @@ module MOM_variables real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points - - integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1 end type accel_diag_ptrs diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index e715edcaf0..cb0f10e6ff 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -290,12 +290,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state) if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state) - do k=1,nz ; do j=js,je ; do I=is-1,ie - ADp%du_dt(I,j,k) = CS%du_dt(I,j,k) - enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je ; do i=is,ie - ADp%dv_dt(i,J,k) = CS%dv_dt(i,J,k) - enddo ; enddo ; enddo !! Diagnostics for terms multiplied by fractional thicknesses diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 0182e14f90..081179eb41 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -127,7 +127,6 @@ module MOM_vert_friction ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1 integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1 integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1 - ! integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1 ! moved to MOM_variables !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure @@ -218,11 +217,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u) - ! [L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v) - ! [L T-2 ~> m2 s-2] - logical :: do_i(SZIB_(G)) logical :: DoStokesMixing @@ -528,27 +522,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & deallocate(h_dv_dt_visc) endif - if (ADp%id_du_dt_visc_rem > 0) then - allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) - du_dt_visc_rem(:,:,:) = 0.0 - do k=1,nz ; do j=js,je ; do I=Isq,Ieq - du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + & - (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k) - enddo ; enddo ; enddo - call post_data(ADp%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag) - deallocate(du_dt_visc_rem) - endif - if (ADp%id_dv_dt_visc_rem > 0) then - allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) - dv_dt_visc_rem(:,:,:) = 0.0 - do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + & - (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k) - enddo ; enddo ; enddo - call post_data(ADp%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag) - deallocate(dv_dt_visc_rem) - endif - end subroutine vertvisc !> Calculate the fraction of momentum originally in a layer that remains @@ -1904,23 +1877,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz) endif - ADp%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, & - 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) - if (ADp%id_du_dt_visc_rem > 0) then - call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) - endif - ADp%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, & - 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) - if (ADp%id_dv_dt_visc_rem > 0) then - call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz) - call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) - call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) - endif - if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) & call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp) From ca0fb8950851b235a18e83ae0a254b6670abb3bc Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Wed, 4 Aug 2021 17:38:57 -0600 Subject: [PATCH 21/31] Fix small inconsistency: CS%visc_rem_v --> CS%ADp%visc_rem_v --- src/core/MOM_dynamics_split_RK2.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 924e9614d4..ebe53fc908 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1160,7 +1160,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) v_BT_accel_visc_rem(:,:,:) = 0.0 do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%visc_rem_v(i,J,k) + v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%visc_rem_v(i,J,k) enddo ; enddo ; enddo call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag) deallocate(v_BT_accel_visc_rem) From a874d6df79e79b84a7706dd7b77a10a97f4c35ea Mon Sep 17 00:00:00 2001 From: NoraLoose Date: Wed, 4 Aug 2021 17:40:44 -0600 Subject: [PATCH 22/31] Don't need inout for ADp anymore --- src/diagnostics/MOM_diagnostics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index cb0f10e6ff..f874d08a12 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -215,7 +215,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & !! [H L2 T-1 ~> m3 s-1 or kg s-1]. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables. - type(accel_diag_ptrs), intent(inout) :: ADp !< structure with pointers to + type(accel_diag_ptrs), intent(in) :: ADp !< structure with pointers to !! accelerations in momentum equation. type(cont_diag_ptrs), intent(in) :: CDp !< structure with pointers to !! terms in continuity equation. From 03d034567b85cdc842abc56ead68e030b4759904 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 18 Aug 2021 12:23:32 -0400 Subject: [PATCH 23/31] Testing: Replace bc with awk; framework support This patch adds three minor updates to the build and test systems: 1. A framework selection bug was fixed. The variable name was missing a $ token. 2. the FRAMEWORK variable was added to .testing/Makefile to facilitate testing of both FMS frameworks. 3. Restart half-periods no longer use `bc` to compute the duration in seconds; we now use awk, which is more commonly available. --- .testing/Makefile | 16 +++++++++++----- ac/configure.ac | 2 +- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index 06b29dc690..de7038eeff 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -114,6 +114,9 @@ CONFIGS := $(wildcard tc*) TESTS = grids layouts restarts nans dims openmps rotations DIMS = t l h z q r +# Set the framework +FRAMEWORK ?= fms1 + # REPRO tests enable reproducibility with optimization, and often do not match # the DEBUG results in older GCCs and vendor compilers, so we can optionally # disable them. @@ -267,7 +270,7 @@ build/%/MOM6: build/%/Makefile build/%/Makefile: ../ac/configure ../ac/Makefile.in $(DEPS)/lib/libFMS.a $(MKMF) $(LIST_PATHS) mkdir -p $(@D) cd $(@D) \ - && $(MOM_ENV) ../../../ac/configure $(MOM_ACFLAGS) \ + && $(MOM_ENV) ../../../ac/configure $(MOM_ACFLAGS) --with-framework=$(FRAMEWORK) \ || (cat config.log && false) @@ -549,7 +552,11 @@ $(eval $(call STAT_RULE,dim.q,symmetric,,Q_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.r,symmetric,,R_RESCALE_POWER=11,,1)) -# Restart tests require significant preprocessing, and are handled separately. +# Generate the half-period input namelist as follows: +# 1. Fetch DAYMAX and TIMEUNIT from MOM_input +# 2. Convert DAYMAX from TIMEUNIT to seconds +# 3. Apply seconds to `ocean_solo_nml` inside input.nml. +# NOTE: Assumes that runtime set by DAYMAX, will fail if set by input.nml work/%/restart/ocean.stats: build/symmetric/MOM6 $(VENV_PATH) rm -rf $(@D) mkdir -p $(@D) @@ -562,14 +569,13 @@ work/%/restart/ocean.stats: build/symmetric/MOM6 $(VENV_PATH) cd work/$*/restart; \ fi mkdir -p $(@D)/RESTART - # Generate the half-period input namelist - # TODO: Assumes that runtime set by DAYMAX, will fail if set by input.nml + # Set the half-period cd $(@D) \ && daymax=$$(grep DAYMAX MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ && timeunit=$$(grep TIMEUNIT MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ && if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \ && printf -v timeunit_int "%.f" "$${timeunit}" \ - && halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \ + && halfperiod=$$(awk -v t=$${daymax} -v dt=$${timeunit} 'BEGIN {printf "%.f", 0.5*t*dt}') \ && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml # Remove any previous archived output rm -f results/$*/std.restart{1,2}.{out,err} diff --git a/ac/configure.ac b/ac/configure.ac index 9cb7147846..3d1af81b05 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -65,7 +65,7 @@ AS_IF([test "x$with_driver" != "x"], MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1 AC_ARG_WITH([framework], AS_HELP_STRING([--with-framework=fms1|fms2], [Select the model framework])) -AS_CASE([with_framework], +AS_CASE(["$with_framework"], [fms1], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1], [fms2], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1] From 8ee65a6ff45efb6e6c5deb780e0b12f87361e6a9 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 23 Aug 2021 11:53:49 -0600 Subject: [PATCH 24/31] Fix units of *_visc_rem terms --- src/core/MOM_dynamics_split_RK2.F90 | 24 +++++++++---------- .../lateral/MOM_hor_visc.F90 | 12 +++++----- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 96a0a5f92f..a168fe1319 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1677,30 +1677,30 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, & - 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, & - 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, & - 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, & - 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, & - 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, & - 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 592fca4d24..61127c745f 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -281,8 +281,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: h_diffv ! h x diffv [H L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m2 s-2] + real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m s-2] + real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m s-2] real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] @@ -2471,15 +2471,15 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) endif CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, & - 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke) endif CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, & - 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', & - conversion=GV%H_to_m*US%L_T2_to_m_s2) + 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) endif From dc205c3a29a1d9304010582d329f03c87f5fec33 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 26 Aug 2021 16:46:43 -0400 Subject: [PATCH 25/31] (*)Fix discretization issues with USE_GME=True Corrected the inconsistently staggered expressions for the total depths used to scale away the GM+E backscatter scheme that is enabled with USE_GME. As a part of these changes, the nominal bathymetric depths have been replaced by the sum of the layer thicknesses, and parentheses were added to some of the GME expressions so that they will be rotationally invariant. Although there are no changes to the answers in the MOM6-examples test suite, answers will change in any cases where USE_GME is true, and GME_H0 is larger than the minimum ocean depth. This commit will address MOM6 issue #1466, which can now be closed. --- .../lateral/MOM_hor_visc.F90 | 68 ++++++++++++------- 1 file changed, 44 insertions(+), 24 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index b4f857dec4..cf1712afc6 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -276,6 +276,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [T-2 ~> s-2] grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2] + GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] + htot, & ! The total thickness of all layers [Z ~> m] boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2] @@ -301,6 +303,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] ! This form guarantees that hq/hu < 4. grad_vel_mag_bt_q, & ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] + GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim] boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: & @@ -338,6 +341,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity ! points where masks are applied [H ~> m or kg m-2]. + real :: h_arith_q ! The arithmetic mean total thickness at q points [Z ~> m] + real :: h_harm_q ! The harmonic mean total thickness at q points [Z ~> m] + real :: I_hq ! The inverse of the arithmetic mean total thickness at q points [Z-1 ~> m-1] + real :: I_GME_h0 ! The inverse of GME tapering scale [Z-1 ~> m-1] real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2] real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6] real :: h_min ! Minimum h at the 4 neighboring velocity points [H ~> m] @@ -501,6 +508,35 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) enddo ; enddo + do j=js-1,je+1 ; do i=is-1,ie+1 + htot(i,j) = 0.0 + enddo ; enddo + do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + htot(i,j) = htot(i,j) + GV%H_to_Z*h(i,j,k) + enddo ; enddo ; enddo + + I_GME_h0 = 1.0 / CS%GME_h0 + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (grad_vel_mag_bt_h(i,j)>0) then + GME_effic_h(i,j) = CS%GME_efficiency * boundary_mask_h(i,j) * & + (MIN(htot(i,j) * I_GME_h0, 1.0)**2) + else + GME_effic_h(i,j) = 0.0 + endif + enddo ; enddo + + do J=js-1,Jeq ; do I=is-1,Ieq + if (grad_vel_mag_bt_q(I,J)>0) then + h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1))) + I_hq = 1.0 / h_arith_q + h_harm_q = 0.25 * h_arith_q * ((htot(i,j)*I_hq + htot(i+1,j+1)*I_hq) + & + (htot(i+1,j)*I_hq + htot(i,j+1)*I_hq)) + GME_effic_q(I,J) = CS%GME_efficiency * boundary_mask_q(I,J) * (MIN(h_harm_q * I_GME_h0, 1.0)**2) + else + GME_effic_q(I,J) = 0.0 + endif + enddo ; enddo + endif ! use_GME !$OMP parallel do default(none) & @@ -509,7 +545,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, & !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & !$OMP use_MEKE_Ku, use_MEKE_Au, boundary_mask_h, boundary_mask_q, & - !$OMP backscat_subround, GME_coeff_limiter, & + !$OMP backscat_subround, GME_coeff_limiter, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & @@ -1388,15 +1424,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call pass_vector(KH_u_GME, KH_v_GME, G%Domain) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (grad_vel_mag_bt_h(i,j)>0) then - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & - (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) - else - GME_coeff = 0.0 - endif - - ! apply mask - GME_coeff = GME_coeff * boundary_mask_h(i,j) + GME_coeff = GME_effic_h(i,j) * 0.25 * & + ((KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff @@ -1405,17 +1434,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq - if (grad_vel_mag_bt_q(I,J)>0) then - !### This expression is not rotationally invariant - bathyT is to the SW of q points, - ! and it needs parentheses in the sum of the 4 diffusivities. - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & - (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) - else - GME_coeff = 0.0 - endif - - ! apply mask - GME_coeff = GME_coeff * boundary_mask_q(I,J) + GME_coeff = GME_effic_q(I,J) * 0.25 * & + ((KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff @@ -1424,8 +1444,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo ! Applying GME diagonal term. This is linear and the arguments can be rescaled. - call smooth_GME(CS,G,GME_flux_h=str_xx_GME) - call smooth_GME(CS,G,GME_flux_q=str_xy_GME) + call smooth_GME(CS, G, GME_flux_h=str_xx_GME) + call smooth_GME(CS, G, GME_flux_q=str_xy_GME) do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j)) @@ -1446,7 +1466,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif - else ! use_GME + else ! .not. use_GME do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo From 12569050c3ab746387b1ecbdc7413f571d220ef6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 26 Aug 2021 17:02:35 -0400 Subject: [PATCH 26/31] (*)Use depths in DOME and dye tracer initialization Initialize the horizontal stripes of tracers in the DOME and dye examples based on the depth from the (flat) sea surface rather than the height above the seafloor. These are mathematically equivalent, but there could be small changes in the tracer distributions. Because these changes only impact passive tracers, all of the physical solutions are bitwise identical. --- src/tracer/DOME_tracer.F90 | 16 +++++++++------- src/tracer/dye_example.F90 | 18 ++++++++++-------- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index c20eda7745..44421c7387 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -75,7 +75,7 @@ function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) character(len=48) :: flux_units ! The units for tracer fluxes, usually ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. character(len=200) :: inputdir - real, pointer :: tr_ptr(:,:,:) => NULL() + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to the tracer field logical :: register_DOME_tracer integer :: isd, ied, jsd, jed, nz, m isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke @@ -166,13 +166,15 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & character(len=48) :: units ! The dimensions of the variable. character(len=48) :: flux_units ! The units for tracer fluxes, usually ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. - real, pointer :: tr_ptr(:,:,:) => NULL() + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to the tracer field real :: PI ! 3.1415926... calculated as 4*atan(1) real :: tr_y ! Initial zonally uniform tracer concentrations. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. - real :: e(SZK_(GV)+1), e_top, e_bot ! Heights [Z ~> m]. - real :: d_tr ! A change in tracer concentraions, in tracer units. + real :: e(SZK_(GV)+1) ! Interface heights relative to the sea surface (negative down) [Z ~> m] + real :: e_top ! Height of the top of the tracer band relative to the sea surface [Z ~> m] + real :: e_bot ! Height of the bottom of the tracer band relative to the sea surface [Z ~> m] + real :: d_tr ! A change in tracer concentrations, in tracer units. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB @@ -215,9 +217,9 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & if (NTR > 7) then do j=js,je ; do i=is,ie - e(nz+1) = -G%bathyT(i,j) - do k=nz,1,-1 - e(K) = e(K+1) + h(i,j,k)*GV%H_to_Z + e(1) = 0.0 + do k=1,nz + e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z do m=7,NTR e_top = (-600.0*real(m-1) + 3000.0) * US%m_to_Z e_bot = (-600.0*real(m-1) + 2700.0) * US%m_to_Z diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index 2919f2d95f..2bf3cd94ed 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -203,9 +203,10 @@ subroutine initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS, sponge_C character(len=48) :: units ! The dimensions of the variable. character(len=48) :: flux_units ! The units for age tracer fluxes, either ! years m3 s-1 or years kg s-1. + real :: z_bot ! Height of the bottom of the layer relative to the sea surface [Z ~> m] + real :: z_center ! Height of the center of the layer relative to the sea surface [Z ~> m] logical :: OK integer :: i, j, k, m - real :: z_bot, z_center if (.not.associated(CS)) return if (CS%ntr < 1) return @@ -221,14 +222,14 @@ subroutine initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS, sponge_C CS%dye_source_minlat(m)=G%geoLatT(i,j) .and. & G%mask2dT(i,j) > 0.0 ) then - z_bot = -G%bathyT(i,j) - do k = GV%ke, 1, -1 + z_bot = 0.0 + do k = 1, GV%ke + z_bot = z_bot - h(i,j,k)*GV%H_to_Z z_center = z_bot + 0.5*h(i,j,k)*GV%H_to_Z if ( z_center > -CS%dye_source_maxdepth(m) .and. & z_center < -CS%dye_source_mindepth(m) ) then CS%tr(i,j,k,m) = 1.0 endif - z_bot = z_bot + h(i,j,k)*GV%H_to_Z enddo endif enddo ; enddo @@ -273,9 +274,10 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US real :: sfc_val ! The surface value for the tracers. real :: Isecs_per_year ! The number of seconds in a year. real :: year ! The time in years. + real :: z_bot ! Height of the bottom of the layer relative to the sea surface [Z ~> m] + real :: z_center ! Height of the center of the layer relative to the sea surface [Z ~> m] integer :: secs, days ! Integer components of the time type. integer :: i, j, k, is, ie, js, je, nz, m - real :: z_bot, z_center is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -305,14 +307,14 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US CS%dye_source_minlat(m)=G%geoLatT(i,j) .and. & G%mask2dT(i,j) > 0.0 ) then - z_bot = -G%bathyT(i,j) - do k=nz,1,-1 + z_bot = 0.0 + do k=1,nz + z_bot = z_bot - h_new(i,j,k)*GV%H_to_Z z_center = z_bot + 0.5*h_new(i,j,k)*GV%H_to_Z if ( z_center > -CS%dye_source_maxdepth(m) .and. & z_center < -CS%dye_source_mindepth(m) ) then CS%tr(i,j,k,m) = 1.0 endif - z_bot = z_bot + h_new(i,j,k)*GV%H_to_Z enddo endif enddo ; enddo From 8cd43302e7810ed93e1e415167954c97ff965084 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 26 Aug 2021 21:18:10 -0400 Subject: [PATCH 27/31] (*+)New internal tide bathymetric parameters Replaced the nominal bathymetric depth with the sum of the layer thicknesses in the calculation of the quadratic drag acting on the internal wave field. Also replaced the hard-coded 1 m scale at which the drag is reduced with a new runtime parameter, INTERNAL_TIDE_DRAG_MIN_DEPTH, which is only logged if the internal tides are used and INTERNAL_TIDE_QUAD_DRAG=True, and made the maximum scale of the topographic roughness as a fraction of the bottom depth into the new runtime parameter INTERNAL_TIDE_ROUGHNESS_FRAC, whose default, 0.1, is the same as the previous hard-coded value. Several of the comments in the MOM_internal_tides module were also clarified or standardized. All answers in the MOM-examples test suite are bitwise identical, but answers could change and there could be changes in the MOM_parameter_doc files in some cases that use the internal tides module. --- .../lateral/MOM_internal_tides.F90 | 57 ++++++++++++------- 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 4d66471408..1bff6287f4 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -96,6 +96,9 @@ module MOM_internal_tides real :: decay_rate !< A constant rate at which internal tide energy is !! lost to the interior ocean internal wave field. real :: cdrag !< The bottom drag coefficient [nondim]. + real :: drag_min_depth !< The minimum total ocean thickness that will be used in the denominator + !! of the quadratic drag terms for internal tides when + !! INTERNAL_TIDE_QUAD_DRAG is true [Z ~> m] logical :: apply_background_drag !< If true, apply a drag due to background processes as a sink. logical :: apply_bottom_drag @@ -185,6 +188,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & tot_En, & ! energy summed over angles, modes, frequencies [R Z3 T-2 ~> J m-2] tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_Froude_loss, tot_allprocesses_loss, & ! energy loss rates summed over angle, freq, and mode [R Z3 T-3 ~> W m-2] + htot, & ! The vertical sum of the layer thicknesses [H ~> m or kg m-2] drag_scale, & ! bottom drag scale [T-1 ~> s-1] itidal_loss_mode, allprocesses_loss_mode ! energy loss rates for a given mode and frequency (summed over angles) [R Z3 T-3 ~> W m-2] @@ -200,7 +204,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & real :: En_initial, Delta_E_check ! Energies for debugging [R Z3 T-2 ~> J m-2] real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! Energy losses for debugging [R Z3 T-3 ~> W m-2] character(len=160) :: mesg ! The text of an error message - integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm + integer :: a, m, fr, i, j, k, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging) type(group_pass_type), save :: pass_test, pass_En type(time_type) :: time_end @@ -360,9 +364,12 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Extract the energy for mixing due to bottom drag------------------------------- if (CS%apply_bottom_drag) then + do j=jsd,jed ; do i=isd,ied ; htot(i,j) = 0.0 ; enddo ; enddo + do k=1,GV%ke ; do j=jsd,jed ; do i=isd,ied + htot(i,j) = htot(i,j) + h(i,j,k) + enddo ; enddo ; enddo do j=jsd,jed ; do i=isd,ied - ! Note the 1 m dimensional scale here. Should this be a parameter? - I_D_here = 1.0 / (max(G%bathyT(i,j), 1.0*US%m_to_Z)) + I_D_here = 1.0 / (max(GV%H_to_Z*htot(i,j), CS%drag_min_depth)) drag_scale(i,j) = CS%cdrag * sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + & tot_En(i,j) * I_rho0 * I_D_here)) * I_D_here enddo ; enddo @@ -2143,20 +2150,20 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) ! Local variables real :: Angle_size ! size of wedges, rad real, allocatable :: angles(:) ! orientations of wedge centers, rad - real, allocatable, dimension(:,:) :: h2 ! topographic roughness scale, m^2 - real :: kappa_itides, kappa_h2_factor - ! characteristic topographic wave number - ! and a scaling factor - real, allocatable :: ridge_temp(:,:) - ! array for temporary storage of flags + real, dimension(:,:), allocatable :: h2 ! topographic roughness scale squared [Z2 ~> m2] + real :: kappa_itides ! characteristic topographic wave number [L-1 ~> m-1] + real, dimension(:,:), allocatable :: ridge_temp ! array for temporary storage of flags ! of cells with double-reflecting ridges - logical :: use_int_tides, use_temperature - real :: period_1 ! The period of the gravest modeled mode [T ~> s] + logical :: use_int_tides, use_temperature + real :: kappa_h2_factor ! A roughness scaling factor [nondim] + real :: RMS_roughness_frac ! The maximum RMS topographic roughness as a fraction of the + ! nominal ocean depth, or a negative value for no limit [nondim] + real :: period_1 ! The period of the gravest modeled mode [T ~> s] integer :: num_angle, num_freq, num_mode, m, fr integer :: isd, ied, jsd, jed, a, id_ang, i, j type(axes_grp) :: axes_ang ! This include declares and sets the variable "version". -#include "version_variable.h" +# include "version_variable.h" character(len=40) :: mdl = "MOM_internal_tides" ! This module's name. character(len=16), dimension(8) :: freq_name character(len=40) :: var_name @@ -2280,16 +2287,20 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "1st-order upwind advection. This scheme is highly "//& "continuity solver. This scheme is highly "//& "diffusive but may be useful for debugging.", default=.false.) - call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", & - CS%apply_background_drag, "If true, the internal tide "//& - "ray-tracing advection uses a background drag term as a sink.",& - default=.false.) + call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", CS%apply_background_drag, & + "If true, the internal tide ray-tracing advection uses a background drag "//& + "term as a sink.", default=.false.) call get_param(param_file, mdl, "INTERNAL_TIDE_QUAD_DRAG", CS%apply_bottom_drag, & "If true, the internal tide ray-tracing advection uses "//& "a quadratic bottom drag term as a sink.", default=.false.) call get_param(param_file, mdl, "INTERNAL_TIDE_WAVE_DRAG", CS%apply_wave_drag, & "If true, apply scattering due to small-scale roughness as a sink.", & default=.false.) + call get_param(param_file, mdl, "INTERNAL_TIDE_DRAG_MIN_DEPTH", CS%drag_min_depth, & + "The minimum total ocean thickness that will be used in the denominator "//& + "of the quadratic drag terms for internal tides.", & + units="m", default=1.0, scale=US%m_to_Z, do_not_log=.not.CS%apply_bottom_drag) + CS%drag_min_depth = MAX(CS%drag_min_depth, GV%H_subroundoff * GV%H_to_Z) call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", CS%apply_Froude_drag, & "If true, apply wave breaking as a sink.", & default=.false.) @@ -2340,10 +2351,18 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) fail_if_missing=.true.) filename = trim(CS%inputdir) // trim(h2_file) call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename) - call MOM_read_data(filename, 'h2', h2, G%domain, scale=US%m_to_Z) + call get_param(param_file, mdl, "INTERNAL_TIDE_ROUGHNESS_FRAC", RMS_roughness_frac, & + "The maximum RMS topographic roughness as a fraction of the nominal ocean depth, "//& + "or a negative value for no limit.", units="nondim", default=0.1) + + call MOM_read_data(filename, 'h2', h2, G%domain, scale=US%m_to_Z**2) do j=G%jsc,G%jec ; do i=G%isc,G%iec - ! Restrict rms topo to 10 percent of column depth. - h2(i,j) = min(0.01*(G%bathyT(i,j))**2, h2(i,j)) + ! Restrict RMS topographic roughness to a fraction (10 percent by default) of the column depth. + if (RMS_roughness_frac >= 0.0) then + h2(i,j) = max(min((RMS_roughness_frac*G%bathyT(i,j))**2, h2(i,j)), 0.0) + else + h2(i,j) = max(h2(i,j), 0.0) + endif ! Compute the fixed part; units are [R L-2 Z3 ~> kg m-2] here ! will be multiplied by N and the squared near-bottom velocity to get into [R Z3 T-3 ~> W m-2] CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*GV%Rho0 * US%L_to_Z*kappa_itides * h2(i,j) From 6979c298419a3a724edacfda506f378a6fab96f1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 26 Aug 2021 21:36:49 -0400 Subject: [PATCH 28/31] +(*)Add option for MEKE to calculate total depth Added the new run-time option, MEKE_FIXED_TOTAL_DEPTH, to use the actual total thickness instead the nominal depth in bathyT in several of the MEKE calculations. Also simplified and corrected a minor dimensional inconsistency in some expressions that effectively set masks when estimating the interface height diffusivities at tracer points for use in MEKE. By default, the answers in the MOM-examples test cases are bitwise identical, but answers could change in some cases due to the proper thickness weighting in the calculation of MEKE%Kh_diff. There are new entries in some MOM_parameter_doc files, so a minor spelling error correction in a MOM_parameter_doc entry was also included in this PR. --- src/parameterizations/lateral/MOM_MEKE.F90 | 34 ++++++++++++------- .../lateral/MOM_thickness_diffuse.F90 | 33 ++++++++++-------- .../vertical/MOM_ALE_sponge.F90 | 4 +-- 3 files changed, 43 insertions(+), 28 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 8779968bcd..8206fd9717 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -82,7 +82,9 @@ module MOM_MEKE real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. - + logical :: fixed_total_depth !< If true, use the nominal bathymetric depth as the estimate of + !! the time-varying ocean depth. Otherwise base the depth on the total + !! ocean mass per unit area. logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. logical :: debug !< If true, write out checksums of data for debugging @@ -283,12 +285,17 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo enddo - !$OMP parallel do default(shared) - do j=js-1,je+1 ; do i=is-1,ie+1 - depth_tot(i,j) = G%bathyT(i,j) - !### Try changing this to: - ! depth_tot(i,j) = mass(i,j) * I_Rho0 - enddo ; enddo + if (CS%fixed_total_depth) then + !$OMP parallel do default(shared) + do j=js-1,je+1 ; do i=is-1,ie+1 + depth_tot(i,j) = G%bathyT(i,j) + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js-1,je+1 ; do i=is-1,ie+1 + depth_tot(i,j) = mass(i,j) * I_Rho0 + enddo ; enddo + endif if (CS%initialize) then call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot) @@ -711,8 +718,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then beta_topo_x = 0. ; beta_topo_y = 0. else - !### Consider different combinations of these estimates of topographic beta, and the use - ! of the water column thickness instead of the bathymetric depth. + !### Consider different combinations of these estimates of topographic beta. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) & / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) & @@ -887,14 +893,13 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, & FatH = 0.25* ( ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) ) + & ( G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1) ) ) ! Coriolis parameter at h points - ! If bathyT is zero, then a division by zero FPE will be raised. In this + ! If depth_tot is zero, then a division by zero FPE will be raised. In this ! case, we apply Adcroft's rule of reciprocals and set the term to zero. ! Since zero-bathymetry cells are masked, this should not affect values. if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then beta_topo_x = 0. ; beta_topo_y = 0. else - !### Consider different combinations of these estimates of topographic beta, and the use - ! of the water column thickness instead of the bathymetric depth. + !### Consider different combinations of these estimates of topographic beta. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) & / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) & @@ -1167,6 +1172,11 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) "If positive, is a fixed length contribution to the expression "//& "for mixing length used in MEKE-derived diffusivity.", & units="m", default=0.0, scale=US%m_to_L) + call get_param(param_file, mdl, "MEKE_FIXED_TOTAL_DEPTH", CS%fixed_total_depth, & + "If true, use the nominal bathymetric depth as the estimate of the "//& + "time-varying ocean depth. Otherwise base the depth on the total ocean mass"//& + "per unit area.", default=.true.) + call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", CS%aDeform, & "If positive, is a coefficient weighting the deformation scale "//& "in the expression for mixing length used in MEKE-derived diffusivity.", & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 3b3d72576c..0c8f25b42e 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -145,6 +145,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp KH_u_CFL ! The maximum stable interface height diffusivity at u grid points [L2 T-1 ~> m2 s-1] real, dimension(SZI_(G), SZJB_(G)) :: & KH_v_CFL ! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1] + real, dimension(SZI_(G), SZJ_(G)) :: & + htot ! The sum of the total layer thicknesses [H ~> m or kg m-2] real :: Khth_Loc_u(SZIB_(G), SZJ_(G)) real :: Khth_Loc_v(SZI_(G), SZJB_(G)) real :: Khth_Loc(SZIB_(G), SZJB_(G)) ! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1] @@ -479,38 +481,41 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp ! thicknesses across u and v faces, converted to 0/1 mask ! layer average of the interface diffusivities KH_u and KH_v do j=js,je ; do I=is-1,ie - hu(I,j) = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect) - if (hu(I,j) /= 0.0) hu(I,j) = 1.0 - !### The same result would be accomplished with the following without a division: - ! hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0 + ! This expression uses harmonic mean thicknesses: + ! hu(I,j) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect) + ! This expression is a 0/1 mask based on depths where there are thick layers: + hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0 KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1)) enddo ; enddo do J=js-1,je ; do i=is,ie - hv(i,J) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect) - if (hv(i,J) /= 0.0) hv(i,J) = 1.0 - !### The same result would be accomplished with the following without a division: - ! hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0 + ! This expression uses harmonic mean thicknesses: + ! hv(i,J) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect) + ! This expression is a 0/1 mask based on depths where there are thick layers: + hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0 KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1)) enddo ; enddo - ! diagnose diffusivity at T-point - !### Because hu and hv are nondimensional here, the denominator is dimensionally inconsistent. + ! diagnose diffusivity at T-points do j=js,je ; do i=is,ie - Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j)) & - +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) & - / (hu(I-1,j)+hu(I,j)+hv(i,J-1)+hv(i,J)+h_neglect) + Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j) + hu(I,j)*KH_u_lay(I,j)) + & + (hv(i,J-1)*KH_v_lay(i,J-1) + hv(i,J)*KH_v_lay(i,J))) / & + ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + 1.0e-20) + ! Use this denominator instead if hu and hv are actual thicknesses rather than a 0/1 mask: + ! ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + h_neglect) enddo ; enddo enddo if (CS%Use_KH_in_MEKE) then MEKE%Kh_diff(:,:) = 0.0 + htot(:,:) = 0.0 do k=1,nz do j=js,je ; do i=is,ie MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) + Kh_t(i,j,k) * h(i,j,k) + htot(i,j) = htot(i,j) + h(i,j,k) enddo ; enddo enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = GV%H_to_Z * MEKE%Kh_diff(i,j) / MAX(1.0*US%m_to_Z, G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0*GV%m_to_H, htot(i,j)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 419b012387..be1d265620 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -218,7 +218,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, & - "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "If true, use the order of arithmetic for horizontal regridding that recovers "//& "the answers from the end of 2018. Otherwise, use rotationally symmetric "//& "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, & @@ -478,7 +478,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, & - "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "If true, use the order of arithmetic for horizontal regridding that recovers "//& "the answers from the end of 2018 and retain a bug in the 3-dimensional mask "//& "returned in certain cases. Otherwise, use rotationally symmetric "//& "forms of the same expressions and initialize the mask properly.", & From f347b10ce4f484eb49edc9cb7595adf6f4defd1d Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 27 Aug 2021 11:13:32 -0400 Subject: [PATCH 29/31] Testing: Disable 2018 answer flags The TC tests are still using 2018 answer reproducibility flags, which fail when aggressive initialization of allocatable arrays is enabled. In at least one case, the `mask_z` variable of the ALE sponge is incompletely initialized when `DEFAULT_2018_ANSWERS` is set. When unset, the array is fully initialized. This patch disables the 2018 answer flags and restores reproducibility on such platforms. --- .testing/tc0/MOM_input | 1 - .testing/tc1/MOM_input | 1 - .testing/tc2/MOM_input | 1 - .testing/tc3/MOM_input | 1 - .testing/tc4/MOM_input | 1 - 5 files changed, 5 deletions(-) diff --git a/.testing/tc0/MOM_input b/.testing/tc0/MOM_input index ff64c55803..e4d1694e72 100644 --- a/.testing/tc0/MOM_input +++ b/.testing/tc0/MOM_input @@ -230,7 +230,6 @@ ENERGYSAVEDAYS = 1.0 DIAG_AS_CHKSUM = True DEBUG = True -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True USE_GM_WORK_BUG = True ! [Boolean] default = True FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False diff --git a/.testing/tc1/MOM_input b/.testing/tc1/MOM_input index 68674f7a86..151c093ff9 100644 --- a/.testing/tc1/MOM_input +++ b/.testing/tc1/MOM_input @@ -575,7 +575,6 @@ ENERGYSAVEDAYS = 0.125 ! [days] default = 3600.0 DIAG_AS_CHKSUM = True DEBUG = True USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True INTERPOLATE_RES_FN = True ! [Boolean] default = True GILL_EQUATORIAL_LD = False ! [Boolean] default = False diff --git a/.testing/tc2/MOM_input b/.testing/tc2/MOM_input index 1818390192..ca84d1c382 100644 --- a/.testing/tc2/MOM_input +++ b/.testing/tc2/MOM_input @@ -610,7 +610,6 @@ DIAG_AS_CHKSUM = True DEBUG = True USE_GM_WORK_BUG = False USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True REMAP_UV_USING_OLD_ALG = True ! [Boolean] default = True USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = False diff --git a/.testing/tc3/MOM_input b/.testing/tc3/MOM_input index 9112898b4c..a034960d1e 100644 --- a/.testing/tc3/MOM_input +++ b/.testing/tc3/MOM_input @@ -469,7 +469,6 @@ ENERGYSAVEDAYS = 3.0 ! [hours] default = 1.44E+04 ! energies of the run and other globally summed diagnostics. DIAG_AS_CHKSUM = True DEBUG = True -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True OBC_RADIATION_MAX = 10.0 ! [nondim] default = 10.0 GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True USE_GM_WORK_BUG = True ! [Boolean] default = True diff --git a/.testing/tc4/MOM_input b/.testing/tc4/MOM_input index 04598a9dc9..e33bf40bf6 100644 --- a/.testing/tc4/MOM_input +++ b/.testing/tc4/MOM_input @@ -409,7 +409,6 @@ DIAG_AS_CHKSUM = True DEBUG = True USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True INTERPOLATE_RES_FN = True ! [Boolean] default = True GILL_EQUATORIAL_LD = False ! [Boolean] default = False From 868c11a2a91603b6673cf31596457507c19f14e7 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Aug 2021 10:40:58 -0600 Subject: [PATCH 30/31] Always initialize drag_rate This patch now always initializes drag_rate, either to zero or to the bottom drag formula. It also changes the logic to always execute the second Strang splitting calculation. Removed unused code, including drag_rate_J15 and some commented code. --- src/parameterizations/lateral/MOM_MEKE.F90 | 41 ++++++++-------------- 1 file changed, 14 insertions(+), 27 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 50fd35bae9..7b92704651 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -130,8 +130,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1]. drag_rate_visc, & ! Near-bottom velocity contribution to bottom dratg [L T-1 ~> m s-1] drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1]. - drag_rate_J15, & ! The MEKE spindown timescale due to bottom drag with the Jansen 2015 scheme. - ! Unfortunately, as written the units seem inconsistent. [T-1 ~> s-1]. del2MEKE, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2]. del4MEKE, & ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2]. LmixScale, & ! Eddy mixing length [L ~> m]. @@ -230,12 +228,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo endif - if (CS%MEKE_Cd_scale == 0.0 .and. .not. CS%visc_drag) then - !$OMP parallel do default(shared) private(ldamping) - do j=js,je ; do i=is,ie - drag_rate(i,j) = 0. ; drag_rate_J15(i,j) = 0. - enddo ; enddo - endif ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow if (CS%visc_drag) then @@ -339,12 +331,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif - ! Increase EKE by a full time-steps worth of source - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j) - enddo ; enddo - if (use_drag_rate) then ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies) !$OMP parallel do default(shared) @@ -352,6 +338,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h drag_rate(i,j) = (US%L_to_Z*Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + drag_rate(i,j) = 0. + enddo ; enddo endif ! First stage of Strang splitting @@ -520,23 +511,19 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h drag_rate(i,j) = (US%L_to_Z*Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) ldamping = 0. - ! notice that the above line ensures a damping only if MEKE is positive, - ! while leaving MEKE unchanged if it is negative - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) - MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) - enddo ; enddo endif + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) ldamping = 0. + ! notice that the above line ensures a damping only if MEKE is positive, + ! while leaving MEKE unchanged if it is negative + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) + MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) + enddo ; enddo endif endif ! MEKE_KH>=0 - ! do j=js,je ; do i=is,ie - ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0) - ! enddo ; enddo - call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_MEKE, G%Domain) call cpu_clock_end(CS%id_clock_pass) From ba6f56d5f002ebf2093ecb8e11b615acd65d827b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Aug 2021 13:33:50 -0600 Subject: [PATCH 31/31] Add code that was deleted unintentionally --- src/parameterizations/lateral/MOM_MEKE.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 66f75c82e3..6ce0e58471 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -331,6 +331,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif + ! Increase EKE by a full time-steps worth of source + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j) + enddo ; enddo + if (use_drag_rate) then ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies) !$OMP parallel do default(shared)