From a61c43af83ef638bbaeb499b806d92718b58b06a Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 16 Jun 2021 14:49:33 -0400 Subject: [PATCH 1/8] Replace MOM_control_struct pointers as locals The MOM_control_struct is declared and passed as a pointer to a CS (usually allocated anonymously) and treated as if it were a pointer, even though there is currently no real advantage to doing so. After the FMS update, the deallocation of this CS was causing a segmentation fault in the PGI compilers. While the underlying cause was never determined, it is likely due to some automated deallocation of the CS contents, whose addressing became scrambled. This problem can be resolved by moving all of the CS contents to stack, so that the contents are automatically removed upon exiting whatever function it was instantiated. Subsequent calls can reference the local (or parent) stack contents. --- .../drivers/FMS_cap/ocean_model_MOM.F90 | 4 +-- .../drivers/mct_cap/mom_ocean_model_mct.F90 | 4 +-- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 4 +-- config_src/drivers/solo_driver/MOM_driver.F90 | 3 +- src/core/MOM.F90 | 34 ++++++++----------- src/framework/MOM_diag_mediator.F90 | 22 ++++++------ 6 files changed, 33 insertions(+), 38 deletions(-) diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index c3e13329f2..50ea6c943d 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -195,8 +195,8 @@ module ocean_model_mod type(unit_scale_type), pointer :: & US => NULL() !< A pointer to a structure containing dimensional !! unit scaling factors. - type(MOM_control_struct), pointer :: & - MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_control_struct) :: MOM_CSp + !< MOM control structure type(ice_shelf_CS), pointer :: & Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This diff --git a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 index 9b40a9e7b4..3bd0e1e28d 100644 --- a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 +++ b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 @@ -195,8 +195,8 @@ module MOM_ocean_model_mct !! about the vertical grid. type(unit_scale_type), pointer :: US => NULL() !< A pointer to a structure containing !! dimensional unit scaling factors. - type(MOM_control_struct), pointer :: & - MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_control_struct) :: MOM_CSp + !< MOM control structure type(ice_shelf_CS), pointer :: & Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1d5de0dd3e..7f91e15a69 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -197,8 +197,8 @@ module MOM_ocean_model_nuopc !! about the vertical grid. type(unit_scale_type), pointer :: US => NULL() !< A pointer to a structure containing !! dimensional unit scaling factors. - type(MOM_control_struct), pointer :: & - MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_control_struct) :: MOM_CSp + !< MOM control structure type(ice_shelf_CS), pointer :: & Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90 index ebb953be93..7dfce01f68 100644 --- a/config_src/drivers/solo_driver/MOM_driver.F90 +++ b/config_src/drivers/solo_driver/MOM_driver.F90 @@ -180,8 +180,7 @@ program MOM_main ! and diffusion equation are read in from files stored from ! a previous integration of the prognostic model - type(MOM_control_struct), pointer :: MOM_CSp => NULL() - !> A pointer to the tracer flow control structure. + type(MOM_control_struct) :: MOM_CSp !> MOM control structure type(tracer_flow_control_CS), pointer :: & tracer_flow_CSp => NULL() !< A pointer to the tracer flow control structure type(surface_forcing_CS), pointer :: surface_forcing_CSp => NULL() diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 05c6fe6c43..91e2ea6afd 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -434,7 +434,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS type(surface), target, intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_int_in !< time interval covered by this run segment [s]. - type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM + type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM type(Wave_parameters_CS), & optional, pointer :: Waves !< An optional pointer to a wave property CS logical, optional, intent(in) :: do_dynamics !< Present and false, do not do updates due @@ -981,7 +981,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & real, intent(in) :: bbl_time_int !< time interval over which updates to the !! bottom boundary layer properties will apply [T ~> s], !! or zero not to update the properties. - type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM + type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM type(time_type), intent(in) :: Time_local !< End time of a segment, as a time type type(wave_parameters_CS), & optional, pointer :: Waves !< Container for wave related parameters; the @@ -1432,7 +1432,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS type(surface), intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_interval !< time interval - type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM + type(MOM_control_struct), intent(inout) :: CS !< control structure from initialize_MOM ! Local pointers type(ocean_grid_type), pointer :: G => NULL() ! Pointer to a structure containing @@ -1630,7 +1630,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse type(directories), intent(out) :: dirs !< structure with directory paths - type(MOM_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure + type(MOM_control_struct), intent(inout), target :: CS !< pointer set in this routine to MOM control structure type(MOM_restart_CS), pointer :: restart_CSp !< pointer set in this routine to the !! restart control structure that will !! be used for MOM. @@ -1730,13 +1730,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & type(ocean_internal_state) :: MOM_internal_state character(len=200) :: area_varname, ice_shelf_file, inputdir, filename - if (associated(CS)) then - call MOM_error(WARNING, "initialize_MOM called with an associated "// & - "control structure.") - return - endif - allocate(CS) - CS%Time => Time id_clock_init = cpu_clock_id('Ocean Initialization', grain=CLOCK_SUBCOMPONENT) @@ -2818,7 +2811,7 @@ end subroutine initialize_MOM subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) type(time_type), intent(in) :: Time !< model time, used in this routine type(directories), intent(in) :: dirs !< structure with directory paths - type(MOM_control_struct), pointer :: CS !< pointer to MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control !! structure that will be used for MOM. ! Local variables @@ -3044,7 +3037,7 @@ end subroutine adjust_ssh_for_p_atm !! setting the appropriate fields in sfc_state. Unused fields !! are set to NULL or are unallocated. subroutine extract_surface_state(CS, sfc_state_in) - type(MOM_control_struct), pointer :: CS !< Master MOM control structure + type(MOM_control_struct), intent(inout), target :: CS !< Master MOM control structure type(surface), target, intent(inout) :: sfc_state_in !< transparent ocean surface state !! structure shared with the calling routine !! data in this structure is intent out. @@ -3471,7 +3464,7 @@ end subroutine rotate_initial_state !> Return true if all phases of step_MOM are at the same point in time. function MOM_state_is_synchronized(CS, adv_dyn) result(in_synch) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure logical, optional, intent(in) :: adv_dyn !< If present and true, only check !! whether the advection is up-to-date with !! the dynamics. @@ -3492,7 +3485,7 @@ end function MOM_state_is_synchronized !> This subroutine offers access to values or pointers to other types from within !! the MOM_control_struct, allowing the MOM_control_struct to be opaque. subroutine get_MOM_state_elements(CS, G, GV, US, C_p, C_p_scaled, use_temp) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout), target :: CS !< MOM control structure type(ocean_grid_type), optional, pointer :: G !< structure containing metrics and grid info type(verticalGrid_type), optional, pointer :: GV !< structure containing vertical grid info type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type @@ -3511,7 +3504,7 @@ end subroutine get_MOM_state_elements !> Find the global integrals of various quantities. subroutine get_ocean_stocks(CS, mass, heat, salt, on_PE_only) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure real, optional, intent(out) :: heat !< The globally integrated integrated ocean heat [J]. real, optional, intent(out) :: salt !< The globally integrated integrated ocean salt [kg]. real, optional, intent(out) :: mass !< The globally integrated integrated ocean mass [kg]. @@ -3528,7 +3521,7 @@ end subroutine get_ocean_stocks !> End of ocean model, including memory deallocation subroutine MOM_end(CS) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure call MOM_sum_output_end(CS%sum_output_CSp) @@ -3604,7 +3597,6 @@ subroutine MOM_end(CS) if (associated(CS%update_OBC_CSp)) call OBC_register_end(CS%update_OBC_CSp) call verticalGridEnd(CS%GV) - call unit_scaling_end(CS%US) call MOM_grid_end(CS%G) if (CS%debug .or. CS%G%symmetric) & @@ -3613,9 +3605,11 @@ subroutine MOM_end(CS) if (CS%rotate_index) & call deallocate_MOM_domain(CS%G%Domain) - call deallocate_MOM_domain(CS%G_in%domain) + ! The MPP domains may be needed by an external coupler, so use `cursory`. + ! TODO: This may create a domain memory leak, and needs investigation. + call deallocate_MOM_domain(CS%G_in%domain, cursory=.true.) - deallocate(CS) + call unit_scaling_end(CS%US) end subroutine MOM_end !> \namespace mom diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index e068d26f5d..4994646086 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3473,16 +3473,18 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) call axes_grp_end(diag_cs%remap_axesCvi(i)) enddo - deallocate(diag_cs%remap_axesZL) - deallocate(diag_cs%remap_axesZi) - deallocate(diag_cs%remap_axesTL) - deallocate(diag_cs%remap_axesTi) - deallocate(diag_cs%remap_axesBL) - deallocate(diag_cs%remap_axesBi) - deallocate(diag_cs%remap_axesCuL) - deallocate(diag_cs%remap_axesCui) - deallocate(diag_cs%remap_axesCvL) - deallocate(diag_cs%remap_axesCvi) + if (diag_cs%num_diag_coords > 0) then + deallocate(diag_cs%remap_axesZL) + deallocate(diag_cs%remap_axesZi) + deallocate(diag_cs%remap_axesTL) + deallocate(diag_cs%remap_axesTi) + deallocate(diag_cs%remap_axesBL) + deallocate(diag_cs%remap_axesBi) + deallocate(diag_cs%remap_axesCuL) + deallocate(diag_cs%remap_axesCui) + deallocate(diag_cs%remap_axesCvL) + deallocate(diag_cs%remap_axesCvi) + endif do dl=2,MAX_DSAMP_LEV if (allocated(diag_cs%dsamp(dl)%remap_axesTL)) & From d7dece9717c629ca92c7dcff2a2c2bb1723920cf Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 1 Jul 2021 20:37:15 -0400 Subject: [PATCH 2/8] +Add wind-stress acceleration diagnostics Added new diagnostics of the acceleration driven by the wind stress accelerations as redistributed by a timestep of the vertical viscosity and not lost to bottom drag within a timestep. This is also in the diagnostics of the accelerations due to the vertical viscosity, but the redistribution can be found from the difference of the two. Also added a diagnostic of the contribution of the wind stresses to kinetic energy, and applied an underflow limiter on both the new acceleration diagnostic and the existing viscous acceleration diagnostic. All solutions are bitwise identical, but there are new diagnostics. --- src/core/MOM_variables.F90 | 4 + src/diagnostics/MOM_diagnostics.F90 | 54 ++++++++--- .../vertical/MOM_vert_friction.F90 | 95 ++++++++++++++----- 3 files changed, 119 insertions(+), 34 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f966ab2ad2..3f98a97052 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -169,6 +169,10 @@ module MOM_variables PFv => NULL(), & !< Meridional acceleration due to pressure forces [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_str => NULL(), & !< Zonal acceleration due to the surface stress (included + !! in du_dt_visc) [L T-2 ~> m s-2] + dv_dt_str => NULL(), & !< Meridional acceleration due to the surface stress (included + !! in dv_dt_visc) [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 e6b01af33d..44b05cc081 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -107,6 +107,7 @@ module MOM_diagnostics !! of this spurious Coriolis source. 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_stress => NULL(), & !< KE source from surface stress (included in KE_visc) [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] @@ -121,8 +122,8 @@ module MOM_diagnostics integer :: id_col_ht = -1, id_dh_dt = -1 integer :: id_KE = -1, id_dKEdt = -1 integer :: id_PE_to_KE = -1, id_KE_BT = -1 - integer :: id_KE_Coradv = -1 - integer :: id_KE_adv = -1, id_KE_visc = -1 + integer :: id_KE_Coradv = -1, id_KE_adv = -1 + integer :: id_KE_visc = -1, id_KE_stress = -1 integer :: id_KE_horvisc = -1, id_KE_dia = -1 integer :: id_uh_Rlay = -1, id_vh_Rlay = -1 integer :: id_uhGM_Rlay = -1, id_vhGM_Rlay = -1 @@ -1060,7 +1061,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_h(i,j) = CS%KE(i,j,k) * CS%dh_dt(i,j,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%dKE_dt(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)) @@ -1078,7 +1079,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%PE_to_KE(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)) @@ -1096,7 +1097,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_BT(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)) @@ -1118,7 +1119,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS * (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) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_CorAdv(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)) @@ -1146,7 +1147,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS * (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) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_adv(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)) @@ -1164,7 +1165,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + 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)) @@ -1173,6 +1174,24 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (CS%id_KE_visc > 0) call post_data(CS%id_KE_visc, CS%KE_visc, CS%diag) endif + if (associated(CS%KE_stress)) 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_str(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_str(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_stress(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_stress > 0) call post_data(CS%id_KE_stress, CS%KE_stress, CS%diag) + endif + if (associated(CS%KE_horvisc)) then do k=1,nz do j=js,je ; do I=Isq,Ieq @@ -1182,7 +1201,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_horvisc(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)) @@ -1203,7 +1222,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_h(i,j) = CS%KE(i,j,k) * (CDp%diapyc_vel(i,j,k) - CDp%diapyc_vel(i,j,k+1)) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_dia(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)) @@ -1894,6 +1913,11 @@ 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_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz) + CS%id_KE_stress = register_diag_field('ocean_model', 'KE_stress', diag%axesTL, Time, & + 'Kinetic Energy Source from Surface Stresses or Body Wind Stress', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_KE_stress>0) call safe_alloc_ptr(CS%KE_stress,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) @@ -2294,7 +2318,7 @@ 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_adv) .or. associated(CS%KE_visc) .or. associated(CS%KE_stress) .or. & associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) & call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) @@ -2323,6 +2347,11 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) endif + if (associated(CS%KE_stress)) then + call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%dv_dt_str,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) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) @@ -2353,6 +2382,7 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp) if (associated(CS%KE_Coradv)) deallocate(CS%KE_Coradv) if (associated(CS%KE_adv)) deallocate(CS%KE_adv) if (associated(CS%KE_visc)) deallocate(CS%KE_visc) + if (associated(CS%KE_stress)) deallocate(CS%KE_stress) if (associated(CS%KE_horvisc)) deallocate(CS%KE_horvisc) if (associated(CS%KE_dia)) deallocate(CS%KE_dia) if (associated(CS%dv_dt)) deallocate(CS%dv_dt) @@ -2368,6 +2398,8 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp) if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu) 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_str)) deallocate(ADp%du_dt_str) + if (associated(ADp%dv_dt_str)) deallocate(ADp%dv_dt_str) 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) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 5b85c5f5f6..1d46f9aee3 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -121,6 +121,7 @@ module MOM_vert_friction !>@{ Diagnostic identifiers integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 + integer :: id_du_dt_str = -1, id_dv_dt_str = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 @@ -207,6 +208,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real :: stress ! The surface stress times the time step, divided ! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1]. + real :: accel_underflow ! An acceleration magnitude that is so small that values that are less + ! than this are diagnosed as 0 [L T-2 ~> m s-2]. real :: zDS, hfr, h_a ! Temporary variables used with direct_stress. real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1]. @@ -236,6 +239,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & h_neglect = GV%H_subroundoff Idt = 1.0 / dt + accel_underflow = CS%vel_underflow * Idt + !Check if Stokes mixing allowed if requested (present and associated) DoStokesMixing=.false. if (CS%StokesMixing) then @@ -265,9 +270,13 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ADp%du_dt_visc(I,j,k) = u(I,j,k) enddo ; enddo ; endif -! One option is to have the wind stress applied as a body force -! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, -! the wind stress is applied as a stress boundary condition. + if (associated(ADp%du_dt_str)) then ; do k=1,nz ; do I=Isq,Ieq + ADp%du_dt_str(I,j,k) = 0.0 + enddo ; enddo ; endif + + ! One option is to have the wind stress applied as a body force + ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, + ! the wind stress is applied as a stress boundary condition. if (CS%direct_stress) then do I=Isq,Ieq ; if (do_i(I)) then surface_stress(I) = 0.0 @@ -277,6 +286,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & h_a = 0.5 * (h(I,j,k) + h(I+1,j,k)) + h_neglect hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a u(I,j,k) = u(I,j,k) + I_Hmix * hfr * stress + if (associated(ADp%du_dt_str)) ADp%du_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt zDS = zDS + h_a ; if (zDS >= Hmix) exit enddo endif ; enddo ! end of i loop @@ -316,6 +326,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_u(I,j,2)) d1(I) = b_denom_1 * b1(I) u(I,j,1) = b1(I) * (CS%h_u(I,j,1) * u(I,j,1) + surface_stress(I)) + if (associated(ADp%du_dt_str)) & + ADp%du_dt_str(I,j,1) = b1(I) * (CS%h_u(I,j,1) * ADp%du_dt_str(I,j,1) + surface_stress(I)*Idt) endif ; enddo do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then c1(I,k) = dt_Z_to_H * CS%a_u(I,j,K) * b1(I) @@ -324,6 +336,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & d1(I) = b_denom_1 * b1(I) u(I,j,k) = (CS%h_u(I,j,k) * u(I,j,k) + & dt_Z_to_H * CS%a_u(I,j,K) * u(I,j,k-1)) * b1(I) + if (associated(ADp%du_dt_str)) & + ADp%du_dt_str(I,j,k) = (CS%h_u(I,j,k) * ADp%du_dt_str(I,j,k) + & + dt_Z_to_H * CS%a_u(I,j,K) * ADp%du_dt_str(I,j,k-1)) * b1(I) endif ; enddo ; enddo ! back substitute to solve for the new velocities @@ -332,8 +347,17 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & u(I,j,k) = u(I,j,k) + c1(I,k+1) * u(I,j,k+1) endif ; enddo ; enddo ! i and k loops + if (associated(ADp%du_dt_str)) then + do i=is,ie ; if (abs(ADp%du_dt_str(I,j,nz)) < accel_underflow) ADp%du_dt_str(I,j,nz) = 0.0 ; enddo + do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then + ADp%du_dt_str(I,j,k) = ADp%du_dt_str(I,j,k) + c1(I,k+1) * ADp%du_dt_str(I,j,k+1) + if (abs(ADp%du_dt_str(I,j,k)) < accel_underflow) ADp%du_dt_str(I,j,k) = 0.0 + endif ; enddo ; enddo + endif + 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 + if (abs(ADp%du_dt_visc(I,j,k)) < accel_underflow) ADp%du_dt_visc(I,j,k) = 0.0 enddo ; enddo ; endif if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq @@ -373,9 +397,13 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ADp%dv_dt_visc(i,J,k) = v(i,J,k) enddo ; enddo ; endif -! One option is to have the wind stress applied as a body force -! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, -! the wind stress is applied as a stress boundary condition. + if (associated(ADp%dv_dt_str)) then ; do k=1,nz ; do i=is,ie + ADp%dv_dt_str(i,J,k) = 0.0 + enddo ; enddo ; endif + + ! One option is to have the wind stress applied as a body force + ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, + ! the wind stress is applied as a stress boundary condition. if (CS%direct_stress) then do i=is,ie ; if (do_i(i)) then surface_stress(i) = 0.0 @@ -385,6 +413,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & h_a = 0.5 * (h(i,J,k) + h(i,J+1,k)) + h_neglect hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a v(i,J,k) = v(i,J,k) + I_Hmix * hfr * stress + if (associated(ADp%dv_dt_str)) ADp%dv_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt zDS = zDS + h_a ; if (zDS >= Hmix) exit enddo endif ; enddo ! end of i loop @@ -401,6 +430,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_v(i,J,2)) d1(i) = b_denom_1 * b1(i) v(i,J,1) = b1(i) * (CS%h_v(i,J,1) * v(i,J,1) + surface_stress(i)) + if (associated(ADp%dv_dt_str)) & + ADp%dv_dt_str(i,J,1) = b1(i) * (CS%h_v(i,J,1) * ADp%dv_dt_str(i,J,1) + surface_stress(i)*Idt) endif ; enddo do k=2,nz ; do i=is,ie ; if (do_i(i)) then c1(i,k) = dt_Z_to_H * CS%a_v(i,J,K) * b1(i) @@ -408,13 +439,25 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_v(i,J,K+1)) d1(i) = b_denom_1 * b1(i) v(i,J,k) = (CS%h_v(i,J,k) * v(i,J,k) + dt_Z_to_H * CS%a_v(i,J,K) * v(i,J,k-1)) * b1(i) + if (associated(ADp%dv_dt_str)) & + ADp%dv_dt_str(i,J,k) = (CS%h_v(i,J,k) * ADp%dv_dt_str(i,J,k) + & + dt_Z_to_H * CS%a_v(i,J,K) * ADp%dv_dt_str(i,J,k-1)) * b1(i) endif ; enddo ; enddo do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then v(i,J,k) = v(i,J,k) + c1(i,k+1) * v(i,J,k+1) endif ; enddo ; enddo ! i and k loops + if (associated(ADp%dv_dt_str)) then + do i=is,ie ; if (abs(ADp%dv_dt_str(i,J,nz)) < accel_underflow) ADp%dv_dt_str(i,J,nz) = 0.0 ; enddo + do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then + ADp%dv_dt_str(i,J,k) = ADp%dv_dt_str(i,J,k) + c1(i,k+1) * ADp%dv_dt_str(i,J,k+1) + if (abs(ADp%dv_dt_str(i,J,k)) < accel_underflow) ADp%dv_dt_str(i,J,k) = 0.0 + endif ; enddo ; enddo + endif + 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 + if (abs(ADp%dv_dt_visc(i,J,k)) < accel_underflow) ADp%dv_dt_visc(i,J,k) = 0.0 enddo ; enddo ; endif if (associated(visc%tauy_shelf)) then ; do i=is,ie @@ -458,7 +501,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & enddo endif -! Offer diagnostic fields for averaging. + ! Offer diagnostic fields for averaging. if (CS%id_du_dt_visc > 0) & call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag) if (CS%id_dv_dt_visc > 0) & @@ -467,6 +510,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & call post_data(CS%id_taux_bot, taux_bot, CS%diag) if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) & call post_data(CS%id_tauy_bot, tauy_bot, CS%diag) + if (CS%id_du_dt_str > 0) & + call post_data(CS%id_du_dt_str, ADp%du_dt_str, CS%diag) + if (CS%id_dv_dt_str > 0) & + call post_data(CS%id_dv_dt_str, ADp%dv_dt_str, CS%diag) ! Diagnostics for terms multiplied by fractional thicknesses @@ -524,10 +571,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & end subroutine vertvisc -!> Calculate the fraction of momentum originally in a layer that remains -!! after a time-step of viscosity, and the fraction of a time-step's -!! worth of barotropic acceleration that a layer experiences after -!! viscosity is applied. +!> Calculate the fraction of momentum originally in a layer that remains in the water column +!! after a time-step of viscosity, equivalently the fraction of a time-step's worth of +!! barotropic acceleration that a layer experiences after viscosity is applied. subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -566,10 +612,8 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo - ! Find the zonal viscous using a modification of a standard tridagonal solver. -!$OMP parallel do default(none) shared(G,Isq,Ieq,CS,nz,visc,dt_Z_to_H,visc_rem_u) & -!$OMP firstprivate(Ray) & -!$OMP private(do_i,b_denom_1,b1,d1,c1) + ! Find the zonal viscous remnant using a modification of a standard tridagonal solver. + !$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1) do j=G%jsc,G%jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo @@ -597,10 +641,8 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) enddo ! end u-component j loop - ! Now find the meridional viscous using a modification. -!$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,G,CS,visc,dt_Z_to_H,visc_rem_v,nz) & -!$OMP firstprivate(Ray) & -!$OMP private(do_i,b_denom_1,b1,d1,c1) + ! Now find the meridional viscous remnant using the robust tridiagonal solver. + !$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1) do J=Jsq,Jeq do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo @@ -1813,13 +1855,20 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', & thickness_units, conversion=GV%H_to_MKS) - CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, & - Time, 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, & + 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_du_dt_visc > 0) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) - CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, & - Time, 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, Time, & + 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_dv_dt_visc > 0) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + CS%id_du_dt_str = register_diag_field('ocean_model', 'du_dt_str', diag%axesCuL, Time, & + 'Zonal Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_du_dt_str > 0) call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz) + CS%id_dv_dt_str = register_diag_field('ocean_model', 'dv_dt_str', diag%axesCvL, Time, & + 'Meridional Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_dv_dt_str > 0) call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz) + CS%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, & Time, 'Zonal Bottom Stress from Ocean to Earth', & 'Pa', conversion=US%RZ_to_kg_m2*US%L_T2_to_m_s2) From d9bdbc3431abdc2fe36d75071fd9ef0ec43c2a06 Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Fri, 9 Jul 2021 12:45:06 -0600 Subject: [PATCH 3/8] fix dim_names assignment in MOM_io:read_var_sizes fixes #1439 --- src/framework/MOM_io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 3ad2c92f41..f8cfb09382 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -714,7 +714,7 @@ subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, d if (status /= NF90_NOERR) call MOM_error(WARNING, trim(hdr) // trim(NF90_STRERROR(status)) //& " Getting dimension length for "//trim(varname)//" in "//trim(filename)) if (present(dim_names)) then - if (n <= size(dim_names)) dim_names = trim(dimname) + if (n <= size(dim_names)) dim_names(n) = trim(dimname) endif enddo deallocate(dimids) From 7c27bfba2bf27c6d78bd13792906c4438731c186 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 9 Jul 2021 18:09:11 -0400 Subject: [PATCH 4/8] +Add query_wave_properties & fix NUOPC wave queries Added a new routine, query_wave_properties, that can be call to get information about the wave properties from the waves control structure, and added a get_param call for SURFBAND_WAVENUMBERS to MOM_wave_interface_init when it is using the options that are typical with NUOPC coupler. These changes should allow the NUOPC coupler to compile and work again, while still keeping the same level of opacity in the wave_parameters_CS. All answers should be bitwise identical, although the order of some entries in the MOM_parameter_doc files may change when coupled with waves is WAVE_METHOD=SURFACE_BANDS and SURFBAND_SOURCE=COUPLER. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 8 ++++-- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 6 ---- src/user/MOM_wave_interface.F90 | 28 +++++++++++++++++++ 3 files changed, 33 insertions(+), 9 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2d79674606..6561f63fd1 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -28,8 +28,9 @@ module MOM_cap_mod use MOM_get_input, only: get_MOM_input, directories use MOM_domains, only: pass_var use MOM_error_handler, only: MOM_error, FATAL, is_root_pe -use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type, get_global_grid_size +use MOM_wave_interface, only: query_wave_properties +use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type use MOM_ocean_model_nuopc, only: ocean_model_init_sfc use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end @@ -696,7 +697,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%frunoff = 0.0 if (ocean_state%use_waves) then - Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands + call query_wave_properties(ocean_state%Waves, NumBands=Ice_ocean_boundary%num_stk_bands) allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & @@ -704,7 +705,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands)) Ice_ocean_boundary%ustk0 = 0.0 Ice_ocean_boundary%vstk0 = 0.0 - Ice_ocean_boundary%stk_wavenumbers = ocean_state%Waves%WaveNum_Cen + call query_wave_properties(ocean_state%Waves, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, & + US=ocean_state%US) Ice_ocean_boundary%ustkb = 0.0 Ice_ocean_boundary%vstkb = 0.0 endif diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1d5de0dd3e..06fa905a23 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -391,12 +391,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! MOM_wave_interface_init is called regardless of the value of USE_WAVES because ! it also initializes statistical waves. call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) - if (OS%use_waves) then - ! I do not know why this is being set here. It seems out of place. -RWH - call get_param(param_file,mdl,"SURFBAND_WAVENUMBERS", OS%Waves%WaveNum_Cen, & - "Central wavenumbers for surface Stokes drift bands.", & - units='rad/m', default=0.12566, scale=OS%US%Z_to_m) - endif if (associated(OS%grid%Domain%maskmap)) then call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index e1e4ab0f77..d51f6ae46a 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -24,6 +24,7 @@ module MOM_wave_interface #include public MOM_wave_interface_init ! Public interface to fully initialize the wave routines. +public query_wave_properties ! Public interface to obtain information from the waves control structure. public Update_Surface_Waves ! Public interface to update wave information at the ! coupler/driver level. public Update_Stokes_Drift ! Public interface to update the Stokes drift profiles @@ -322,6 +323,9 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) "STK_BAND_COUPLER is the number of Stokes drift bands in the coupler. "// & "This has to be consistent with the number of Stokes drift bands in WW3, "//& "or the model will fail.",units='', default=1) + call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & + "Central wavenumbers for surface Stokes drift bands.", & + units='rad/m', default=0.12566, scale=US%Z_to_m) allocate( CS%WaveNum_Cen(CS%NumBands) ) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) @@ -428,6 +432,30 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) end subroutine MOM_wave_interface_init +!> This interface provides the caller with information from the waves control structure. +subroutine query_wave_properties(CS, NumBands, WaveNumbers, US) + type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure + integer, optional, intent(out) :: NumBands !< If present, this returns the number of + !!< wavenumber partitions in the wave discretization + real, dimension(:), optional, intent(out) :: Wavenumbers !< If present this returns the characteristic + !! wavenumbers of the wave discretization [m-1 or Z-1 ~> m-1] + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type that is used to undo + !! the dimensional scaling of the output variables, if present + integer :: n + + if (present(NumBands)) NumBands = CS%NumBands + if (present(Wavenumbers)) then + if (size(Wavenumbers) < CS%NumBands) call MOM_error(FATAL, "query_wave_properties called "//& + "with a Wavenumbers array that is smaller than the number of bands.") + if (present(US)) then + do n=1,CS%NumBands ; Wavenumbers(n) = US%m_to_Z * CS%WaveNum_Cen(n) ; enddo + else + do n=1,CS%NumBands ; Wavenumbers(n) = CS%WaveNum_Cen(n) ; enddo + endif + endif + +end subroutine query_wave_properties + !> Subroutine that handles updating of surface wave/Stokes drift related properties subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS, forces) type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure From f6c1fc73a700d26b4fc84014ad69b93f956524ec Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 12 Jul 2021 18:47:14 +0000 Subject: [PATCH 5/8] Allocate Waves%WaveNum_Cen before reading into it Moved the recently added call to read SURFBAND_WAVENUMBERS into an array when SURFBAND_SOURCE=COUPLER down by a line to follow the allocate call for that array. This change should avert a memory access problem that would otherwise arise when exercising this newly added code. All answers are bitwise identical in any case that ran with the previous version, and they should reproduce answers from before this PR as a whole, although obviously this code is not as well tested as would be ideal, based on the fact that this bug was in the last commit, which also passed testing. --- src/user/MOM_wave_interface.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index d51f6ae46a..b612ef4270 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -323,10 +323,10 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) "STK_BAND_COUPLER is the number of Stokes drift bands in the coupler. "// & "This has to be consistent with the number of Stokes drift bands in WW3, "//& "or the model will fail.",units='', default=1) + allocate( CS%WaveNum_Cen(CS%NumBands) ) call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & "Central wavenumbers for surface Stokes drift bands.", & units='rad/m', default=0.12566, scale=US%Z_to_m) - allocate( CS%WaveNum_Cen(CS%NumBands) ) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) CS%WaveNum_Cen(:) = 0.0 From 01b51e825a8312b0730394f565641f4e38189953 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Jul 2021 10:23:37 -0400 Subject: [PATCH 6/8] +Add query_ocean_state Added the new routine query_ocean_state to mom_ocean_model_nuopc to allow for the wave properties to be obtained by the mom_cap without having to rely on the elements of an otherwise opaque type being public. I believe that all elements of the ocean_state_type could now be declared as private, and that the version of MOM6 with the NUOPC coupler should now work as intended. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 16 +++++----- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 29 +++++++++++++++++-- 2 files changed, 36 insertions(+), 9 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 6561f63fd1..394cf05285 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -29,12 +29,11 @@ module MOM_cap_mod use MOM_domains, only: pass_var use MOM_error_handler, only: MOM_error, FATAL, is_root_pe use MOM_grid, only: ocean_grid_type, get_global_grid_size -use MOM_wave_interface, only: query_wave_properties use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type use MOM_ocean_model_nuopc, only: ocean_model_init_sfc use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end -use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh +use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh, query_ocean_state use MOM_cap_time, only: AlarmInit use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor use MOM_cap_methods, only: med2mod_areacor, state_diagnose @@ -422,6 +421,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) character(len=64) :: logmsg logical :: isPresent, isPresentDiro, isPresentLogfile, isSet logical :: existflag + logical :: use_waves ! If true, the wave modules are active. integer :: userRc integer :: localPet integer :: localPeCount @@ -696,8 +696,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%lrunoff = 0.0 Ice_ocean_boundary%frunoff = 0.0 - if (ocean_state%use_waves) then - call query_wave_properties(ocean_state%Waves, NumBands=Ice_ocean_boundary%num_stk_bands) + call query_ocean_state(ocean_state, use_waves=use_waves) + if (use_waves) then + call query_ocean_state(ocean_state, NumWaveBands=Ice_ocean_boundary%num_stk_bands) allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & @@ -705,11 +706,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands)) Ice_ocean_boundary%ustk0 = 0.0 Ice_ocean_boundary%vstk0 = 0.0 - call query_wave_properties(ocean_state%Waves, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, & - US=ocean_state%US) + call query_ocean_state(ocean_state, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, unscale=.true.) Ice_ocean_boundary%ustkb = 0.0 Ice_ocean_boundary%vstkb = 0.0 endif + ! Consider adding this: + ! if (.not.use_waves) Ice_ocean_boundary%num_stk_bands = 0 ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state call ESMF_GridCompSetInternalState(gcomp, ocean_internalstate, rc) @@ -754,7 +756,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !These are not currently used and changing requires a nuopc dictionary change !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") - if (ocean_state%use_waves) then + if (use_waves) then if (Ice_ocean_boundary%num_stk_bands > 3) then call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this") endif diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 06fa905a23..03f2b1ed9d 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -58,7 +58,7 @@ module MOM_ocean_model_nuopc use mpp_mod, only : mpp_chksum use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_wave_interface, only : wave_parameters_CS, MOM_wave_interface_init -use MOM_wave_interface, only : Update_Surface_Waves +use MOM_wave_interface, only : Update_Surface_Waves, query_wave_properties use MOM_surface_forcing_nuopc, only : surface_forcing_init, convert_IOB_to_fluxes use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS @@ -80,7 +80,7 @@ module MOM_ocean_model_nuopc public ocean_model_restart public ice_ocn_bnd_type_chksum public ocean_public_type_chksum -public get_ocean_grid +public get_ocean_grid, query_ocean_state public get_eps_omesh !> This type is used for communication with other components via the FMS coupler. @@ -1000,6 +1000,31 @@ subroutine ocean_model_flux_init(OS, verbosity) end subroutine ocean_model_flux_init +!> This interface allows certain properties that are stored in the ocean_state_type to be +!! obtained. +subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale) + type(ocean_state_type), intent(in) :: OS !< The structure with the complete ocean state + logical, optional, intent(out) :: use_waves !< Indicates whether surface waves are in use + integer, optional, intent(out) :: NumWaveBands !< If present, this gives the number of + !! wavenumber partitions in the wave discretization + real, dimension(:), optional, intent(out) :: Wavenumbers !< If present, this gives the characteristic + !! wavenumbers of the wave discretization [m-1 or Z-1 ~> m-1] + logical, optional, intent(in) :: unscale !< If present and true, undo any dimensional + !! rescaling and return dimensional values in MKS units + + logical :: undo_scaling + undo_scaling = .false. ; if (present(unscale)) undo_scaling = unscale + + if (present(use_waves)) use_waves = OS%use_waves + if (present(NumWaveBands)) call query_wave_properties(OS%Waves, NumBands=NumWaveBands) + if (present(Wavenumbers) .and. undo_scaling) then + call query_wave_properties(OS%Waves, WaveNumbers=WaveNumbers, US=OS%US) + elseif (present(Wavenumbers)) then + call query_wave_properties(OS%Waves, WaveNumbers=WaveNumbers) + endif + +end subroutine query_ocean_state + !> Ocean_stock_pe - returns the integrated stocks of heat, water, etc. for conservation checks. !! Because of the way FMS is coded, only the root PE has the integrated amount, !! while all other PEs get 0. From 2ff2419eb086a6eac8d02322ea2dabd734de556f Mon Sep 17 00:00:00 2001 From: jiandewang Date: Sun, 18 Jul 2021 13:03:38 -0400 Subject: [PATCH 7/8] initialize CS%WaveNum_Cen before read in this parameter --- src/user/MOM_wave_interface.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index b612ef4270..38aa6b13a5 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -324,15 +324,15 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) "This has to be consistent with the number of Stokes drift bands in WW3, "//& "or the model will fail.",units='', default=1) allocate( CS%WaveNum_Cen(CS%NumBands) ) - call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & - "Central wavenumbers for surface Stokes drift bands.", & - units='rad/m', default=0.12566, scale=US%Z_to_m) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) CS%WaveNum_Cen(:) = 0.0 CS%STKx0(:,:,:) = 0.0 CS%STKy0(:,:,:) = 0.0 CS%PartitionMode = 0 + call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & + "Central wavenumbers for surface Stokes drift bands.", & + units='rad/m', default=0.12566, scale=US%Z_to_m) case (INPUT_STRING)! A method to input the Stokes band (globally uniform) CS%DataSource = INPUT call get_param(param_file,mdl,"SURFBAND_NB",CS%NumBands, & From 91dabbac277eb5b5a399d96bad9d482bfe67420e Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 26 Jul 2021 13:56:42 -0800 Subject: [PATCH 8/8] (*)Solution to UV sponge restarts (#1434) * Adding a pass_var to surface h This is required for the u,v sponges to be invariant of tiling. I don't know why, but the problem only showed up for me in a narrow channel in the Bering domain. * Matt's suggestion for not sponging at the grid edge. * Finish the masking out of edge in uv sponge. - Without this change, the edges don't reproduce on restart due to the h values outside being nonsense. The non-merge commits in this squashed PR that were not already in dev/gfdl are: https://github.com/NOAA-GFDL/MOM6/pull/1434/commits/2bfa4bce5aaecdaacc8606485b0c2fa164efb8f1 "Adding a pass_var to surface h" https://github.com/NOAA-GFDL/MOM6/pull/1434/commits/95a770ac76882b5c6da94e5ad3c3e5bdd986c7a7 "Matt's suggestion for not sponging at the grid edge." and https://github.com/NOAA-GFDL/MOM6/pull/1434/commits/53dfdc71dbb317317a393509c1c4250b8fb3a5a8 "Finish the masking out of edge in uv sponge." --- .../vertical/MOM_ALE_sponge.F90 | 41 ++++++++++--------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 31d2ab5a76..e122452368 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -16,7 +16,7 @@ module MOM_ALE_sponge use MOM_coms, only : sum_across_PEs use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field use MOM_diag_mediator, only : diag_ctrl -use MOM_domains, only : pass_var +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_error, FATAL, NOTE, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -709,7 +709,6 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, integer :: nz_data !< the number of vertical levels in this input field character(len=256) :: mesg ! String for error messages ! Local variables for ALE remapping - real, dimension(:), allocatable :: tmpT1d real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m]. type(remapping_CS) :: remapCS ! Remapping parameters and work arrays @@ -1019,30 +1018,31 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& answers_2018=CS%hor_regrid_answers_2018) - call pass_var(sp_val,G%Domain) + call pass_var(sp_val, G%Domain) + call pass_var(mask_z, G%Domain) do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB sp_val_u(I,j,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i+1,j,1:nz_data)) - mask_u(I,j,1:nz_data) = max(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data)) + mask_u(I,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data)) enddo ; enddo allocate( hsrc(nz_data) ) - allocate( tmpT1d(nz_data) ) do c=1,CS%num_col_u ! c is an index for the next 3 lines but a multiplier for the rest of the loop ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i_u(c) ; j = CS%col_j_u(c) - CS%Ref_val_u%p(1:nz_data,c) = sp_val_u(i,j,1:nz_data) + if (mask_u(i,j,1) == 1.0) then + CS%Ref_val_u%p(1:nz_data,c) = sp_val_u(i,j,1:nz_data) + else + CS%Ref_val_u%p(1:nz_data,c) = 0.0 + endif ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 do k=1,nz_data if (mask_u(i,j,k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) - tmpT1d(k) = sp_val_u(i,j,k) elseif (k>1) then zBottomOfCell = -G%bathyT(i,j) - tmpT1d(k) = tmpT1d(k-1) else ! This next block should only ever be reached over land - tmpT1d(k) = -99.9 endif hsrc(k) = zTopOfCell - zBottomOfCell if (hsrc(k)>0.) nPoints = nPoints + 1 @@ -1052,7 +1052,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) CS%Ref_val_u%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) enddo - deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc, tmpT1d) + deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc) nz_data = CS%Ref_val_v%nz_data allocate(sp_val( G%isd:G%ied,G%jsd:G%jed,1:nz_data)) allocate(sp_val_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data)) @@ -1066,30 +1066,31 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& answers_2018=CS%hor_regrid_answers_2018) - call pass_var(sp_val,G%Domain) + call pass_var(sp_val, G%Domain) + call pass_var(mask_z, G%Domain) do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec sp_val_v(i,J,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i,j+1,1:nz_data)) - mask_v(i,J,1:nz_data) = max(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data)) + mask_v(i,J,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data)) enddo ; enddo !call pass_var(mask_z,G%Domain) allocate( hsrc(nz_data) ) - allocate( tmpT1d(nz_data) ) do c=1,CS%num_col_v ! c is an index for the next 3 lines but a multiplier for the rest of the loop ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i_v(c) ; j = CS%col_j_v(c) - CS%Ref_val_v%p(1:nz_data,c) = sp_val_v(i,j,1:nz_data) + if (mask_v(i,j,1) == 1.0) then + CS%Ref_val_v%p(1:nz_data,c) = sp_val_v(i,j,1:nz_data) + else + CS%Ref_val_v%p(1:nz_data,c) = 0.0 + endif ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 do k=1,nz_data if (mask_v(i,j,k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) - tmpT1d(k) = sp_val_v(i,j,k) elseif (k>1) then zBottomOfCell = -G%bathyT(i,j) - tmpT1d(k) = tmpT1d(k-1) else ! This next block should only ever be reached over land - tmpT1d(k) = -99.9 endif hsrc(k) = zTopOfCell - zBottomOfCell if (hsrc(k)>0.) nPoints = nPoints + 1 @@ -1099,7 +1100,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) CS%Ref_val_v%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) enddo - deallocate(sp_val, sp_val_v, mask_v, mask_z, hsrc, tmpT1d) + deallocate(sp_val, sp_val_v, mask_v, mask_z, hsrc) endif call pass_var(h,G%Domain)