From 916700dac14d4d48f3c6f5536d9f10e53fba609e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 26 Jun 2023 17:24:41 -0400 Subject: [PATCH] +Add tv%valid_SpV_halo to debug non-Boussinesq mode Added the integer valid_SpV_halo to the thermo_var_ptrs type to indicate whether the SpV_array has been updated and its valid halo size, to facilitate error detection and debugging in non-Boussinesq mode. Tv%valid_SpV_halo is set to the halo size in calc_derived_thermo or after a halo update is done to tv%SpV_avg, and it is set to a negative value right after calls that change temperatures and salinities (such as by ALE remapping) unless there is a call to calc_derived_thermo. Tests for the validity of tv%SpV_avg are added to the routines behind thickness_to_dz, with fatal errors issued if invalid arrays would be used, but more tests could perhaps be used in any parameterization routines where tv%SpV_avg is used directly. Handling the updates to tv%SpV_avg this way helps to avoid unnecessary calls to calc_derived_thermo, which in turn has equation of state calls that can be expensive, while also providing essential verification of new code related to the non-Boussinesq code. These tests can probably be commented out or removed for efficiency once there is a full suite of regression tests for the fully non-Boussinesq mode of MOM6. In addition, a new optional debug argument was added to calc_derived_thermo which can be used to triggers checksums for the variables used to calculate tv%SpV_avg. One call to calc_derived_thermo was also added just before the initialization call to ALE_regrid that will be needed with the next commit, but does not change answers yet. All answers are bitwise identical, but there is a new element in a transparent and widely used type and a new optional argument to a public interface. --- src/ALE/MOM_ALE.F90 | 6 ++++ src/core/MOM.F90 | 24 ++++++++++---- src/core/MOM_interface_heights.F90 | 51 +++++++++++++++++++++++++----- src/core/MOM_variables.F90 | 2 ++ src/tracer/MOM_offline_main.F90 | 1 + 5 files changed, 70 insertions(+), 14 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index a341fd1835..8a82ed10dd 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -531,6 +531,7 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) ! Remap all variables from old grid h onto new grid h_new call ALE_remap_tracers(CS, G, GV, h, h_new, Reg, debug=CS%show_call_tree) + if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_inputs)") ! Reintegrate mass transports from Zstar to the offline vertical coordinate @@ -570,6 +571,8 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) h(i,j,k) = h_new(i,j,k) enddo ; enddo ; enddo + if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. + if (CS%show_call_tree) call callTree_leave("ALE_offline_inputs()") end subroutine ALE_offline_inputs @@ -672,6 +675,9 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d ! save total dzregrid for diags if needed? if (present(dzRegrid)) dzRegrid(:,:,:) = dzIntTotal(:,:,:) + + if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. + end subroutine ALE_regrid_accelerated !> This routine takes care of remapping all tracer variables between the old and the diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index af8481fd1c..d72bb20b79 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -669,8 +669,9 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS dt_therm = dt ; ntstep = 1 if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf CS%tv%p_surf => NULL() - if (associated(fluxes%p_surf)) then - if (CS%use_p_surf_in_EOS) CS%tv%p_surf => fluxes%p_surf + if (CS%use_p_surf_in_EOS .and. associated(fluxes%p_surf)) then + CS%tv%p_surf => fluxes%p_surf + if (allocated(CS%tv%SpV_avg)) call pass_var(fluxes%p_surf, G%Domain, clock=id_clock_pass) endif if (CS%UseWaves) call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass) endif @@ -1108,6 +1109,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif if (CS%interface_filter) then + if (allocated(CS%tv%SpV_avg)) call pass_var(CS%tv%SpV_avg, G%Domain, clock=id_clock_pass) + CS%tv%valid_SpV_halo = min(G%Domain%nihalo, G%Domain%njhalo) call cpu_clock_begin(id_clock_int_filter) call interface_filter(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, & CS%CDp, CS%interface_filter_CSp) @@ -1243,6 +1246,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif if (CS%interface_filter) then + if (allocated(CS%tv%SpV_avg)) call pass_var(CS%tv%SpV_avg, G%Domain, clock=id_clock_pass) + CS%tv%valid_SpV_halo = min(G%Domain%nihalo, G%Domain%njhalo) call cpu_clock_begin(id_clock_int_filter) call interface_filter(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, & CS%CDp, CS%interface_filter_CSp) @@ -1390,6 +1395,8 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) if (associated(CS%tv%T)) then call extract_diabatic_member(CS%diabatic_CSp, diabatic_halo=halo_sz) + ! The bottom boundary layer calculation may need halo values of SpV_avg, including the corners. + if (allocated(CS%tv%SpV_avg)) halo_sz = max(halo_sz, 1) if (halo_sz > 0) then call create_group_pass(pass_T_S, CS%tv%T, G%Domain, To_All, halo=halo_sz) call create_group_pass(pass_T_S, CS%tv%S, G%Domain, To_All, halo=halo_sz) @@ -1405,7 +1412,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) ! Update derived thermodynamic quantities. if (allocated(CS%tv%SpV_avg)) then - call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz) + call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz, debug=CS%debug) endif endif @@ -1557,6 +1564,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! Remap all variables from the old grid h onto the new grid h_new call ALE_remap_tracers(CS%ALE_CSp, G, GV, h, h_new, CS%tracer_Reg, showCallTree, dtdia, PCM_cell) call ALE_remap_velocities(CS%ALE_CSp, G, GV, h, h_new, u, v, CS%OBC, dzRegrid, showCallTree, dtdia) + if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. if (CS%remap_aux_vars) then if (CS%split) & @@ -1589,7 +1597,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! Update derived thermodynamic quantities. if (allocated(tv%SpV_avg)) then - call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil) + call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil, debug=CS%debug) endif if (CS%debug .and. CS%use_ALE_algorithm) then @@ -1645,7 +1653,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! Update derived thermodynamic quantities. if (allocated(tv%SpV_avg)) then - call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil) + call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil, debug=CS%debug) endif endif @@ -1818,6 +1826,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS ! are used are intended to ensure that in the case where transports don't quite conserve, ! the offline layer thicknesses do not drift too far away from the online model. call ALE_remap_tracers(CS%ALE_CSp, G, GV, CS%h, h_new, CS%tracer_Reg, debug=CS%debug) + if (allocated(CS%tv%SpV_avg)) CS%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. ! Update the tracer grid. do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 @@ -2845,6 +2854,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Allocate any derived equation of state fields. if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0) + CS%tv%valid_SpV_halo = -1 ! This array does not yet have any valid data. endif if (use_ice_shelf .and. CS%debug) then @@ -2893,6 +2903,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)") call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h ) + if (allocated(CS%tv%SpV_avg)) call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=1) call pre_ALE_adjustments(G, GV, US, CS%h, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%u, CS%v) call callTree_waypoint("Calling ALE_regrid() to remap initial conditions (initialize_MOM)") @@ -2909,6 +2920,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Remap all variables from the old grid h onto the new grid h_new call ALE_remap_tracers(CS%ALE_CSp, G, GV, CS%h, h_new, CS%tracer_Reg, CS%debug, PCM_cell=PCM_cell) call ALE_remap_velocities(CS%ALE_CSp, G, GV, CS%h, h_new, CS%u, CS%v, CS%OBC, dzRegrid, debug=CS%debug) + if (allocated(CS%tv%SpV_avg)) CS%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. ! Replace the old grid with new one. All remapping must be done at this point. !$OMP parallel do default(shared) @@ -3135,7 +3147,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Update derived thermodynamic quantities. if (allocated(CS%tv%SpV_avg)) then - call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil) + call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil, debug=CS%debug) endif if (associated(CS%visc%Kv_shear)) & diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index befeb1c2ad..dfd1048b82 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -4,13 +4,14 @@ module MOM_interface_heights ! This file is part of MOM6. See LICENSE.md for the license. use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol +use MOM_debugging, only : hchksum use MOM_error_handler, only : MOM_error, FATAL -use MOM_EOS, only : calculate_density, EOS_type, EOS_domain -use MOM_file_parser, only : log_version -use MOM_grid, only : ocean_grid_type -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : calculate_density, average_specific_vol, EOS_type, EOS_domain +use MOM_file_parser, only : log_version +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -262,7 +263,7 @@ end subroutine find_eta_2d !> Calculate derived thermodynamic quantities for re-use later. -subroutine calc_derived_thermo(tv, h, G, GV, US, halo) +subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug) 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 @@ -271,13 +272,16 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo) !! which will be set here. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. - integer, optional, intent(in) :: halo !< Width of halo within which to + integer, optional, intent(in) :: halo !< Width of halo within which to !! calculate thicknesses + logical, optional, intent(in) :: debug !< If present and true, write debugging checksums ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa] real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa] + logical :: do_debug ! If true, write checksums for debugging. integer :: i, j, k, is, ie, js, je, halos, nz + do_debug = .false. ; if (present(debug)) do_debug = debug halos = 0 ; if (present(halo)) halos = max(0,halo) is = G%isc-halos ; ie = G%iec+halos ; js = G%jsc-halos ; je = G%jec+halos ; nz = GV%ke @@ -296,6 +300,15 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo) p_t(i,j) = p_t(i,j) + dp(i,j) enddo ; enddo ; endif enddo + tv%valid_SpV_halo = halos + + if (do_debug) then + call hchksum(h, "derived_thermo h", G%HI, haloshift=halos, scale=GV%H_to_MKS) + if (associated(tv%p_surf)) call hchksum(tv%p_surf, "derived_thermo p_surf", G%HI, & + haloshift=halos, scale=US%RL2_T2_to_Pa) + call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, scale=US%C_to_degC) + call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, scale=US%S_to_ppt) + endif endif end subroutine calc_derived_thermo @@ -493,12 +506,23 @@ subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size) integer, optional, intent(in) :: halo_size !< Width of halo within which to !! calculate thicknesses ! Local variables + character(len=128) :: mesg ! A string for error messages integer :: i, j, k, is, ie, js, je, halo, nz halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then + if (tv%valid_SpV_halo < 0) then + mesg = "invalid values of SpV_avg." + else + write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') & + tv%valid_SpV_halo, halo + endif + call MOM_error(FATAL, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg)) + endif + do k=1,nz ; do j=js,je ; do i=is,ie dz(i,j,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) enddo ; enddo ; enddo @@ -529,12 +553,23 @@ subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size) integer, optional, intent(in) :: halo_size !< Width of halo within which to !! calculate thicknesses ! Local variables + character(len=128) :: mesg ! A string for error messages integer :: i, k, is, ie, halo, nz halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) is = G%isc-halo ; ie = G%iec+halo ; nz = GV%ke if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then + if (tv%valid_SpV_halo < 0) then + mesg = "invalid values of SpV_avg." + else + write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') & + tv%valid_SpV_halo, halo + endif + call MOM_error(FATAL, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg)) + endif + do k=1,nz ; do i=is,ie dz(i,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) enddo ; enddo diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index bec93376af..0c4a42e8c6 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -95,6 +95,8 @@ module MOM_variables real :: min_salinity !< The minimum value of salinity when BOUND_SALINITY=True [S ~> ppt]. real, allocatable, dimension(:,:,:) :: SpV_avg !< The layer averaged in situ specific volume [R-1 ~> m3 kg-1]. + integer :: valid_SpV_halo = -1 !< If positive, the valid halo size for SpV_avg, or if negative + !! SpV_avg is not currently set. ! These arrays are accumulated fluxes for communication with other components. real, dimension(:,:), pointer :: frazil => NULL() diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 40dced9b20..37ddf1ae0c 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -364,6 +364,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C ! Remap all variables from the old grid h_new onto the new grid h_post_remap call ALE_remap_tracers(CS%ALE_CSp, G, GV, h_new, h_post_remap, CS%tracer_Reg, & CS%debug, dt=CS%dt_offline) + if (allocated(CS%tv%SpV_avg)) CS%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 h_new(i,j,k) = h_post_remap(i,j,k)