From a22f69190527804ee03c2e67d08bbc546d2f6c39 Mon Sep 17 00:00:00 2001 From: WenhaoChen89 <96131003+WenhaoChen89@users.noreply.github.com> Date: Fri, 19 Aug 2022 14:50:16 -0400 Subject: [PATCH 01/15] Fix ALE sponge diagnostics (#188) * Change dumbbell initialization * Change in Dumbbell Layer Mode * Fix sponge diagnostics * Fix ALE Sponge Diagnostics * A minor style change removing spaces around = in optional. function arguments. * Fix ALE sponge diagnostics * character declaration fix --- .../MOM_state_initialization.F90 | 12 ++-- .../vertical/MOM_ALE_sponge.F90 | 59 +++++++++++++++---- src/tracer/RGC_tracer.F90 | 2 +- src/user/DOME2d_initialization.F90 | 10 ++-- src/user/ISOMIP_initialization.F90 | 11 ++-- src/user/RGC_initialization.F90 | 6 +- src/user/dense_water_initialization.F90 | 6 +- src/user/dumbbell_initialization.F90 | 7 ++- 8 files changed, 79 insertions(+), 34 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index b8e74e3c45..58a6f9ed8d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2135,9 +2135,11 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t if (use_temperature) then allocate(tmp_tr(isd:ied,jsd:jed,nz_data)) call MOM_read_data(filename, potemp_var, tmp_tr(:,:,:), G%Domain, scale=US%degC_to_C) - call set_up_ALE_sponge_field(tmp_tr, G, GV, tv%T, ALE_CSp) + call set_up_ALE_sponge_field(tmp_tr, G, GV, tv%T, ALE_CSp, 'temp', & + sp_long_name='temperature', sp_unit='degC s-1') call MOM_read_data(filename, salin_var, tmp_tr(:,:,:), G%Domain, scale=US%ppt_to_S) - call set_up_ALE_sponge_field(tmp_tr, G, GV, tv%S, ALE_CSp) + call set_up_ALE_sponge_field(tmp_tr, G, GV, tv%S, ALE_CSp, 'salt', & + sp_long_name='salinity', sp_unit='g kg-1 s-1') deallocate(tmp_tr) endif if (sponge_uv) then @@ -2160,8 +2162,10 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t endif ! The remaining calls to set_up_sponge_field can be in any order. if ( use_temperature) then - call set_up_ALE_sponge_field(filename, potemp_var, Time, G, GV, US, tv%T, ALE_CSp, scale=US%degC_to_C) - call set_up_ALE_sponge_field(filename, salin_var, Time, G, GV, US, tv%S, ALE_CSp, scale=US%ppt_to_S) + call set_up_ALE_sponge_field(filename, potemp_var, Time, G, GV, US, tv%T, ALE_CSp, & + 'temp', sp_long_name='temperature', sp_unit='degC s-1', scale=US%degC_to_C) + call set_up_ALE_sponge_field(filename, salin_var, Time, G, GV, US, tv%S, ALE_CSp, & + 'salt', sp_long_name='salinity', sp_unit='g kg-1 s-1', scale=US%ppt_to_S) endif if (sponge_uv) then filename = trim(inputdir)//trim(state_uv_file) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 1631a76dd6..9f5241bb9a 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -81,6 +81,9 @@ module MOM_ALE_sponge real :: scale = 1.0 !< A multiplicative factor by which to rescale input data real, dimension(:,:), pointer :: p => NULL() !< pointer the data. real, dimension(:,:), pointer :: h => NULL() !< pointer the data grid. + character(len=:), allocatable :: name !< The name of the input field + character(len=:), allocatable :: long_name !< The long name of the input field + character(len=:), allocatable :: unit !< The unit of the input field end type p2d !> ALE sponge control structure @@ -134,7 +137,7 @@ module MOM_ALE_sponge logical :: tripolar_N !< grid is folded at its north edge !>@{ Diagnostic IDs - integer, dimension(2) :: id_sp_tendency !< Diagnostic ids for temperature and salinity + integer, dimension(MAX_FIELDS_) :: id_sp_tendency !< Diagnostic ids for tracers !! tendency due to sponges integer :: id_sp_u_tendency !< Diagnostic id for zonal momentum tendency due to !! Rayleigh damping @@ -666,15 +669,19 @@ subroutine init_ALE_sponge_diags(Time, G, diag, CS, US) !! output. type(ALE_sponge_CS), intent(inout) :: CS !< ALE sponge control structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + ! Local Variables + integer :: m CS%diag => diag - CS%id_sp_tendency(1) = -1 - CS%id_sp_tendency(1) = register_diag_field('ocean_model', 'sp_tendency_temp', diag%axesTL, Time, & - 'Time tendency due to temperature restoring', 'degC s-1', conversion=US%s_to_T) - CS%id_sp_tendency(2) = -1 - CS%id_sp_tendency(2) = register_diag_field('ocean_model', 'sp_tendency_salt', diag%axesTL, Time, & - 'Time tendency due to salinity restoring', 'g kg-1 s-1', conversion=US%s_to_T) + do m=1,CS%fldno + CS%id_sp_tendency(m) = -1 + CS%id_sp_tendency(m) = register_diag_field('ocean_model', & + 'sp_tendency_' // CS%Ref_val(m)%name, diag%axesTL, Time, & + 'Time tendency due to restoring ' // CS%Ref_val(m)%long_name, & + CS%Ref_val(m)%unit, conversion=US%s_to_T) + enddo + CS%id_sp_u_tendency = -1 CS%id_sp_u_tendency = register_diag_field('ocean_model', 'sp_tendency_u', diag%axesCuL, Time, & 'Zonal acceleration due to sponges', 'm s-2', conversion=US%L_T2_to_m_s2) @@ -686,7 +693,8 @@ end subroutine init_ALE_sponge_diags !> This subroutine stores the reference profile at h points for the variable !! whose address is given by f_ptr. -subroutine set_up_ALE_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, scale) +subroutine set_up_ALE_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, & + sp_name, sp_long_name, sp_unit, scale) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(ALE_sponge_CS), pointer :: CS !< ALE sponge control structure (in/out). @@ -695,16 +703,27 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, scale) !! arbitrary number of layers. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & target, intent(in) :: f_ptr !< Pointer to the field to be damped + character(len=*), intent(in) :: sp_name !< The name of the tracer field + character(len=*), optional, & + intent(in) :: sp_long_name !< The long name of the tracer field + !! if not given, use the sp_name + character(len=*), optional, & + intent(in) :: sp_unit !< The unit of the tracer field + !! if not given, use the none real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any !! contributions due to dimensional rescaling. The default is 1. real :: scale_fac ! A factor by which to scale sp_val before storing it. integer :: k, col character(len=256) :: mesg ! String for error messages + character(len=256) :: long_name ! The long name of the tracer field + character(len=256) :: unit ! The unit of the tracer field if (.not.associated(CS)) return scale_fac = 1.0 ; if (present(scale)) scale_fac = scale + long_name = sp_name; if (present(sp_long_name)) long_name = sp_long_name + unit = 'none'; if (present(sp_unit)) unit = sp_unit CS%fldno = CS%fldno + 1 if (CS%fldno > MAX_FIELDS_) then @@ -716,6 +735,9 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, scale) ! stores the reference profile CS%Ref_val(CS%fldno)%nz_data = CS%nz_data + CS%Ref_val(CS%fldno)%name = sp_name + CS%Ref_val(CS%fldno)%long_name = long_name + CS%Ref_val(CS%fldno)%unit = unit allocate(CS%Ref_val(CS%fldno)%p(CS%nz_data,CS%num_col), source=0.0) do col=1,CS%num_col do k=1,CS%nz_data @@ -729,7 +751,8 @@ end subroutine set_up_ALE_sponge_field_fixed !> This subroutine stores the reference profile at h points for the variable !! whose address is given by filename and fieldname. -subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, f_ptr, CS, scale) +subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, f_ptr, CS, & + sp_name, sp_long_name, sp_unit, scale) character(len=*), intent(in) :: filename !< The name of the file with the !! time varying field data character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -741,6 +764,13 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & target, intent(in) :: f_ptr !< Pointer to the field to be damped (in). type(ALE_sponge_CS), pointer :: CS !< Sponge control structure (in/out). + character(len=*), intent(in) :: sp_name !< The name of the tracer field + character(len=*), optional, & + intent(in) :: sp_long_name !< The long name of the tracer field + !! if not given, use the sp_name + character(len=*), optional, & + intent(in) :: sp_unit !< The unit of the tracer field + !! if not given, use 'none' real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any !! contributions due to dimensional rescaling. The default is 1. @@ -749,6 +779,11 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, integer, dimension(4) :: fld_sz integer :: nz_data !< the number of vertical levels in this input field character(len=256) :: mesg ! String for error messages + character(len=256) :: long_name ! The long name of the tracer field + character(len=256) :: unit ! The unit of the tracer field + long_name = sp_name; if (present(sp_long_name)) long_name = sp_long_name + unit = 'none'; if (present(sp_unit)) unit = sp_unit + ! Local variables for ALE remapping if (.not.associated(CS)) return @@ -768,6 +803,9 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, else CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname) endif + CS%Ref_val(CS%fldno)%name = sp_name + CS%Ref_val(CS%fldno)%long_name = long_name + CS%Ref_val(CS%fldno)%unit = unit fld_sz(1:4) = -1 call get_external_field_info(CS%Ref_val(CS%fldno)%id, size=fld_sz) nz_data = fld_sz(3) @@ -1290,7 +1328,8 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, GV, turns, param_file) call rotate_array(sp_val_in, turns, sp_val) ! NOTE: This points sp_val with the unrotated field. See note below. - call set_up_ALE_sponge_field(sp_val, G, GV, sp_ptr, sponge) + call set_up_ALE_sponge_field(sp_val, G, GV, sp_ptr, sponge, & + sponge_in%Ref_val(n)%name, sp_long_name=sponge_in%Ref_val(n)%long_name, sp_unit=sponge_in%Ref_val(n)%unit) deallocate(sp_val_in) else diff --git a/src/tracer/RGC_tracer.F90 b/src/tracer/RGC_tracer.F90 index 6801269245..9ceadc602d 100644 --- a/src/tracer/RGC_tracer.F90 +++ b/src/tracer/RGC_tracer.F90 @@ -231,7 +231,7 @@ subroutine initialize_RGC_tracer(restart, day, G, GV, h, diag, OBC, CS, & do m=1,1 ! This is needed to force the compiler not to do a copy in the sponge calls. tr_ptr => CS%tr(:,:,:,m) - call set_up_ALE_sponge_field(temp, G, GV, tr_ptr, sponge_CSp) + call set_up_ALE_sponge_field(temp, G, GV, tr_ptr, sponge_CSp, 'RGC_tracer') enddo deallocate(temp) endif diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index d0ed88c128..a997cde26b 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -473,12 +473,10 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A enddo enddo ; enddo - if ( associated(tv%T) ) then - call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp) - endif - if ( associated(tv%S) ) then - call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp) - endif + if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp, 'temp', & + sp_long_name='temperature', sp_unit='degC s-1') + if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', & + sp_long_name='salinity', sp_unit='g kg-1 s-1') else diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 5e91eaa86a..aaededaa8c 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -626,12 +626,11 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, ! momentum is typically not damped within the sponge. ! ! The remaining calls to set_up_sponge_field can be in any order. ! - if ( associated(tv%T) ) then - call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp) - endif - if ( associated(tv%S) ) then - call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp) - endif + if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp, 'temp', & + sp_long_name='temperature', sp_unit='degC s-1') + if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', & + sp_long_name='salinity', sp_unit='g kg-1 s-1') + else ! layer mode ! 1) Read eta, salt and temp from IC file diff --git a/src/user/RGC_initialization.F90 b/src/user/RGC_initialization.F90 index d9c1846a0e..b56e0b895a 100644 --- a/src/user/RGC_initialization.F90 +++ b/src/user/RGC_initialization.F90 @@ -168,8 +168,10 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, use_ALE, C call initialize_ALE_sponge(Idamp, G, GV, PF, ACSp, h, nz) ! The remaining calls to set_up_sponge_field can be in any order. - if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp) - if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp) + if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp, 'temp', & + sp_long_name='temperature', sp_unit='degC s-1') + if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', & + sp_long_name='salinity', sp_unit='g kg-1 s-1') if (sponge_uv) then U1(:,:,:) = 0.0 ; V1(:,:,:) = 0.0 diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90 index fa44a78604..a81c400256 100644 --- a/src/user/dense_water_initialization.F90 +++ b/src/user/dense_water_initialization.F90 @@ -279,8 +279,10 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file, enddo enddo - if (associated(tv%T)) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp) - if (associated(tv%S)) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp) + if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp, 'temp', & + sp_long_name='temperature', sp_unit='degC s-1') + if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', & + sp_long_name='salinity', sp_unit='g kg-1 s-1') else call MOM_error(FATAL, "dense_water_initialize_sponges: trying to use non ALE sponge") endif diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index e4ce7e77f5..762477b6c4 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -446,12 +446,13 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_fil enddo endif enddo ; enddo - if (associated(tv%S)) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp) + if (associated(tv%S)) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', & + sp_long_name='salinity', sp_unit='g kg-1 s-1') else do j=G%jsc,G%jec ; do i=G%isc,G%iec eta(i,j,1) = 0.0 do k=2,nz - eta(i,j,k) = eta(i,j,k-1)- GV%H_to_Z * h_in(i,j,k-1) + eta(i,j,k) = eta(i,j,k-1) - GV%H_to_Z * h_in(i,j,k-1) enddo eta(i,j,nz+1) = -depth_tot(i,j) do k=1,nz @@ -469,4 +470,4 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_fil end subroutine dumbbell_initialize_sponges -end module dumbbell_initialization +end module dumbbell_initialization \ No newline at end of file From 9c5f77e9ef72c50b23adeeae7f3c556e4f4079cb Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 8 Aug 2022 20:02:58 -0400 Subject: [PATCH 02/15] Fix ice shelf initialization - Ice shelf needs to be initialized prior to ocean state initialization. This fixes any cases with an ice shelf using the FMScap. --- config_src/drivers/FMS_cap/ocean_model_MOM.F90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index b6bb14fc01..a12ab35240 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -53,6 +53,7 @@ module ocean_model_mod use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only: Update_Surface_Waves @@ -274,9 +275,13 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas if (.not.OS%is_ocean_pe) return OS%Time = Time_in ; OS%Time_dyn = Time_in + ! Call initialize MOM with an optional Ice Shelf CS which, if present triggers + ! initialization of ice shelf parameters and arrays. + call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, & OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & - diag_ptr=OS%diag, count_calls=.true., waves_CSp=OS%Waves) + diag_ptr=OS%diag, count_calls=.true., ice_shelf_CSp=OS%ice_shelf_CSp, & + waves_CSp=OS%Waves) call get_MOM_state_elements(OS%MOM_CSp, G=OS%grid, GV=OS%GV, US=OS%US, C_p=OS%C_p, & C_p_scaled=OS%fluxes%C_p, use_temp=use_temperature) @@ -372,9 +377,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas endif if (OS%use_ice_shelf) then - call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & - OS%diag, OS%forces, OS%fluxes) + call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes) + call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces) endif + if (OS%icebergs_alter_ocean) then call marine_ice_init(OS%Time, OS%grid, param_file, OS%diag, OS%marine_ice_CSp) if (.not. OS%use_ice_shelf) & From e3c71b070dce3da60d2b43cf55728da41a0a30ae Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 6 Aug 2022 11:28:34 -0400 Subject: [PATCH 03/15] +Add runtime parameter CHANNEL_DRAG_MAX_BBL_THICK Added the new runtime parameter CHANNEL_DRAG_MAX_BBL_THICK, to allow the CHANNEL_DRAG parameterization to use the full dynamically determined depth of the bottom boundary layer, or to use a fixed maximum thickness, regardless of the settings of other parameters, although for now the default value is set based on the settings of USE_JACKSON_PARAM and DRAG_AS_BODY_FORCE in order to reproduce existing solutions by default. There are new entries in some MOM_parameter_doc files. All answers in the MOM6-examples test suite are bitwise identical, and this test suite has been tested and works with revised values of this parameter. --- .../vertical/MOM_set_viscosity.F90 | 53 ++++++++++++++----- 1 file changed, 39 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 9bd995633f..448ae079f4 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -74,25 +74,27 @@ module MOM_set_visc !! the properties of the bottom boundary layer. logical :: linear_drag !< If true, the drag law is cdrag*`DRAG_BG_VEL`*u. !! Runtime parameter `LINEAR_DRAG`. - logical :: Channel_drag !< If true, the drag is exerted directly on each - !! layer according to what fraction of the bottom - !! they overlie. + logical :: Channel_drag !< If true, the drag is exerted directly on each layer + !! according to what fraction of the bottom they overlie. + real :: Chan_drag_max_vol !< The maximum bottom boundary layer volume within which the + !! channel drag is applied, normalized by the the full cell area, + !! or a negative value to apply no maximum [H ~> m or kg m-2]. logical :: correct_BBL_bounds !< If true, uses the correct bounds on the BBL thickness and !! viscosity so that the bottom layer feels the intended drag. logical :: RiNo_mix !< If true, use Richardson number dependent mixing. logical :: dynamic_viscous_ML !< If true, use a bulk Richardson number criterion to !! determine the mixed layer thickness for viscosity. real :: bulk_Ri_ML !< The bulk mixed layer used to determine the - !! thickness of the viscous mixed layer. Nondim. + !! thickness of the viscous mixed layer [nondim] real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. real :: ustar_min !< A minimum value of ustar to avoid numerical !! problems [Z T-1 ~> m s-1]. If the value is small enough, !! this should not affect the solution. real :: TKE_decay !< The ratio of the natural Ekman depth to the TKE - !! decay scale, nondimensional. - real :: omega_frac !< 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). + !! decay scale [nondim] + real :: omega_frac !< 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) [nondim] integer :: answer_date !< The vintage of the order of arithmetic and expressions in the set !! viscosity calculations. Values below 20190101 recover the answers !! from the end of 2018, while higher values use updated and more robust @@ -233,6 +235,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: a2x48_apb3, Iapb, Ibma_2 ! Combinations of a and slope [H-1 ~> m-1 or m2 kg-1]. ! All of the following "volumes" have units of thickness because they are normalized ! by the full horizontal area of a velocity cell. + real :: Vol_bbl_chan ! The volume of the bottom boundary layer as used in the channel + ! drag parameterization, normalized by the full horizontal area + ! of the velocity cell [H ~> m or kg m-2]. real :: Vol_open ! The cell volume above which it is open [H ~> m or kg m-2]. real :: Vol_direct ! With less than Vol_direct [H ~> m or kg m-2], there is a direct ! solution of a cubic equation for L. @@ -733,6 +738,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (bbl_thick < CS%BBL_thick_min) bbl_thick = CS%BBL_thick_min endif + ! Store the normalized bottom boundary layer volume. + if (CS%Channel_drag) Vol_bbl_chan = bbl_thick + ! If there is Richardson number dependent mixing, that determines ! the vertical extent of the bottom boundary layer, and there is no ! need to set that scale here. In fact, viscously reducing the @@ -746,10 +754,13 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (CS%body_force_drag) bbl_thick = h_bbl_drag(i) if (CS%Channel_drag) then - ! The drag within the bottommost bbl_thick is applied as a part of + ! The drag within the bottommost Vol_bbl_chan is applied as a part of ! an enhanced bottom viscosity, while above this the drag is applied ! directly to the layers in question as a Rayleigh drag term. + ! Restrict the volume over which the channel drag is applied. + if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol) + !### The harmonic mean edge depths here are not invariant to offsets! if (m==1) then D_vel = D_u(I,j) @@ -931,8 +942,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! Determine the drag contributing to the bottom boundary layer ! and the Rayleigh drag that acts on each layer. if (L(K) > L(K+1)) then - if (vol_below < bbl_thick) then - BBL_frac = (1.0-vol_below/bbl_thick)**2 + if (vol_below < Vol_bbl_chan) then + BBL_frac = (1.0-vol_below/Vol_bbl_chan)**2 BBL_visc_frac = BBL_visc_frac + BBL_frac*(L(K) - L(K+1)) else BBL_frac = 0.0 @@ -1957,6 +1968,8 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS real :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate that ! is used in place of the absolute value of the local Coriolis ! parameter in the denominator of some expressions [nondim] + real :: Chan_max_thick_dflt ! The default value for CHANNEL_DRAG_MAX_THICK [m] + real :: Z_rescale ! A rescaling factor for heights from the representation in ! a restart file to the internal representation in this run. real :: I_T_rescale ! A rescaling factor for time from the internal representation in this run @@ -2154,9 +2167,6 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS "The thickness over which near-surface velocities are "//& "averaged for the drag law under an ice shelf. By "//& "default this is the same as HBBL", units="m", default=CS%Hbbl, scale=GV%m_to_H) - ! These unit conversions are out outside the get_param calls because the are also defaults. - CS%Hbbl = CS%Hbbl * GV%m_to_H ! Rescale - CS%BBL_thick_min = CS%BBL_thick_min * GV%m_to_H ! Rescale call get_param(param_file, mdl, "KV", Kv_background, & "The background kinematic viscosity in the interior. "//& @@ -2195,9 +2205,24 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS if (CS%c_Smag < 0.0) CS%c_Smag = 0.15 endif + Chan_max_thick_dflt = -1.0 + if (CS%RiNo_mix) Chan_max_thick_dflt = 0.5*CS%Hbbl + if (CS%body_force_drag) Chan_max_thick_dflt = CS%Hbbl + call get_param(param_file, mdl, "CHANNEL_DRAG_MAX_BBL_THICK", CS%Chan_drag_max_vol, & + "The maximum bottom boundary layer thickness over which the channel drag is "//& + "exerted, or a negative value for no fixed limit, instead basing the BBL "//& + "thickness on the bottom stress, rotation and stratification. The default is "//& + "proportional to HBBL if USE_JACKSON_PARAM or DRAG_AS_BODY_FORCE is true.", & + units="m", default=Chan_max_thick_dflt, scale=GV%m_to_H, & + do_not_log=.not.CS%Channel_drag) + call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, & default=.false., do_not_log=.true.) + ! These unit conversions are out outside the get_param calls because they are also defaults. + CS%Hbbl = CS%Hbbl * GV%m_to_H ! Rescale + CS%BBL_thick_min = CS%BBL_thick_min * GV%m_to_H ! Rescale + if (CS%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then ! This is necessary for reproduciblity across restarts in non-symmetric mode. call pass_var(visc%Kv_shear_Bu, G%Domain, position=CORNER, complete=.true.) From 3db3882c9114d6766d7bd40579dc149c9c930f53 Mon Sep 17 00:00:00 2001 From: He Wang Date: Fri, 15 Apr 2022 21:15:11 -0400 Subject: [PATCH 04/15] Change porous_barrier_ptrs to allocatable * Porous barrier variables, which are grouped in porous_barrier_ptrs, are changed from pointers to allocatables. * Copies (aliases) of these variables in MOM_CS are removed. * The name porous_barrier_ptrs is yet to be changed to minimize the number of files needed to be edited. --- src/core/MOM.F90 | 21 +++++++-------------- src/core/MOM_variables.F90 | 15 ++++++++------- 2 files changed, 15 insertions(+), 21 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 78170064ff..94ce1cb39f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -412,12 +412,6 @@ module MOM !! increments and priors type(dbcomms_CS_type) :: dbcomms_CS !< Control structure for database client used for online ML/AI type(porous_barrier_ptrs) :: pbv !< porous barrier fractional cell metrics - real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: por_face_areaU !< fractional open area of U-faces [nondim] - real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: por_face_areaV !< fractional open area of V-faces [nondim] - real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: por_layer_widthU !< fractional open width - !! of U-faces [nondim] - real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: por_layer_widthV !< fractional open width - !! of V-faces [nondim] type(particles), pointer :: particles => NULL() ! NULL() !< a pointer to the stochastics control structure end type MOM_control_struct @@ -2478,12 +2472,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0 !allocate porous topography variables - ALLOC_(CS%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%por_face_areaU(:,:,:) = 1.0 - ALLOC_(CS%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%por_face_areaV(:,:,:) = 1.0 - ALLOC_(CS%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%por_layer_widthU(:,:,:) = 1.0 - ALLOC_(CS%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%por_layer_widthV(:,:,:) = 1.0 - CS%pbv%por_face_areaU => CS%por_face_areaU; CS%pbv%por_face_areaV=> CS%por_face_areaV - CS%pbv%por_layer_widthU => CS%por_layer_widthU; CS%pbv%por_layer_widthV => CS%por_layer_widthV + ALLOC_(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0 + ALLOC_(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%pbv%por_face_areaV(:,:,:) = 1.0 + ALLOC_(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%pbv%por_layer_widthU(:,:,:) = 1.0 + ALLOC_(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%pbv%por_layer_widthV(:,:,:) = 1.0 + ! Use the Wright equation of state by default, unless otherwise specified ! Note: this line and the following block ought to be in a separate ! initialization routine for tv. @@ -3775,8 +3768,8 @@ subroutine MOM_end(CS) if (CS%use_ALE_algorithm) call ALE_end(CS%ALE_CSp) !deallocate porous topography variables - DEALLOC_(CS%por_face_areaU) ; DEALLOC_(CS%por_face_areaV) - DEALLOC_(CS%por_layer_widthU) ; DEALLOC_(CS%por_layer_widthV) + DEALLOC_(CS%pbv%por_face_areaU) ; DEALLOC_(CS%pbv%por_face_areaV) + DEALLOC_(CS%pbv%por_layer_widthU) ; DEALLOC_(CS%pbv%por_layer_widthV) ! NOTE: Allocated in PressureForce_FV_Bouss if (associated(CS%tv%varT)) deallocate(CS%tv%varT) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index a6f9d79fe6..ef90b0b3c5 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -309,16 +309,17 @@ module MOM_variables type(group_pass_type) :: pass_FA_uv !< Structure for face area group halo updates end type BT_cont_type - -!> pointers to grids modifying cell metric at porous barriers +!> Container for grids modifying cell metric at porous barriers +! TODO: rename porous_barrier_ptrs to porous_barrier_type type, public :: porous_barrier_ptrs - real, pointer, dimension(:,:,:) :: por_face_areaU => NULL() !< fractional open area of U-faces [nondim] - real, pointer, dimension(:,:,:) :: por_face_areaV => NULL() !< fractional open area of V-faces [nondim] - real, pointer, dimension(:,:,:) :: por_layer_widthU => NULL() !< fractional open width of U-faces [nondim] - real, pointer, dimension(:,:,:) :: por_layer_widthV => NULL() !< fractional open width of V-faces [nondim] + ! Each of the following fields has nz layers. + real, allocatable :: por_face_areaU(:,:,:) !< fractional open area of U-faces [nondim] + real, allocatable :: por_face_areaV(:,:,:) !< fractional open area of V-faces [nondim] + ! Each of the following fields is found at nz+1 interfaces. + real, allocatable :: por_layer_widthU(:,:,:) !< fractional open width of U-faces [nondim] + real, allocatable :: por_layer_widthV(:,:,:) !< fractional open width of V-faces [nondim] end type porous_barrier_ptrs - contains !> Allocates the fields for the surface (return) properties of From e5cfaa990c756b4a3a0710ac8afddaffd2ed8d1e Mon Sep 17 00:00:00 2001 From: He Wang Date: Sat, 16 Apr 2022 15:57:07 -0400 Subject: [PATCH 05/15] Add a control structure for porous barrier * The interface porous_widths is deleted and subroutine por_widths is renamed as porous_widths. * porous_barrier_CS is added to control input and diagnostics of the module * Diagnostics for both interface and layer averages weights are added in subroutine porous_widths. * An _init subroutine is added to facilitate reading parameters and registering diagnostics * checksum debugs are added within subroutine porous_widths. --- src/core/MOM.F90 | 15 ++++--- src/core/MOM_porous_barriers.F90 | 68 ++++++++++++++++++++++++++++---- 2 files changed, 70 insertions(+), 13 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 94ce1cb39f..7e50be1de6 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -103,6 +103,7 @@ module MOM use MOM_open_boundary, only : open_boundary_register_restarts use MOM_open_boundary, only : update_segment_tracer_reservoirs use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init +use MOM_porous_barriers, only : porous_barrier_CS, porous_widths, porous_barriers_init use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS use MOM_set_visc, only : set_visc_init, set_visc_end @@ -141,8 +142,6 @@ module MOM use MOM_wave_interface, only : wave_parameters_CS, waves_end, waves_register_restarts use MOM_wave_interface, only : Update_Stokes_Drift -use MOM_porous_barriers, only : porous_widths - ! Database client used for machine-learning interface use MOM_database_comms, only : dbcomms_CS_type, database_comms_init, dbclient_type @@ -403,6 +402,8 @@ module MOM !< Pointer to the MOM diagnostics control structure type(offline_transport_CS), pointer :: offline_CSp => NULL() !< Pointer to the offline tracer transport control structure + type(porous_barrier_CS) :: por_bar_CS + !< Control structure for porous barrier logical :: ensemble_ocean !< if true, this run is part of a !! larger ensemble for the purpose of data assimilation @@ -1089,8 +1090,10 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call diag_update_remap_grids(CS%diag) endif - !update porous barrier fractional cell metrics - call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv) + ! Update porous barrier fractional cell metrics + call enable_averages(dt, Time_local, CS%diag) + call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv, CS%por_bar_CS) + call disable_averaging(CS%diag) ! The bottom boundary layer properties need to be recalculated. if (bbl_time_int > 0.0) then @@ -1417,7 +1420,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! and set_viscous_BBL is called as a part of the dynamic stepping. call cpu_clock_begin(id_clock_BBL_visc) !update porous barrier fractional cell metrics - call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv) + call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv, CS%por_bar_CS) call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM_thermo)") @@ -2815,6 +2818,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & new_sim = is_new_run(restart_CSp) call MOM_stoch_eos_init(G,Time,param_file,CS%stoch_eos_CS,restart_CSp,diag) + call porous_barriers_init(Time, US, param_file, diag, CS%por_bar_CS) + if (CS%split) then allocate(eta(SZI_(G),SZJ_(G)), source=0.0) call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, & diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index e807f19484..139baf0320 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -10,22 +10,30 @@ module MOM_porous_barriers use MOM_variables, only : thermo_var_ptrs, porous_barrier_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_interface_heights, only : find_eta +use MOM_time_manager, only : time_type +use MOM_diag_mediator, only : register_diag_field, diag_ctrl, post_data +use MOM_file_parser, only : param_file_type, get_param +use MOM_unit_scaling, only : unit_scale_type +use MOM_debugging, only : hchksum, uvchksum implicit none ; private -#include +public porous_widths, porous_barriers_init -public porous_widths +#include -!> Calculates curve fit from D_min, D_max, D_avg -interface porous_widths - module procedure por_widths, calc_por_layer -end interface porous_widths +type, public :: porous_barrier_CS; private + logical :: initialized = .false. !< True if this control structure has been initialized. + type(diag_ctrl), pointer :: diag => Null() !< A structure to regulate diagnostic output timing + logical :: debug !< If true, write verbose checksums for debugging purposes. + real :: mask_depth !< The depth below which porous barrier is not applied. + integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, id_por_face_areaU = -1, id_por_face_areaV = -1 +end type porous_barrier_CS contains !> subroutine to assign cell face areas and layer widths for porous topography -subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) +subroutine porous_widths(h, tv, G, GV, US, eta, pbv, CS, eta_bt, halo_size, eta_to_m) !eta_bt, halo_size, eta_to_m not currently used !variables needed to call find_eta type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -46,6 +54,7 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) real, optional, intent(in) :: eta_to_m !< The conversion factor from !! the units of eta to m; by default this is US%Z_to_m. type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_CS), intent(in) :: CS !local variables integer i, j, k, nk, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -54,6 +63,10 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) A_layer_prev, & ! integral of fractional open width from bottom to previous layer [Z ~> m] eta_s, & ! layer height used for fit [Z ~> m] eta_prev ! interface height of previous layer [Z ~> m] + + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_Porous_barrier: Module must be initialized before it is used.") + isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed IsdB = G%IsdB; IedB = G%IedB; JsdB = G%JsdB; JedB = G%JedB @@ -118,7 +131,19 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) endif enddo; enddo -end subroutine por_widths + if (CS%debug) then + call hchksum(eta, "Interface height used by porous barrier", G%HI, haloshift=0, scale=GV%H_to_m) + call uvchksum("Porous barrier weights at the layer-interface: por_layer_width[UV]", & + pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, haloshift=0) + call uvchksum("Porous barrier layer-averaged weights: por_face_area[UV]", & + pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0) + endif + + if (CS%id_por_layer_widthU > 0) call post_data(CS%id_por_layer_widthU, pbv%por_layer_widthU, CS%diag) + if (CS%id_por_layer_widthV > 0) call post_data(CS%id_por_layer_widthV, pbv%por_layer_widthV, CS%diag) + if (CS%id_por_face_areaU > 0) call post_data(CS%id_por_face_areaU, pbv%por_face_areaU, CS%diag) + if (CS%id_por_face_areaV > 0) call post_data(CS%id_por_face_areaV, pbv%por_face_areaV, CS%diag) +end subroutine porous_widths !> subroutine to calculate the profile fit for a single layer in a column subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) @@ -165,4 +190,31 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) end subroutine calc_por_layer +subroutine porous_barriers_init(Time, US, param_file, diag, CS) + type(porous_barrier_CS), intent(inout) :: CS + type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse + type(time_type), intent(in) :: Time !< Current model time + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(unit_scale_type), intent(in) :: US + + character(len=40) :: mdl = "MOM_porous_barriers" ! This module's name. + CS%initialized = .true. + CS%diag => diag + + call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) + call get_param(param_file, mdl, "POROUS_BARRIER_MASKING_DEPTH", CS%mask_depth, & + "The depth below which porous barrier is not applied. "//& + "This criterion is tested against TOPO_AT_VEL_VARNAME_U_AVE and TOPO_AT_VEL_VARNAME_V_AVE.", & + units="m", default=0.0, scale=US%m_to_Z) + + CS%id_por_layer_widthU = register_diag_field('ocean_model', 'por_layer_widthU', diag%axesCui, Time, & + 'Porous barrier open width fraction (at the layer interfaces) of the u-faces', 'nondim') + CS%id_por_layer_widthV = register_diag_field('ocean_model', 'por_layer_widthV', diag%axesCvi, Time, & + 'Porous barrier open width fraction (at the layer interfaces) of the v-faces', 'nondim') + CS%id_por_face_areaU = register_diag_field('ocean_model', 'por_face_areaU', diag%axesCuL, Time, & + 'Porous barrier open area fraction (layer averaged) of U-faces', 'nondim') + CS%id_por_face_areaV = register_diag_field('ocean_model', 'por_face_areaV', diag%axesCvL, Time, & + 'Porous barrier open area fraction (layer averaged) of V-faces', 'nondim') +end subroutine + end module MOM_porous_barriers From 00b48f13897bfd5a5ab7d1ae31ed5f864200afd1 Mon Sep 17 00:00:00 2001 From: He Wang Date: Wed, 20 Apr 2022 16:17:12 -0400 Subject: [PATCH 06/15] Fix porous_widths loop range This commit primarily fixes indexing bugs in subroutine porous_widths. * Loop range is subroutine porous_widths is changed from data domain to computation domain. * find_eta call is now with a proper halo to cover all velocity points. * Halo update for porous barrier fields is added in step_MOM_*. Other changes: * mask_depth, a component of porous_barrier_CS is now used to as the criterion for masking. This removes the limit that the cell edge depth has to be below sea surface. * Output variable eta_cor was unused and is now changed to a local. * Some unused indexing variables are removed. --- src/core/MOM.F90 | 15 +++++------ src/core/MOM_porous_barriers.F90 | 43 ++++++++++++++++---------------- 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7e50be1de6..ea5baf6116 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1051,8 +1051,6 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_por ! layer interface heights - !! for porous topo. [Z ~> m or 1/eta_to_m] G => CS%G ; GV => CS%GV ; US => CS%US ; IDs => CS%IDs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -1092,8 +1090,12 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & ! Update porous barrier fractional cell metrics call enable_averages(dt, Time_local, CS%diag) - call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv, CS%por_bar_CS) + call porous_widths(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) call disable_averaging(CS%diag) + call pass_vector(CS%pbv%por_face_areaU, CS%pbv%por_face_areaV, & + G%Domain, direction=To_All+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) + ! call pass_vector(CS%pbv%por_layer_widthU, CS%pbv%por_layer_widthV, & + ! G%Domain, direction=To_All+SCALAR_PAIR clock=id_clock_pass, halo=CS%cont_stencil) ! The bottom boundary layer properties need to be recalculated. if (bbl_time_int > 0.0) then @@ -1381,9 +1383,6 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & integer :: halo_sz ! The size of a halo where data must be valid. integer :: is, ie, js, je, nz - real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_por ! layer interface heights - !! for porous topo. [Z ~> m or 1/eta_to_m] - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("step_MOM_thermo(), MOM.F90") @@ -1420,7 +1419,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! and set_viscous_BBL is called as a part of the dynamic stepping. call cpu_clock_begin(id_clock_BBL_visc) !update porous barrier fractional cell metrics - call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv, CS%por_bar_CS) + call porous_widths(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) + call pass_vector(CS%pbv%por_layer_widthU, CS%pbv%por_layer_widthV, & + G%Domain, direction=To_ALL+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM_thermo)") diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 139baf0320..3de15d6336 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -33,17 +33,15 @@ module MOM_porous_barriers contains !> subroutine to assign cell face areas and layer widths for porous topography -subroutine porous_widths(h, tv, G, GV, US, eta, pbv, CS, eta_bt, halo_size, eta_to_m) +subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) !eta_bt, halo_size, eta_to_m not currently used !variables needed to call find_eta 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] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), 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. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta !< layer interface heights - !! [Z ~> m] or 1/eta_to_m m). real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic !! variable that gives the "correct" free surface height (Boussinesq) or total water !! column mass per unit area (non-Boussinesq). This is used to dilate the layer. @@ -57,30 +55,31 @@ subroutine porous_widths(h, tv, G, GV, US, eta, pbv, CS, eta_bt, halo_size, eta_ type(porous_barrier_CS), intent(in) :: CS !local variables - integer i, j, k, nk, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real w_layer, & ! fractional open width of layer interface [nondim] - A_layer, & ! integral of fractional open width from bottom to current layer[Z ~> m] - A_layer_prev, & ! integral of fractional open width from bottom to previous layer [Z ~> m] - eta_s, & ! layer height used for fit [Z ~> m] - eta_prev ! interface height of previous layer [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1):: eta !< layer interface heights [Z ~> m or 1/eta_to_m]. + integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq + real :: dmask + real :: w_layer, & ! fractional open width of layer interface [nondim] + A_layer, & ! integral of fractional open width from bottom to current layer[Z ~> m] + A_layer_prev, & ! integral of fractional open width from bottom to previous layer [Z ~> m] + eta_s, & ! layer height used for fit [Z ~> m] + eta_prev ! interface height of previous layer [Z ~> m] if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") - isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed - IsdB = G%IsdB; IedB = G%IedB; JsdB = G%JsdB; JedB = G%JedB + is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke + Isq = G%IscB; Ieq = G%IecB; Jsq = G%JscB; Jeq = G%JecB - !eta is zero at surface and decreases downward - - nk = SZK_(G) + dmask = CS%mask_depth + !eta is zero at surface and decreases downward !currently no treatment for using optional find_eta arguments if present - call find_eta(h, tv, G, GV, US, eta) + call find_eta(h, tv, G, GV, US, eta, halo_size=1) - do j=jsd,jed; do I=IsdB,IedB - if (G%porous_DavgU(I,j) < 0.) then + do j=js,je; do I=Isq,Ieq + if (G%porous_DavgU(I,j) < dmask) then do K = nk+1,1,-1 - eta_s = max(eta(I,j,K), eta(I+1,j,K)) !take shallower layer height + eta_s = max(eta(i,j,K), eta(i+1,j,K)) !take shallower layer height if (eta_s <= G%porous_DminU(I,j)) then pbv%por_layer_widthU(I,j,K) = 0.0 A_layer_prev = 0.0 @@ -104,10 +103,10 @@ subroutine porous_widths(h, tv, G, GV, US, eta, pbv, CS, eta_bt, halo_size, eta_ endif enddo; enddo - do J=JsdB,JedB; do i=isd,ied - if (G%porous_DavgV(i,J) < 0.) then + do J=Jsq,Jeq; do i=is,ie + if (G%porous_DavgV(i,J) < dmask) then do K = nk+1,1,-1 - eta_s = max(eta(i,J,K), eta(i,J+1,K)) !take shallower layer height + eta_s = max(eta(i,j,K), eta(i,j+1,K)) !take shallower layer height if (eta_s <= G%porous_DminV(i,J)) then pbv%por_layer_widthV(i,J,K) = 0.0 A_layer_prev = 0.0 From 6e91d0278e483eb42b37f24027ca6f64cb1dc9f2 Mon Sep 17 00:00:00 2001 From: He Wang Date: Wed, 20 Apr 2022 21:47:54 -0400 Subject: [PATCH 07/15] Add a method to read a full set of porous barrier parameters * subroutine set_subgrid_topo_at_vel_from_file is added to read max, min and avg depth at the edges from file. * The subroutine is only called when a new runtime parameter SUBGRID_TOPO_AT_VEL is True. Default is false. --- src/core/MOM_porous_barriers.F90 | 1 + .../MOM_fixed_initialization.F90 | 9 +++ .../MOM_shared_initialization.F90 | 73 +++++++++++++++++++ 3 files changed, 83 insertions(+) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 3de15d6336..112ed682c0 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -205,6 +205,7 @@ subroutine porous_barriers_init(Time, US, param_file, diag, CS) "The depth below which porous barrier is not applied. "//& "This criterion is tested against TOPO_AT_VEL_VARNAME_U_AVE and TOPO_AT_VEL_VARNAME_V_AVE.", & units="m", default=0.0, scale=US%m_to_Z) + CS%mask_depth = -CS%mask_depth CS%id_por_layer_widthU = register_diag_field('ocean_model', 'por_layer_widthU', diag%axesCui, Time, & 'Porous barrier open width fraction (at the layer interfaces) of the u-faces', 'nondim') diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index f0fb1d23f9..88c6377abc 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -23,6 +23,7 @@ module MOM_fixed_initialization use MOM_shared_initialization, only : set_rotation_planetary, set_rotation_beta_plane, initialize_grid_rotation_angle use MOM_shared_initialization, only : reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list use MOM_shared_initialization, only : read_face_length_list, set_velocity_depth_max, set_velocity_depth_min +use MOM_shared_initialization, only : set_subgrid_topo_at_vel_from_file use MOM_shared_initialization, only : compute_global_grid_integrals, write_ocean_geometry_file use MOM_unit_scaling, only : unit_scale_type @@ -62,6 +63,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) ! Local character(len=200) :: inputdir ! The directory where NetCDF input files are. character(len=200) :: config + logical :: read_porous_file character(len=40) :: mdl = "MOM_fixed_initialization" ! This module's name. logical :: debug ! This include declares and sets the variable "version". @@ -142,6 +144,13 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) end select endif + ! Read sub-grid scale topography parameters at velocity points used for porous barrier calculation + call get_param(PF, mdl, "SUBGRID_TOPO_AT_VEL", read_porous_file, & + "If true, use variables from TOPO_AT_VEL_FILE as parameters for porous barrier.", & + default=.False.) + if (read_porous_file) & + call set_subgrid_topo_at_vel_from_file(G, PF, US) + ! Calculate the value of the Coriolis parameter at the latitude ! ! of the q grid points [T-1 ~> s-1]. call MOM_initialize_rotation(G%CoriolisBu, G, PF, US=US) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 52f47f9581..04eab60569 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -27,6 +27,7 @@ module MOM_shared_initialization public set_rotation_planetary, set_rotation_beta_plane, initialize_grid_rotation_angle public reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list public read_face_length_list, set_velocity_depth_max, set_velocity_depth_min +public set_subgrid_topo_at_vel_from_file public compute_global_grid_integrals, write_ocean_geometry_file ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional @@ -1165,6 +1166,78 @@ subroutine read_face_length_list(iounit, filename, num_lines, lines) end subroutine read_face_length_list ! ----------------------------------------------------------------------------- +! ----------------------------------------------------------------------------- +!> Read from a file the maximum, minimum and average bathymetry at velocity points, +!! for the use of porous barrier. +!! Note that we assume the depth values in the sub-grid bathymetry file of the same +!! convention as in-cell bathymetry file, i.e. positive below the sea surface and +!! increasing downward. This is the opposite of the convention in subroutine +!! porous_widths. Therefore, all signs of the variable are reverted here. +subroutine set_subgrid_topo_at_vel_from_file(G, param_file, US) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + ! Local variables + character(len=200) :: filename, topo_file, inputdir ! Strings for file/path + character(len=200) :: varname_uhi, varname_ulo, varname_uav, & + varname_vhi, varname_vlo, varname_vav ! Variable names in file + character(len=40) :: mdl = "set_subgrid_topo_at_vel_from_file" ! This subroutine's name. + integer :: i, j + + call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") + + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + call get_param(param_file, mdl, "TOPO_AT_VEL_FILE", topo_file, & + "The file from which the bathymetry parameters at the velocity points are read. "//& + "While the names of the parameters reflect their physical locations, i.e. HIGH is above LOW, "//& + "their signs follow the model's convention, which is positive below the sea surface", & + default="topog_edge.nc") + call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_HIGH", varname_uhi, & + "The variable name of the highest bathymetry at the u-cells in TOPO_AT_VEL_FILE.", & + default="depthu_hi") + call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_LOW", varname_ulo, & + "The variable name of the lowest bathymetry at the u-cells in TOPO_AT_VEL_FILE.", & + default="depthu_lo") + call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_AVE", varname_uav, & + "The variable name of the average bathymetry at the u-cells in TOPO_AT_VEL_FILE.", & + default="depthu_av") + call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_HIGH", varname_vhi, & + "The variable name of the highest bathymetry at the v-cells in TOPO_AT_VEL_FILE.", & + default="depthv_hi") + call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_LOW", varname_vlo, & + "The variable name of the lowest bathymetry at the v-cells in TOPO_AT_VEL_FILE.", & + default="depthv_lo") + call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_AVE", varname_vav, & + "The variable name of the average bathymetry at the v-cells in TOPO_AT_VEL_FILE.", & + default="depthv_av") + + filename = trim(inputdir)//trim(topo_file) + call log_param(param_file, mdl, "INPUTDIR/TOPO_AT_VEL_FILE", filename) + + if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + " set_subgrid_topo_at_vel_from_file: Unable to open "//trim(filename)) + + call MOM_read_vector(filename, trim(varname_uhi), trim(varname_vhi), & + G%porous_DmaxU, G%porous_DmaxV, G%Domain, stagger=CGRID_NE, scale=US%m_to_Z) + call MOM_read_vector(filename, trim(varname_ulo), trim(varname_vlo), & + G%porous_DminU, G%porous_DminV, G%Domain, stagger=CGRID_NE, scale=US%m_to_Z) + call MOM_read_vector(filename, trim(varname_uav), trim(varname_vav), & + G%porous_DavgU, G%porous_DavgV, G%Domain, stagger=CGRID_NE, scale=US%m_to_Z) + + ! The signs of the depth parameters need to be reverted to comply with subroutine calc_por_layer, + ! which assumes depth is negative below the sea surface. + G%porous_DmaxU = -G%porous_DmaxU; G%porous_DminU = -G%porous_DminU; G%porous_DavgU = -G%porous_DavgU + G%porous_DmaxV = -G%porous_DmaxV; G%porous_DminV = -G%porous_DminV; G%porous_DavgV = -G%porous_DavgV + + call pass_vector(G%porous_DmaxU, G%porous_DmaxV, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + call pass_vector(G%porous_DminU, G%porous_DminV, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + call pass_vector(G%porous_DavgU, G%porous_DavgV, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + + call callTree_leave(trim(mdl)//'()') +end subroutine set_subgrid_topo_at_vel_from_file +! ----------------------------------------------------------------------------- + ! ----------------------------------------------------------------------------- !> Set the bathymetry at velocity points to be the maximum of the depths at the !! neighoring tracer points. From 0e25c5fa3aa238793c7102d043f4e1d91a85bf67 Mon Sep 17 00:00:00 2001 From: He Wang Date: Thu, 28 Apr 2022 11:09:48 -0400 Subject: [PATCH 08/15] Add timer for porous barrier * id_clock is added (as a private variable) to time porous barrier calculations. * Some small format fixes --- src/core/MOM_porous_barriers.F90 | 65 ++++++++++++++++++-------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 112ed682c0..2841f18248 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -4,17 +4,18 @@ module MOM_porous_barriers ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_error_handler, only : MOM_error, FATAL -use MOM_grid, only : ocean_grid_type -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs, porous_barrier_ptrs -use MOM_verticalGrid, only : verticalGrid_type +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE +use MOM_error_handler, only : MOM_error, FATAL +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs, porous_barrier_ptrs +use MOM_verticalGrid, only : verticalGrid_type use MOM_interface_heights, only : find_eta -use MOM_time_manager, only : time_type -use MOM_diag_mediator, only : register_diag_field, diag_ctrl, post_data -use MOM_file_parser, only : param_file_type, get_param -use MOM_unit_scaling, only : unit_scale_type -use MOM_debugging, only : hchksum, uvchksum +use MOM_time_manager, only : time_type +use MOM_diag_mediator, only : register_diag_field, diag_ctrl, post_data +use MOM_file_parser, only : param_file_type, get_param +use MOM_unit_scaling, only : unit_scale_type +use MOM_debugging, only : hchksum, uvchksum implicit none ; private @@ -30,6 +31,8 @@ module MOM_porous_barriers integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, id_por_face_areaU = -1, id_por_face_areaV = -1 end type porous_barrier_CS +integer :: id_clock_porous_barrier !< CPU clock for porous barrier + contains !> subroutine to assign cell face areas and layer widths for porous topography @@ -67,6 +70,8 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") + call cpu_clock_begin(id_clock_porous_barrier) + is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke Isq = G%IscB; Ieq = G%IecB; Jsq = G%JscB; Jeq = G%JecB @@ -142,24 +147,26 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) if (CS%id_por_layer_widthV > 0) call post_data(CS%id_por_layer_widthV, pbv%por_layer_widthV, CS%diag) if (CS%id_por_face_areaU > 0) call post_data(CS%id_por_face_areaU, pbv%por_face_areaU, CS%diag) if (CS%id_por_face_areaV > 0) call post_data(CS%id_por_face_areaV, pbv%por_face_areaV, CS%diag) + + call cpu_clock_end(id_clock_porous_barrier) end subroutine porous_widths !> subroutine to calculate the profile fit for a single layer in a column subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) - - real, intent(in) :: D_min !< minimum topographic height [Z ~> m] - real, intent(in) :: D_max !< maximum topographic height [Z ~> m] - real, intent(in) :: D_avg !< mean topographic height [Z ~> m] - real, intent(in) :: eta_layer !< height of interface [Z ~> m] - real, intent(out) :: w_layer !< frac. open interface width of current layer [nondim] - real, intent(out) :: A_layer !< frac. open face area of current layer [Z ~> m] - !local variables - real m, a, & !convenience constant for fit [nondim] - zeta, & !normalized vertical coordinate [nondim] - psi, & !fractional width of layer between D_min and D_max [nondim] - psi_int !integral of psi from 0 to zeta - - !three parameter fit from Adcroft 2013 + real, intent(in) :: D_min !< minimum topographic height [Z ~> m] + real, intent(in) :: D_max !< maximum topographic height [Z ~> m] + real, intent(in) :: D_avg !< mean topographic height [Z ~> m] + real, intent(in) :: eta_layer !< height of interface [Z ~> m] + real, intent(out) :: w_layer !< frac. open interface width of current layer [nondim] + real, intent(out) :: A_layer !< frac. open face area of current layer [Z ~> m] + + ! local variables + real :: m, a, & ! convenience constant for fit [nondim] + zeta, & ! normalized vertical coordinate [nondim] + psi, & ! fractional width of layer between D_min and D_max [nondim] + psi_int ! integral of psi from 0 to zeta + + ! three parameter fit from Adcroft 2013 m = (D_avg - D_min)/(D_max - D_min) a = (1. - m)/m @@ -185,18 +192,18 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) w_layer = psi A_layer = (D_max - D_min)*psi_int endif - - end subroutine calc_por_layer subroutine porous_barriers_init(Time, US, param_file, diag, CS) - type(porous_barrier_CS), intent(inout) :: CS + type(porous_barrier_CS), intent(inout) :: CS !< Module control structure type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse type(time_type), intent(in) :: Time !< Current model time type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure - type(unit_scale_type), intent(in) :: US + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + ! local variables character(len=40) :: mdl = "MOM_porous_barriers" ! This module's name. + CS%initialized = .true. CS%diag => diag @@ -215,6 +222,8 @@ subroutine porous_barriers_init(Time, US, param_file, diag, CS) 'Porous barrier open area fraction (layer averaged) of U-faces', 'nondim') CS%id_por_face_areaV = register_diag_field('ocean_model', 'por_face_areaV', diag%axesCvL, Time, & 'Porous barrier open area fraction (layer averaged) of V-faces', 'nondim') + + id_clock_porous_barrier = cpu_clock_id('(Ocean porous barrier)', grain=CLOCK_MODULE) end subroutine end module MOM_porous_barriers From bd4ae089f1073d53c22413b312cac4cf53673060 Mon Sep 17 00:00:00 2001 From: He Wang Date: Sun, 1 May 2022 23:33:59 -0400 Subject: [PATCH 09/15] Add options for various eta interpolation methods * A runtime parameter PORBAR_ETA_INTERP is added to control different methods for calculating interface height at the velocity points from adjacent tracer cells. * Two small thickness variables are added and scaled the unit of eta. This is to assit the harmonic mean calculation, but eventually the if-test eta_s - eta_prev>0 should be relaced by Angstrom. * Runtime parameter POROUS_BARRIER_MASKING_DEPTH is renamed to make it a liitlle bit shorter. * log_version call is added to separate out porous_barriers module in parameter file. --- src/core/MOM_porous_barriers.F90 | 82 +++++++++++++++++++++++++++++--- 1 file changed, 75 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 2841f18248..c1f8bf68ed 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -13,7 +13,7 @@ module MOM_porous_barriers use MOM_interface_heights, only : find_eta use MOM_time_manager, only : time_type use MOM_diag_mediator, only : register_diag_field, diag_ctrl, post_data -use MOM_file_parser, only : param_file_type, get_param +use MOM_file_parser, only : param_file_type, get_param, log_version use MOM_unit_scaling, only : unit_scale_type use MOM_debugging, only : hchksum, uvchksum @@ -28,11 +28,25 @@ module MOM_porous_barriers type(diag_ctrl), pointer :: diag => Null() !< A structure to regulate diagnostic output timing logical :: debug !< If true, write verbose checksums for debugging purposes. real :: mask_depth !< The depth below which porous barrier is not applied. + integer :: eta_interp !< An integer indicating how the interface heights at the velocity points + !! are calculated. Valid values are given by the parameters defined below: + !! MAX, MIN, ARITHMETIC and HARMONIC. integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, id_por_face_areaU = -1, id_por_face_areaV = -1 end type porous_barrier_CS integer :: id_clock_porous_barrier !< CPU clock for porous barrier +!>@{ Enumeration values for eta interpolation schemes +integer, parameter :: ETA_INTERP_MAX = 1 +integer, parameter :: ETA_INTERP_MIN = 2 +integer, parameter :: ETA_INTERP_ARITH = 3 +integer, parameter :: ETA_INTERP_HARM = 4 +character(len=20), parameter :: ETA_INTERP_MAX_STRING = "MAX" +character(len=20), parameter :: ETA_INTERP_MIN_STRING = "MIN" +character(len=20), parameter :: ETA_INTERP_ARITH_STRING = "ARITHMETIC" +character(len=20), parameter :: ETA_INTERP_HARM_STRING = "HARMONIC" +!>@} + contains !> subroutine to assign cell face areas and layer widths for porous topography @@ -66,6 +80,9 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) A_layer_prev, & ! integral of fractional open width from bottom to previous layer [Z ~> m] eta_s, & ! layer height used for fit [Z ~> m] eta_prev ! interface height of previous layer [Z ~> m] + real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. + real :: h_neglect, & ! ! Negligible thicknesses, often [Z ~> m] + h_min ! ! The minimum layer thickness, often [Z ~> m] if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") @@ -81,10 +98,26 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) !currently no treatment for using optional find_eta arguments if present call find_eta(h, tv, G, GV, US, eta, halo_size=1) + Z_to_eta = 1.0 ; if (present(eta_to_m)) Z_to_eta = US%Z_to_m / eta_to_m + H_to_eta = GV%H_to_m * US%m_to_Z * Z_to_eta + h_neglect = GV%H_subroundoff * H_to_eta + h_min = GV%Angstrom_H * H_to_eta + do j=js,je; do I=Isq,Ieq if (G%porous_DavgU(I,j) < dmask) then do K = nk+1,1,-1 - eta_s = max(eta(i,j,K), eta(i+1,j,K)) !take shallower layer height + select case (CS%eta_interp) + case (ETA_INTERP_MAX) !take shallower layer height + eta_s = max(eta(i,j,K), eta(i+1,j,K)) + case (ETA_INTERP_MIN) !take deeper layer height + eta_s = min(eta(i,j,K), eta(i+1,j,K)) + case (ETA_INTERP_ARITH) !take arithmetic mean layer height + eta_s = 0.5 * (eta(i,j,K) + eta(i+1,j,K)) + case (ETA_INTERP_HARM) !take harmonic mean layer height + eta_s = 2.0 * eta(i,j,K) * eta(i+1,j,K) / (eta(i,j,K) + eta(i+1,j,K) + h_neglect) + case default + call MOM_error(FATAL, "porous_widths: invalid value for eta interpolation method.") + end select if (eta_s <= G%porous_DminU(I,j)) then pbv%por_layer_widthU(I,j,K) = 0.0 A_layer_prev = 0.0 @@ -111,7 +144,18 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) do J=Jsq,Jeq; do i=is,ie if (G%porous_DavgV(i,J) < dmask) then do K = nk+1,1,-1 - eta_s = max(eta(i,j,K), eta(i,j+1,K)) !take shallower layer height + select case (CS%eta_interp) + case (ETA_INTERP_MAX) !take shallower layer height + eta_s = max(eta(i,j,K), eta(i,j+1,K)) + case (ETA_INTERP_MIN) !take deeper layer height + eta_s = min(eta(i,j,K), eta(i,j+1,K)) + case (ETA_INTERP_ARITH) !take arithmetic mean layer height + eta_s = 0.5 * (eta(i,j,K) + eta(i,j+1,K)) + case (ETA_INTERP_HARM) !take harmonic mean layer height + eta_s = 2.0 * eta(i,j,K) * eta(i,j+1,K) / (eta(i,j,K) + eta(i,j+1,K) + h_neglect) + case default + call MOM_error(FATAL, "porous_widths: invalid value for eta interpolation method.") + end select if (eta_s <= G%porous_DminV(i,J)) then pbv%por_layer_widthV(i,J,K) = 0.0 A_layer_prev = 0.0 @@ -201,18 +245,42 @@ subroutine porous_barriers_init(Time, US, param_file, diag, CS) type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + !> This include declares and sets the variable "version". +# include "version_variable.h" ! local variables - character(len=40) :: mdl = "MOM_porous_barriers" ! This module's name. + character(len=40) :: mdl = "MOM_porous_barriers" ! This module's name. + character(len=20) :: interp_method ! String storing eta interpolation method CS%initialized = .true. CS%diag => diag + call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.false., & + debugging=.false.) call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) - call get_param(param_file, mdl, "POROUS_BARRIER_MASKING_DEPTH", CS%mask_depth, & + call get_param(param_file, mdl, "PORBAR_MASKING_DEPTH", CS%mask_depth, & "The depth below which porous barrier is not applied. "//& - "This criterion is tested against TOPO_AT_VEL_VARNAME_U_AVE and TOPO_AT_VEL_VARNAME_V_AVE.", & - units="m", default=0.0, scale=US%m_to_Z) + "The effective average depths at the velocity cells are used "//& + "to test against this criterion.", units="m", default=0.0, & + scale=US%m_to_Z) CS%mask_depth = -CS%mask_depth + call get_param(param_file, mdl, "PORBAR_ETA_INTERP", interp_method, & + "A string describing the method that decicdes how the "//& + "interface heights at the velocity points are calculated. "//& + "Valid values are:\n"//& + "\t MAX (the default) - maximum of the adjacent cells \n"//& + "\t MIN - minimum of the adjacent cells \n"//& + "\t ARITHMETIC - arithmetic mean of the adjacent cells \n"//& + "\t HARMOINIC - harmonic mean of the adjacent cells \n", & + default=ETA_INTERP_MAX_STRING) ! do_not_log=.not.CS%use_por_bar) + select case (interp_method) + case (ETA_INTERP_MAX_STRING) ; CS%eta_interp = ETA_INTERP_MAX + case (ETA_INTERP_MIN_STRING) ; CS%eta_interp = ETA_INTERP_MIN + case (ETA_INTERP_ARITH_STRING) ; CS%eta_interp = ETA_INTERP_ARITH + case (ETA_INTERP_HARM_STRING) ; CS%eta_interp = ETA_INTERP_HARM + case default + call MOM_error(FATAL, "porous_barriers_init: Unrecognized setting "// & + "#define PORBAR_ETA_INTERP "//trim(interp_method)//" found in input file.") + end select CS%id_por_layer_widthU = register_diag_field('ocean_model', 'por_layer_widthU', diag%axesCui, Time, & 'Porous barrier open width fraction (at the layer interfaces) of the u-faces', 'nondim') From 8243071e9a9e0ec3ba7dee57e0f4bf26dcd2627d Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 2 May 2022 17:22:52 -0400 Subject: [PATCH 10/15] Re-arrange looping orders in porous_widths The main loops in porous_widths are simpified and cleaned up. * calc_por_layer iss slightly modified to reduced the if-blocks in the main loops. * k-loops are moved out. * To assist this new structure, two 2-D arrays are added to the stack. * A new function eta_at_uv is added to treat different interpolations. This is currently done at every point within the loop, and it is probably adding too much an overhead. It is better pre-calculate eta_[uv] all together, which would increase the stack size. --- src/core/MOM_porous_barriers.F90 | 176 ++++++++++++++----------------- 1 file changed, 82 insertions(+), 94 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index c1f8bf68ed..d1e9d65f15 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -72,17 +72,18 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) type(porous_barrier_CS), intent(in) :: CS !local variables - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1):: eta !< layer interface heights [Z ~> m or 1/eta_to_m]. - integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq - real :: dmask - real :: w_layer, & ! fractional open width of layer interface [nondim] - A_layer, & ! integral of fractional open width from bottom to current layer[Z ~> m] - A_layer_prev, & ! integral of fractional open width from bottom to previous layer [Z ~> m] - eta_s, & ! layer height used for fit [Z ~> m] - eta_prev ! interface height of previous layer [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1):: eta ! Layer interface heights [Z ~> m or 1/eta_to_m]. + real, dimension(SZIB_(G),SZJB_(G)) :: A_layer_prev ! Integral of fractional open width from the bottom + ! to the previous layer at u or v points [Z ~> m] + real, dimension(SZIB_(G),SZJB_(G)) :: eta_prev ! Layter interface height of the previous layer + ! at u or v points [Z ~> m] + real :: A_layer, & ! Integral of fractional open width from bottom to current layer [Z ~> m] + eta_s ! Layer height used for fit [Z ~> m] real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. real :: h_neglect, & ! ! Negligible thicknesses, often [Z ~> m] h_min ! ! The minimum layer thickness, often [Z ~> m] + real :: dmask ! The depth below which porous barrier is not applied. + integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") @@ -103,81 +104,49 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) h_neglect = GV%H_subroundoff * H_to_eta h_min = GV%Angstrom_H * H_to_eta - do j=js,je; do I=Isq,Ieq - if (G%porous_DavgU(I,j) < dmask) then - do K = nk+1,1,-1 - select case (CS%eta_interp) - case (ETA_INTERP_MAX) !take shallower layer height - eta_s = max(eta(i,j,K), eta(i+1,j,K)) - case (ETA_INTERP_MIN) !take deeper layer height - eta_s = min(eta(i,j,K), eta(i+1,j,K)) - case (ETA_INTERP_ARITH) !take arithmetic mean layer height - eta_s = 0.5 * (eta(i,j,K) + eta(i+1,j,K)) - case (ETA_INTERP_HARM) !take harmonic mean layer height - eta_s = 2.0 * eta(i,j,K) * eta(i+1,j,K) / (eta(i,j,K) + eta(i+1,j,K) + h_neglect) - case default - call MOM_error(FATAL, "porous_widths: invalid value for eta interpolation method.") - end select - if (eta_s <= G%porous_DminU(I,j)) then - pbv%por_layer_widthU(I,j,K) = 0.0 - A_layer_prev = 0.0 - if (K < nk+1) then - pbv%por_face_areaU(I,j,k) = 0.0; endif - else - call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), & - G%porous_DavgU(I,j), eta_s, w_layer, A_layer) - pbv%por_layer_widthU(I,j,K) = w_layer - if (k <= nk) then - if ((eta_s - eta_prev) > 0.0) then - pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/& - (eta_s-eta_prev) - else - pbv%por_face_areaU(I,j,k) = 0.0; endif - endif - eta_prev = eta_s - A_layer_prev = A_layer - endif - enddo + ! u-points + do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + eta_prev(I,j) = eta_at_uv(eta(i,j,nk+1), eta(i+1,j,nk+1), CS%eta_interp, h_neglect) + ! eta_prev(I,j) = max(eta(i,j,nk+1), eta(i+1,j,nk+1)) + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_prev(I,j), pbv%por_layer_widthU(I,j,nk+1), A_layer_prev(I,j)) + endif ; enddo ; enddo + + do k = nk,1,-1; do j=js,je; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + eta_s = eta_at_uv(eta(i,j,K), eta(i+1,j,K), CS%eta_interp, h_neglect) + ! eta_s = max(eta(i,j,K), eta(i+1,j,K)) + if ((eta_s - eta_prev(I,j)) > 0.0) then + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_s, pbv%por_layer_widthU(I,j,K), A_layer) + pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev(I,j)) / (eta_s - eta_prev(I,j)) + else + pbv%por_face_areaU(I,j,k) = 0.0 endif - enddo; enddo - - do J=Jsq,Jeq; do i=is,ie - if (G%porous_DavgV(i,J) < dmask) then - do K = nk+1,1,-1 - select case (CS%eta_interp) - case (ETA_INTERP_MAX) !take shallower layer height - eta_s = max(eta(i,j,K), eta(i,j+1,K)) - case (ETA_INTERP_MIN) !take deeper layer height - eta_s = min(eta(i,j,K), eta(i,j+1,K)) - case (ETA_INTERP_ARITH) !take arithmetic mean layer height - eta_s = 0.5 * (eta(i,j,K) + eta(i,j+1,K)) - case (ETA_INTERP_HARM) !take harmonic mean layer height - eta_s = 2.0 * eta(i,j,K) * eta(i,j+1,K) / (eta(i,j,K) + eta(i,j+1,K) + h_neglect) - case default - call MOM_error(FATAL, "porous_widths: invalid value for eta interpolation method.") - end select - if (eta_s <= G%porous_DminV(i,J)) then - pbv%por_layer_widthV(i,J,K) = 0.0 - A_layer_prev = 0.0 - if (K < nk+1) then - pbv%por_face_areaV(i,J,k) = 0.0; endif - else - call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), & - G%porous_DavgV(i,J), eta_s, w_layer, A_layer) - pbv%por_layer_widthV(i,J,K) = w_layer - if (k <= nk) then - if ((eta_s - eta_prev) > 0.0) then - pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/& - (eta_s-eta_prev) - else - pbv%por_face_areaU(I,j,k) = 0.0; endif - endif - eta_prev = eta_s - A_layer_prev = A_layer - endif - enddo + eta_prev(I,j) = eta_s + A_layer_prev(I,j) = A_layer + endif ; enddo ; enddo ; enddo + + ! v-points + do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + eta_prev(i,J) = eta_at_uv(eta(i,j,nk+1), eta(i,j+1,nk+1), CS%eta_interp, h_neglect) + ! eta_prev(i,J) = max(eta(i,j,nk+1), eta(i,j+1,nk+1)) + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_prev(i,J), pbv%por_layer_widthV(i,J,nk+1), A_layer_prev(i,J)) + endif ; enddo ; enddo + + do k = nk,1,-1; do J=Jsq,Jeq ; do i=is,ie ;if (G%porous_DavgV(i,J) < dmask) then + eta_s = eta_at_uv(eta(i,j,K), eta(i,j+1,K), CS%eta_interp, h_neglect) + ! eta_s = max(eta(i,j,K), eta(i,j+1,K)) + if ((eta_s - eta_prev(i,J)) > 0.0) then + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_s, pbv%por_layer_widthV(i,J,K), A_layer) + pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev(i,J)) / (eta_s - eta_prev(i,J)) + else + pbv%por_face_areaV(i,J,k) = 0.0 endif - enddo; enddo + eta_prev(i,J) = eta_s + A_layer_prev(i,J) = A_layer + endif ; enddo ; enddo ; enddo if (CS%debug) then call hchksum(eta, "Interface height used by porous barrier", G%HI, haloshift=0, scale=GV%H_to_m) @@ -195,7 +164,8 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) call cpu_clock_end(id_clock_porous_barrier) end subroutine porous_widths -!> subroutine to calculate the profile fit for a single layer in a column +!> subroutine to calculate the profile fit (the three parameter fit from Adcroft 2013) +! for a single layer in a column subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) real, intent(in) :: D_min !< minimum topographic height [Z ~> m] real, intent(in) :: D_max !< maximum topographic height [Z ~> m] @@ -210,34 +180,52 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) psi, & ! fractional width of layer between D_min and D_max [nondim] psi_int ! integral of psi from 0 to zeta - ! three parameter fit from Adcroft 2013 - m = (D_avg - D_min)/(D_max - D_min) - a = (1. - m)/m - - zeta = (eta_layer - D_min)/(D_max - D_min) - if (eta_layer <= D_min) then w_layer = 0.0 A_layer = 0.0 - elseif (eta_layer >= D_max) then + elseif (eta_layer > D_max) then w_layer = 1.0 A_layer = eta_layer - D_avg else + m = (D_avg - D_min) / (D_max - D_min) + a = (1.0 - m) / m + zeta = (eta_layer - D_min) / (D_max - D_min) if (m < 0.5) then - psi = zeta**(1./a) - psi_int = (1.-m)*zeta**(1./(1.-m)) + psi = zeta**(1.0 / a) + psi_int = (1.0 - m) * zeta**(1.0 / (1.0 - m)) elseif (m == 0.5) then psi = zeta - psi_int = 0.5*zeta*zeta + psi_int = 0.5 * zeta * zeta else - psi = 1. - (1. - zeta)**a - psi_int = zeta - m + m*((1-zeta)**(1/m)) + psi = 1.0 - (1.0 - zeta)**a + psi_int = zeta - m + m * ((1.0 - zeta)**(1 / m)) endif w_layer = psi A_layer = (D_max - D_min)*psi_int endif end subroutine calc_por_layer +function eta_at_uv(e1, e2, interp, h_neglect) result(eatuv) + real, intent(in) :: e1, e2 ! Interface heights at the adjacent tracer cells [Z ~> m] + real, intent(in) :: h_neglect ! Negligible thicknesses, often [Z ~> m] + integer, intent(in) :: interp ! Interpolation method coded by an integer + real :: eatuv + + select case (interp) + case (ETA_INTERP_MAX) ! The shallower interface height + eatuv = max(e1, e2) + case (ETA_INTERP_MIN) ! The deeper interface height + eatuv = min(e1, e2) + case (ETA_INTERP_ARITH) ! Arithmetic mean + eatuv = 0.5 * (e1 + e2) + case (ETA_INTERP_HARM) ! Harmonic mean + eatuv = 2.0 * e1 * e2 / (e1 + e2 + h_neglect) + case default + call MOM_error(FATAL, "porous_widths::eta_at_uv: "//& + "invalid value for eta interpolation method.") + end select +end function eta_at_uv + subroutine porous_barriers_init(Time, US, param_file, diag, CS) type(porous_barrier_CS), intent(inout) :: CS !< Module control structure type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse From 43502d4d65b1e7bf71b64109cc434bed18b456db Mon Sep 17 00:00:00 2001 From: He Wang Date: Tue, 3 May 2022 10:32:41 -0400 Subject: [PATCH 11/15] Split porous barrier interface and layer widths The averaged weight over a layer (por_face_area[UV]) and the weight at the interface (por_layer_width[UV]) are now calculated separately, as the only two updates in step_mom are in fact using only one of them. This reduces some unnecessary calculations. * Two new subroutines replaced the original `porous_widths`. They could be combined with interface if we remove porous_barrier_ptrs type, and use variables (of different nz) as output. * There is a place holder switch do_next in the two calc_por subs. This is used to further reduce calculations for all layers above D_max. Will be tested and implemented in the next commit. * A new subroutine calc_eta_at_uv is added to calculate eta at uv. This simplifies the code, but also increases stack and some performance overhead. --- src/core/MOM.F90 | 7 +- src/core/MOM_porous_barriers.F90 | 321 +++++++++++++++++++++---------- 2 files changed, 224 insertions(+), 104 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ea5baf6116..7865c7aa04 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -103,7 +103,8 @@ module MOM use MOM_open_boundary, only : open_boundary_register_restarts use MOM_open_boundary, only : update_segment_tracer_reservoirs use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init -use MOM_porous_barriers, only : porous_barrier_CS, porous_widths, porous_barriers_init +use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init +use MOM_porous_barriers, only : porous_barrier_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS use MOM_set_visc, only : set_visc_init, set_visc_end @@ -1090,7 +1091,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & ! Update porous barrier fractional cell metrics call enable_averages(dt, Time_local, CS%diag) - call porous_widths(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) + call porous_widths_layer(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) call disable_averaging(CS%diag) call pass_vector(CS%pbv%por_face_areaU, CS%pbv%por_face_areaV, & G%Domain, direction=To_All+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) @@ -1419,7 +1420,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! and set_viscous_BBL is called as a part of the dynamic stepping. call cpu_clock_begin(id_clock_BBL_visc) !update porous barrier fractional cell metrics - call porous_widths(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) + call porous_widths_interface(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) call pass_vector(CS%pbv%por_layer_widthU, CS%pbv%por_layer_widthV, & G%Domain, direction=To_ALL+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index d1e9d65f15..dff9325757 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -19,19 +19,20 @@ module MOM_porous_barriers implicit none ; private -public porous_widths, porous_barriers_init +public porous_widths_layer, porous_widths_interface, porous_barriers_init #include type, public :: porous_barrier_CS; private logical :: initialized = .false. !< True if this control structure has been initialized. type(diag_ctrl), pointer :: diag => Null() !< A structure to regulate diagnostic output timing - logical :: debug !< If true, write verbose checksums for debugging purposes. - real :: mask_depth !< The depth below which porous barrier is not applied. + logical :: debug !< If true, write verbose checksums for debugging purposes. + real :: mask_depth !< The depth below which porous barrier is not applied. integer :: eta_interp !< An integer indicating how the interface heights at the velocity points !! are calculated. Valid values are given by the parameters defined below: !! MAX, MIN, ARITHMETIC and HARMONIC. - integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, id_por_face_areaU = -1, id_por_face_areaV = -1 + integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, & + id_por_face_areaU = -1, id_por_face_areaV = -1 end type porous_barrier_CS integer :: id_clock_porous_barrier !< CPU clock for porous barrier @@ -49,41 +50,34 @@ module MOM_porous_barriers contains -!> subroutine to assign cell face areas and layer widths for porous topography -subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) +!> subroutine to assign porous barrier widths averaged over a layer +subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) !eta_bt, halo_size, eta_to_m not currently used !variables needed to call find_eta - 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_(GV)), 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. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic - !! variable that gives the "correct" free surface height (Boussinesq) or total water - !! column mass per unit area (non-Boussinesq). This is used to dilate the layer. - !! thicknesses when calculating interfaceheights [H ~> m or kg m-2]. - integer, optional, intent(in) :: halo_size !< width of halo points on - !! which to calculate eta. - - real, optional, intent(in) :: eta_to_m !< The conversion factor from - !! the units of eta to m; by default this is US%Z_to_m. - type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics - type(porous_barrier_CS), intent(in) :: CS + 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_(GV)), 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. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable + !! used to dilate the layer thicknesses + !! [H ~> m or kg m-2]. + type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_CS), intent(in) :: CS !< Control structure for porous barrier !local variables - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1):: eta ! Layer interface heights [Z ~> m or 1/eta_to_m]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_u ! Layer interface heights at u points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_v ! Layer interface heights at v points [Z ~> m] real, dimension(SZIB_(G),SZJB_(G)) :: A_layer_prev ! Integral of fractional open width from the bottom ! to the previous layer at u or v points [Z ~> m] - real, dimension(SZIB_(G),SZJB_(G)) :: eta_prev ! Layter interface height of the previous layer - ! at u or v points [Z ~> m] - real :: A_layer, & ! Integral of fractional open width from bottom to current layer [Z ~> m] - eta_s ! Layer height used for fit [Z ~> m] + real :: A_layer ! Integral of fractional open width from bottom to current layer [Z ~> m] real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. - real :: h_neglect, & ! ! Negligible thicknesses, often [Z ~> m] + real :: h_neglect, & ! Negligible thicknesses, often [Z ~> m] h_min ! ! The minimum layer thickness, often [Z ~> m] - real :: dmask ! The depth below which porous barrier is not applied. + real :: dmask ! The depth below which porous barrier is not applied [Z ~> m] integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq + logical :: do_next if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") @@ -95,136 +89,261 @@ subroutine porous_widths(h, tv, G, GV, US, pbv, CS, eta_bt, halo_size, eta_to_m) dmask = CS%mask_depth - !eta is zero at surface and decreases downward - !currently no treatment for using optional find_eta arguments if present - call find_eta(h, tv, G, GV, US, eta, halo_size=1) + call calc_eta_at_uv(eta_u, eta_v, CS%eta_interp, dmask, h, tv, G, GV, US) - Z_to_eta = 1.0 ; if (present(eta_to_m)) Z_to_eta = US%Z_to_m / eta_to_m + Z_to_eta = 1.0 H_to_eta = GV%H_to_m * US%m_to_Z * Z_to_eta - h_neglect = GV%H_subroundoff * H_to_eta h_min = GV%Angstrom_H * H_to_eta ! u-points do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then - eta_prev(I,j) = eta_at_uv(eta(i,j,nk+1), eta(i+1,j,nk+1), CS%eta_interp, h_neglect) - ! eta_prev(I,j) = max(eta(i,j,nk+1), eta(i+1,j,nk+1)) call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_prev(I,j), pbv%por_layer_widthU(I,j,nk+1), A_layer_prev(I,j)) + eta_u(I,j,nk+1), A_layer_prev(I,j), do_next) endif ; enddo ; enddo - do k = nk,1,-1; do j=js,je; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then - eta_s = eta_at_uv(eta(i,j,K), eta(i+1,j,K), CS%eta_interp, h_neglect) - ! eta_s = max(eta(i,j,K), eta(i+1,j,K)) - if ((eta_s - eta_prev(I,j)) > 0.0) then + do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + if ((eta_u(I,j,K) - eta_u(I,j,K+1)) > 0.0) then call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_s, pbv%por_layer_widthU(I,j,K), A_layer) - pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev(I,j)) / (eta_s - eta_prev(I,j)) + eta_u(I,j,K), A_layer, do_next) + pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1)) else pbv%por_face_areaU(I,j,k) = 0.0 endif - eta_prev(I,j) = eta_s A_layer_prev(I,j) = A_layer endif ; enddo ; enddo ; enddo ! v-points do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then - eta_prev(i,J) = eta_at_uv(eta(i,j,nk+1), eta(i,j+1,nk+1), CS%eta_interp, h_neglect) - ! eta_prev(i,J) = max(eta(i,j,nk+1), eta(i,j+1,nk+1)) call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_prev(i,J), pbv%por_layer_widthV(i,J,nk+1), A_layer_prev(i,J)) + eta_v(i,J,nk+1), A_layer_prev(i,J), do_next) endif ; enddo ; enddo - do k = nk,1,-1; do J=Jsq,Jeq ; do i=is,ie ;if (G%porous_DavgV(i,J) < dmask) then - eta_s = eta_at_uv(eta(i,j,K), eta(i,j+1,K), CS%eta_interp, h_neglect) - ! eta_s = max(eta(i,j,K), eta(i,j+1,K)) - if ((eta_s - eta_prev(i,J)) > 0.0) then + do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + if ((eta_v(i,J,K) - eta_v(i,J,K+1)) > 0.0) then call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_s, pbv%por_layer_widthV(i,J,K), A_layer) - pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev(i,J)) / (eta_s - eta_prev(i,J)) + eta_v(i,J,K), A_layer, do_next) + pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1)) else pbv%por_face_areaV(i,J,k) = 0.0 endif - eta_prev(i,J) = eta_s A_layer_prev(i,J) = A_layer endif ; enddo ; enddo ; enddo if (CS%debug) then - call hchksum(eta, "Interface height used by porous barrier", G%HI, haloshift=0, scale=GV%H_to_m) - call uvchksum("Porous barrier weights at the layer-interface: por_layer_width[UV]", & - pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, haloshift=0) + call uvchksum("Interface height used by porous barrier for layer weights", & + eta_u, eta_v, G%HI, haloshift=0) call uvchksum("Porous barrier layer-averaged weights: por_face_area[UV]", & - pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0) + pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0) endif - if (CS%id_por_layer_widthU > 0) call post_data(CS%id_por_layer_widthU, pbv%por_layer_widthU, CS%diag) - if (CS%id_por_layer_widthV > 0) call post_data(CS%id_por_layer_widthV, pbv%por_layer_widthV, CS%diag) if (CS%id_por_face_areaU > 0) call post_data(CS%id_por_face_areaU, pbv%por_face_areaU, CS%diag) if (CS%id_por_face_areaV > 0) call post_data(CS%id_por_face_areaV, pbv%por_face_areaV, CS%diag) call cpu_clock_end(id_clock_porous_barrier) -end subroutine porous_widths +end subroutine porous_widths_layer + +!> subroutine to assign porous barrier widths at the layer interfaces +subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) + !eta_bt, halo_size, eta_to_m not currently used + !variables needed to call find_eta + 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_(GV)), 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. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable + !! used to dilate the layer thicknesses + !! [H ~> m or kg m-2]. + type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_CS), intent(in) :: CS !< Control structure for porous barrier + + !local variables + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_u ! Layer interface height at u points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_v ! Layer interface height at v points [Z ~> m] + real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. + real :: h_neglect ! Negligible thicknesses, often [Z ~> m] + real :: dmask ! The depth below which porous barrier is not applied [Z ~> m] + integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq + logical :: do_next + + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_Porous_barrier: Module must be initialized before it is used.") + + call cpu_clock_begin(id_clock_porous_barrier) + + is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke + Isq = G%IscB; Ieq = G%IecB; Jsq = G%JscB; Jeq = G%JecB + + dmask = CS%mask_depth + + call calc_eta_at_uv(eta_u, eta_v, CS%eta_interp, dmask, h, tv, G, GV, US) + + ! u-points + do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_next) + endif ; enddo ; enddo ; enddo + + ! v-points + do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_next) + endif ; enddo ; enddo ; enddo + + if (CS%debug) then + call uvchksum("Interface height used by porous barrier for interface weights", & + eta_u, eta_v, G%HI, haloshift=0) + call uvchksum("Porous barrier weights at the layer-interface: por_layer_width[UV]", & + pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, haloshift=0) + endif + + if (CS%id_por_layer_widthU > 0) call post_data(CS%id_por_layer_widthU, pbv%por_layer_widthU, CS%diag) + if (CS%id_por_layer_widthV > 0) call post_data(CS%id_por_layer_widthV, pbv%por_layer_widthV, CS%diag) + + call cpu_clock_end(id_clock_porous_barrier) +end subroutine porous_widths_interface + +subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) + !variables needed to call find_eta + 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_(GV)), 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. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable + !! used to dilate the layer thicknesses + !! [H ~> m or kg m-2]. + real, intent(in) :: dmask !< The depth below which + !! porous barrier is not applied [Z ~> m] + integer, intent(in) :: interp !< eta interpolation method + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(out) :: eta_u !< Layer interface heights at u points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta_v !< Layer interface heights at v points [Z ~> m] + + ! local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Layer interface heights [Z ~> m or 1/eta_to_m]. + real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. + real :: h_neglect ! Negligible thicknesses, often [Z ~> m] + integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq + + is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke + Isq = G%IscB; Ieq = G%IecB; Jsq = G%JscB; Jeq = G%JecB + + ! currently no treatment for using optional find_eta arguments if present + call find_eta(h, tv, G, GV, US, eta, halo_size=1) + + Z_to_eta = 1.0 + H_to_eta = GV%H_to_m * US%m_to_Z * Z_to_eta + h_neglect = GV%H_subroundoff * H_to_eta + + select case (interp) + case (ETA_INTERP_MAX) ! The shallower interface height + do K=1,nk+1 + do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + eta_u(I,j,K) = max(eta(i,j,K), eta(i+1,j,K)) + endif ; enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + eta_v(i,J,K) = max(eta(i,j,K), eta(i,j+1,K)) + endif ; enddo ; enddo + enddo + case (ETA_INTERP_MIN) ! The deeper interface height + do K=1,nk+1 + do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + eta_u(I,j,K) = min(eta(i,j,K), eta(i+1,j,K)) + endif ; enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + eta_v(i,J,K) = min(eta(i,j,K), eta(i,j+1,K)) + endif ; enddo ; enddo + enddo + case (ETA_INTERP_ARITH) ! Arithmetic mean + do K=1,nk+1 + do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + eta_u(I,j,K) = 0.5 * (eta(i,j,K) + eta(i+1,j,K)) + endif ; enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + eta_v(i,J,K) = 0.5 * (eta(i,j,K) + eta(i,j+1,K)) + endif ; enddo ; enddo + enddo + case (ETA_INTERP_HARM) ! Harmonic mean + do K=1,nk+1 + do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + eta_u(I,j,K) = 2.0 * (eta(i,j,K) * eta(i+1,j,K)) / (eta(i,j,K) + eta(i+1,j,K) + h_neglect) + endif ; enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + eta_v(i,J,K) = 2.0 * (eta(i,j,K) * eta(i,j+1,K)) / (eta(i,j,K) + eta(i,j+1,K) + h_neglect) + endif ; enddo ; enddo + enddo + case default + call MOM_error(FATAL, "porous_widths::calc_eta_at_uv: "//& + "invalid value for eta interpolation method.") + end select +end subroutine calc_eta_at_uv !> subroutine to calculate the profile fit (the three parameter fit from Adcroft 2013) ! for a single layer in a column -subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) - real, intent(in) :: D_min !< minimum topographic height [Z ~> m] - real, intent(in) :: D_max !< maximum topographic height [Z ~> m] +subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, A_layer, do_next) + real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m] + real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m] real, intent(in) :: D_avg !< mean topographic height [Z ~> m] real, intent(in) :: eta_layer !< height of interface [Z ~> m] - real, intent(out) :: w_layer !< frac. open interface width of current layer [nondim] - real, intent(out) :: A_layer !< frac. open face area of current layer [Z ~> m] + real, intent(out) :: A_layer !< frac. open face area of below eta_layer [Z ~> m] + logical, intent(out) :: do_next !< False if eta_layer > D_max ! local variables - real :: m, a, & ! convenience constant for fit [nondim] - zeta, & ! normalized vertical coordinate [nondim] - psi, & ! fractional width of layer between D_min and D_max [nondim] - psi_int ! integral of psi from 0 to zeta + real :: m, & ! convenience constant for fit [nondim] + zeta ! normalized vertical coordinate [nondim] + do_next = .True. if (eta_layer <= D_min) then - w_layer = 0.0 A_layer = 0.0 elseif (eta_layer > D_max) then - w_layer = 1.0 A_layer = eta_layer - D_avg + do_next = .False. else m = (D_avg - D_min) / (D_max - D_min) - a = (1.0 - m) / m zeta = (eta_layer - D_min) / (D_max - D_min) if (m < 0.5) then - psi = zeta**(1.0 / a) - psi_int = (1.0 - m) * zeta**(1.0 / (1.0 - m)) + A_layer = (D_max - D_min) * ((1.0 - m) * zeta**(1.0 / (1.0 - m))) elseif (m == 0.5) then - psi = zeta - psi_int = 0.5 * zeta * zeta + A_layer = (D_max - D_min) * (0.5 * zeta * zeta) else - psi = 1.0 - (1.0 - zeta)**a - psi_int = zeta - m + m * ((1.0 - zeta)**(1 / m)) + A_layer = (D_max - D_min) * (zeta - m + m * ((1.0 - zeta)**(1.0 / m))) endif - w_layer = psi - A_layer = (D_max - D_min)*psi_int endif end subroutine calc_por_layer -function eta_at_uv(e1, e2, interp, h_neglect) result(eatuv) - real, intent(in) :: e1, e2 ! Interface heights at the adjacent tracer cells [Z ~> m] - real, intent(in) :: h_neglect ! Negligible thicknesses, often [Z ~> m] - integer, intent(in) :: interp ! Interpolation method coded by an integer - real :: eatuv +subroutine calc_por_interface(D_min, D_max, D_avg, eta_layer, w_layer, do_next) + real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m] + real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m] + real, intent(in) :: D_avg !< mean topographic height [Z ~> m] + real, intent(in) :: eta_layer !< height of interface [Z ~> m] + real, intent(out) :: w_layer !< frac. open interface width at eta_layer [nondim] + logical, intent(out) :: do_next !< False if eta_layer > D_max - select case (interp) - case (ETA_INTERP_MAX) ! The shallower interface height - eatuv = max(e1, e2) - case (ETA_INTERP_MIN) ! The deeper interface height - eatuv = min(e1, e2) - case (ETA_INTERP_ARITH) ! Arithmetic mean - eatuv = 0.5 * (e1 + e2) - case (ETA_INTERP_HARM) ! Harmonic mean - eatuv = 2.0 * e1 * e2 / (e1 + e2 + h_neglect) - case default - call MOM_error(FATAL, "porous_widths::eta_at_uv: "//& - "invalid value for eta interpolation method.") - end select -end function eta_at_uv + ! local variables + real :: m, a, & ! convenience constant for fit [nondim] + zeta ! normalized vertical coordinate [nondim] + + do_next = .True. + if (eta_layer <= D_min) then + w_layer = 0.0 + elseif (eta_layer > D_max) then + w_layer = 1.0 + do_next = .False. + else + m = (D_avg - D_min) / (D_max - D_min) + a = (1.0 - m) / m + zeta = (eta_layer - D_min) / (D_max - D_min) + if (m < 0.5) then + w_layer = zeta**(1.0 / a) + elseif (m == 0.5) then + w_layer = zeta + else + w_layer = 1.0 - (1.0 - zeta)**a + endif + endif +end subroutine calc_por_interface subroutine porous_barriers_init(Time, US, param_file, diag, CS) type(porous_barrier_CS), intent(inout) :: CS !< Module control structure From 016f042603215917b00ca21252c6ec0d3dd47f71 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 23 May 2022 16:38:37 -0400 Subject: [PATCH 12/15] Reduce unnecessary porous barrier calculations * A new runtime parameter USE_POROUS_BARRIER is added to control this feature. The default is True at the moment to assist potential test. It should be changed to false in the future. * A boolean array is added to avoid unnecessary porous barrier calculations above the shallowest points. * The layer-averaged weights are bounded by 1.0 explicitly, to avoid the cases it goes beyond due to roundoff errors. * For very thin layers (Angstrom), layer-averaged weights are set to zeros for simplicity. * Perhaps it is more reasonable to use the inferace weight for these cases. * It might be useful to make the thin layer definition a runtime parameter. * A bug related the size of eta_[uv] is fixed. This commit is potentially answer-changing when the porouss barrier module is used. --- src/core/MOM.F90 | 31 +++++++++------ src/core/MOM_porous_barriers.F90 | 68 ++++++++++++++++++++------------ 2 files changed, 62 insertions(+), 37 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7865c7aa04..8d2556f185 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -413,6 +413,7 @@ module MOM !! ensemble model state vectors and data assimilation !! increments and priors type(dbcomms_CS_type) :: dbcomms_CS !< Control structure for database client used for online ML/AI + logical :: use_porbar type(porous_barrier_ptrs) :: pbv !< porous barrier fractional cell metrics type(particles), pointer :: particles => NULL() ! NULL() !< a pointer to the stochastics control structure @@ -1090,13 +1091,13 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif ! Update porous barrier fractional cell metrics - call enable_averages(dt, Time_local, CS%diag) - call porous_widths_layer(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) - call disable_averaging(CS%diag) - call pass_vector(CS%pbv%por_face_areaU, CS%pbv%por_face_areaV, & - G%Domain, direction=To_All+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) - ! call pass_vector(CS%pbv%por_layer_widthU, CS%pbv%por_layer_widthV, & - ! G%Domain, direction=To_All+SCALAR_PAIR clock=id_clock_pass, halo=CS%cont_stencil) + if (CS%use_porbar) then + call enable_averages(dt, Time_local, CS%diag) + call porous_widths_layer(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) + call disable_averaging(CS%diag) + call pass_vector(CS%pbv%por_face_areaU, CS%pbv%por_face_areaV, & + G%Domain, direction=To_All+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) + endif ! The bottom boundary layer properties need to be recalculated. if (bbl_time_int > 0.0) then @@ -1420,9 +1421,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! and set_viscous_BBL is called as a part of the dynamic stepping. call cpu_clock_begin(id_clock_BBL_visc) !update porous barrier fractional cell metrics - call porous_widths_interface(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) - call pass_vector(CS%pbv%por_layer_widthU, CS%pbv%por_layer_widthV, & - G%Domain, direction=To_ALL+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) + if (CS%use_porbar) then + call porous_widths_interface(h, CS%tv, G, GV, US, CS%pbv, CS%por_bar_CS) + call pass_vector(CS%pbv%por_layer_widthU, CS%pbv%por_layer_widthV, & + G%Domain, direction=To_ALL+SCALAR_PAIR, clock=id_clock_pass, halo=CS%cont_stencil) + endif call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM_thermo)") @@ -1997,6 +2000,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "This is only used if THICKNESSDIFFUSE is true.", & default=.false.) if (.not.CS%thickness_diffuse) CS%thickness_diffuse_first = .false. + call get_param(param_file, "MOM", "USE_POROUS_BARRIER", CS%use_porbar, & + "If true, use porous barrier to constrain the widths "//& + " and face areas at the edges of the grid cells. ", & + default=.true.) ! The default should be false after tests. call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, & "If true, there are separate values for the basin depths "//& "at velocity points. Otherwise the effects of topography "//& @@ -2820,7 +2827,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & new_sim = is_new_run(restart_CSp) call MOM_stoch_eos_init(G,Time,param_file,CS%stoch_eos_CS,restart_CSp,diag) - call porous_barriers_init(Time, US, param_file, diag, CS%por_bar_CS) + + if (CS%use_porbar) & + call porous_barriers_init(Time, US, param_file, diag, CS%por_bar_CS) if (CS%split) then allocate(eta(SZI_(G),SZJ_(G)), source=0.0) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index dff9325757..b96c33d3f3 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -67,17 +67,18 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) type(porous_barrier_CS), intent(in) :: CS !< Control structure for porous barrier !local variables - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_u ! Layer interface heights at u points [Z ~> m] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_v ! Layer interface heights at v points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_u ! Layer interface heights at u points [Z ~> m] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_v ! Layer interface heights at v points [Z ~> m] real, dimension(SZIB_(G),SZJB_(G)) :: A_layer_prev ! Integral of fractional open width from the bottom ! to the previous layer at u or v points [Z ~> m] + logical, dimension(SZIB_(G),SZJB_(G)) :: do_I ! Booleans for calculation at u or v points + ! updated while moving up layers real :: A_layer ! Integral of fractional open width from bottom to current layer [Z ~> m] real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. real :: h_neglect, & ! Negligible thicknesses, often [Z ~> m] h_min ! ! The minimum layer thickness, often [Z ~> m] real :: dmask ! The depth below which porous barrier is not applied [Z ~> m] integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq - logical :: do_next if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") @@ -96,35 +97,39 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) h_min = GV%Angstrom_H * H_to_eta ! u-points + do j=js,je ; do I=Isq,Ieq ; do_I(I,j) = .False. ; enddo ; enddo + do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,nk+1), A_layer_prev(I,j), do_next) + eta_u(I,j,nk+1), A_layer_prev(I,j), do_I(I,j)) endif ; enddo ; enddo - do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then - if ((eta_u(I,j,K) - eta_u(I,j,K+1)) > 0.0) then - call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,K), A_layer, do_next) - pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1)) + do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), A_layer, do_I(I,j)) + if (eta_u(I,j,K) - (eta_u(I,j,K+1)+h_min) > 0.0) then + pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1))) else - pbv%por_face_areaU(I,j,k) = 0.0 + pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice endif A_layer_prev(I,j) = A_layer endif ; enddo ; enddo ; enddo ! v-points + do J=Jsq,Jeq ; do i=is,ie; do_I(i,J) = .False. ; enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,nk+1), A_layer_prev(i,J), do_next) + eta_v(i,J,nk+1), A_layer_prev(i,J), do_I(i,J)) endif ; enddo ; enddo - do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then - if ((eta_v(i,J,K) - eta_v(i,J,K+1)) > 0.0) then - call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,K), A_layer, do_next) - pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1)) + do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), A_layer, do_I(i,J)) + if (eta_v(i,J,K) - (eta_v(i,J,K+1)+h_min) > 0.0) then + pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1))) else - pbv%por_face_areaV(i,J,k) = 0.0 + pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice endif A_layer_prev(i,J) = A_layer endif ; enddo ; enddo ; enddo @@ -159,13 +164,14 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) type(porous_barrier_CS), intent(in) :: CS !< Control structure for porous barrier !local variables - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_u ! Layer interface height at u points [Z ~> m] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_v ! Layer interface height at v points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_u ! Layer interface height at u points [Z ~> m] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_v ! Layer interface height at v points [Z ~> m] + logical, dimension(SZIB_(G),SZJB_(G)) :: do_I ! Booleans for calculation at u or v points + ! updated while moving up layers real :: Z_to_eta, H_to_eta ! Unit conversion factors for eta. real :: h_neglect ! Negligible thicknesses, often [Z ~> m] real :: dmask ! The depth below which porous barrier is not applied [Z ~> m] integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq - logical :: do_next if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_Porous_barrier: Module must be initialized before it is used.") @@ -180,15 +186,25 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) call calc_eta_at_uv(eta_u, eta_v, CS%eta_interp, dmask, h, tv, G, GV, US) ! u-points - do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + do j=js,je ; do I=Isq,Ieq + do_I(I,j) = .False. + if (G%porous_DavgU(I,j) < dmask) do_I(I,j) = .True. + enddo ; enddo + + do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_next) + eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) endif ; enddo ; enddo ; enddo ! v-points - do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + do J=Jsq,Jeq ; do i=is,ie + do_I(i,J) = .False. + if (G%porous_DavgV(i,J) < dmask) do_I(i,J) = .True. + enddo ; enddo + + do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_next) + eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) endif ; enddo ; enddo ; enddo if (CS%debug) then @@ -218,8 +234,8 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) real, intent(in) :: dmask !< The depth below which !! porous barrier is not applied [Z ~> m] integer, intent(in) :: interp !< eta interpolation method - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(out) :: eta_u !< Layer interface heights at u points [Z ~> m] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta_v !< Layer interface heights at v points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta_u !< Layer interface heights at u points [Z ~> m] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(out) :: eta_v !< Layer interface heights at v points [Z ~> m] ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Layer interface heights [Z ~> m or 1/eta_to_m]. From 91131053c3ccb18b17d77240d1f1346a2d119e58 Mon Sep 17 00:00:00 2001 From: He Wang Date: Sat, 6 Aug 2022 17:28:52 -0400 Subject: [PATCH 13/15] Add PORBAR_ANSWER_DATE flag * The recently introduced ANSWER_DATE functionality is used to maintain a version that should be bit-for-bit reproducible for previous porous barrier related runs. * Rename porous_barrier_ptrs to porous_barrier_type * Some small documentation edits --- src/core/MOM.F90 | 6 +- src/core/MOM_CoriolisAdv.F90 | 4 +- src/core/MOM_continuity.F90 | 4 +- src/core/MOM_continuity_PPM.F90 | 4 +- src/core/MOM_dynamics_split_RK2.F90 | 6 +- src/core/MOM_dynamics_unsplit.F90 | 4 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 4 +- src/core/MOM_grid.F90 | 8 +- src/core/MOM_porous_barriers.F90 | 199 ++++++++++++------ src/core/MOM_variables.F90 | 6 +- src/framework/MOM_dyn_horgrid.F90 | 8 +- .../MOM_shared_initialization.F90 | 16 +- .../vertical/MOM_set_viscosity.F90 | 4 +- 13 files changed, 169 insertions(+), 104 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 8d2556f185..9f1fd99a26 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -134,7 +134,7 @@ module MOM use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state -use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_ptrs +use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_type use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state use MOM_variables, only : rotate_surface_state use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd @@ -414,7 +414,7 @@ module MOM !! increments and priors type(dbcomms_CS_type) :: dbcomms_CS !< Control structure for database client used for online ML/AI logical :: use_porbar - type(porous_barrier_ptrs) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type) :: pbv !< porous barrier fractional cell metrics type(particles), pointer :: particles => NULL() ! NULL() !< a pointer to the stochastics control structure end type MOM_control_struct @@ -2002,7 +2002,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (.not.CS%thickness_diffuse) CS%thickness_diffuse_first = .false. call get_param(param_file, "MOM", "USE_POROUS_BARRIER", CS%use_porbar, & "If true, use porous barrier to constrain the widths "//& - " and face areas at the edges of the grid cells. ", & + "and face areas at the edges of the grid cells. ", & default=.true.) ! The default should be false after tests. call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, & "If true, there are separate values for the basin depths "//& diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 6aacc479af..154db3eaa3 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -16,7 +16,7 @@ module MOM_CoriolisAdv use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : accel_diag_ptrs, porous_barrier_ptrs +use MOM_variables, only : accel_diag_ptrs, porous_barrier_type use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS @@ -140,7 +140,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv, Wav type(accel_diag_ptrs), intent(inout) :: AD !< Storage for acceleration diagnostics type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(CoriolisAdv_CS), intent(in) :: CS !< Control structure for MOM_CoriolisAdv - type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics type(Wave_parameters_CS), optional, pointer :: Waves !< An optional pointer to Stokes drift CS ! Local variables diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 0852d10cd2..541dcde66a 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -13,7 +13,7 @@ module MOM_continuity use MOM_grid, only : ocean_grid_type use MOM_open_boundary, only : ocean_OBC_type use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : BT_cont_type, porous_barrier_ptrs +use MOM_variables, only : BT_cont_type, porous_barrier_type use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -61,7 +61,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, v type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_CS), intent(in) :: CS !< Control structure for mom_continuity. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. - type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics real, dimension(SZIB_(G),SZJ_(G)), & optional, intent(in) :: uhbt !< The vertically summed volume !! flux through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1]. diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 402a6921ae..59d119f5d4 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -11,7 +11,7 @@ module MOM_continuity_PPM use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : BT_cont_type, porous_barrier_ptrs +use MOM_variables, only : BT_cont_type, porous_barrier_type use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -90,7 +90,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), intent(in) :: CS !< Module's control structure. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. - type(porous_barrier_ptrs), intent(in) :: pbv !< pointers to porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< pointers to porous barrier fractional cell metrics real, dimension(SZIB_(G),SZJ_(G)), & optional, intent(in) :: uhbt !< The summed volume flux through zonal faces !! [H L2 T-1 ~> m3 s-1 or kg s-1]. diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index c011d18c44..f06c692eaf 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -3,7 +3,7 @@ module MOM_dynamics_split_RK2 ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_type use MOM_variables, only : BT_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs use MOM_forcing_type, only : mech_forcing @@ -297,7 +297,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s type(MEKE_type), intent(inout) :: MEKE !< MEKE fields type(thickness_diffuse_CS), intent(inout) :: thickness_diffuse_CSp !< Pointer to a structure containing !! interface height diffusivities - type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -1152,7 +1152,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !! the number of times the velocity is !! truncated (this should be 0). logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step - type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics integer, intent(out) :: cont_stencil !< The stencil for thickness !! from the continuity solver. diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index a7517ccc4f..bc20c30a0f 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -50,7 +50,7 @@ module MOM_dynamics_unsplit !* * !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_type use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs use MOM_forcing_type, only : mech_forcing use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum @@ -217,7 +217,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & !! initialize_dyn_unsplit. type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure type(MEKE_type), intent(inout) :: MEKE !< MEKE fields - type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 9acb2b5c83..957306eb3d 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -48,7 +48,7 @@ module MOM_dynamics_unsplit_RK2 !* * !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_type use MOM_variables, only : ocean_internal_state, accel_diag_ptrs, cont_diag_ptrs use MOM_forcing_type, only : mech_forcing use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum @@ -230,7 +230,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, type(MEKE_type), intent(inout) :: MEKE !< MEKE fields !! fields related to the Mesoscale !! Eddy Kinetic Energy. - type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av ! Averaged layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted layer thicknesses [H ~> m or kg m-2] diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index f73fcc33af..384e4a6d2b 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -113,13 +113,13 @@ module MOM_grid areaCv !< The areas of the v-grid cells [L2 ~> m2]. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & - porous_DminU, & !< minimum topographic height of U-face [Z ~> m] - porous_DmaxU, & !< maximum topographic height of U-face [Z ~> m] + porous_DminU, & !< minimum topographic height (deepest) of U-face [Z ~> m] + porous_DmaxU, & !< maximum topographic height (shallowest) of U-face [Z ~> m] porous_DavgU !< average topographic height of U-face [Z ~> m] real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & - porous_DminV, & !< minimum topographic height of V-face [Z ~> m] - porous_DmaxV, & !< maximum topographic height of V-face [Z ~> m] + porous_DminV, & !< minimum topographic height (deepest) of V-face [Z ~> m] + porous_DmaxV, & !< maximum topographic height (shallowest) of V-face [Z ~> m] porous_DavgV !< average topographic height of V-face [Z ~> m] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index b96c33d3f3..39a2de5b7f 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -8,7 +8,7 @@ module MOM_porous_barriers use MOM_error_handler, only : MOM_error, FATAL use MOM_grid, only : ocean_grid_type use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs, porous_barrier_ptrs +use MOM_variables, only : thermo_var_ptrs, porous_barrier_type use MOM_verticalGrid, only : verticalGrid_type use MOM_interface_heights, only : find_eta use MOM_time_manager, only : time_type @@ -25,12 +25,16 @@ module MOM_porous_barriers type, public :: porous_barrier_CS; private logical :: initialized = .false. !< True if this control structure has been initialized. - type(diag_ctrl), pointer :: diag => Null() !< A structure to regulate diagnostic output timing - logical :: debug !< If true, write verbose checksums for debugging purposes. - real :: mask_depth !< The depth below which porous barrier is not applied. - integer :: eta_interp !< An integer indicating how the interface heights at the velocity points - !! are calculated. Valid values are given by the parameters defined below: - !! MAX, MIN, ARITHMETIC and HARMONIC. + type(diag_ctrl), pointer :: & + diag => Null() !< A structure to regulate diagnostic output timing + logical :: debug !< If true, write verbose checksums for debugging purposes. + real :: mask_depth !< The depth shallower than which porous barrier is not applied. + integer :: eta_interp !< An integer indicating how the interface heights at the velocity + !! points are calculated. Valid values are given by the parameters + !! defined below: MAX, MIN, ARITHMETIC and HARMONIC. + integer :: answer_date !< The vintage of the porous barrier weight function calculations. + !! Values below 20220806 recover the old answers in which the layer + !! averaged weights are not strictly limited by an upper-bound of 1.0 . integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, & id_por_face_areaU = -1, id_por_face_areaV = -1 end type porous_barrier_CS @@ -52,8 +56,7 @@ module MOM_porous_barriers !> subroutine to assign porous barrier widths averaged over a layer subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) - !eta_bt, halo_size, eta_to_m not currently used - !variables needed to call find_eta + ! Note: eta_bt is not currently used 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 @@ -63,7 +66,7 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable !! used to dilate the layer thicknesses !! [H ~> m or kg m-2]. - type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(inout) :: pbv !< porous barrier fractional cell metrics type(porous_barrier_CS), intent(in) :: CS !< Control structure for porous barrier !local variables @@ -88,7 +91,11 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke Isq = G%IscB; Ieq = G%IecB; Jsq = G%JscB; Jeq = G%JecB - dmask = CS%mask_depth + if (CS%answer_date < 20220806) then + dmask = 0.0 + else + dmask = CS%mask_depth + endif call calc_eta_at_uv(eta_u, eta_v, CS%eta_interp, dmask, h, tv, G, GV, US) @@ -104,16 +111,29 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) eta_u(I,j,nk+1), A_layer_prev(I,j), do_I(I,j)) endif ; enddo ; enddo - do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then - call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,K), A_layer, do_I(I,j)) - if (eta_u(I,j,K) - (eta_u(I,j,K+1)+h_min) > 0.0) then - pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1))) - else - pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice - endif - A_layer_prev(I,j) = A_layer - endif ; enddo ; enddo ; enddo + if (CS%answer_date < 20220806) then + do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), A_layer, do_I(I,j)) + if (eta_u(I,j,K) - eta_u(I,j,K+1) > 0.0) then + pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1)) + else + pbv%por_face_areaU(I,j,k) = 0.0 + endif + A_layer_prev(I,j) = A_layer + endif ; enddo ; enddo ; enddo + else + do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), A_layer, do_I(I,j)) + if (eta_u(I,j,K) - (eta_u(I,j,K+1)+h_min) > 0.0) then + pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1))) + else + pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice + endif + A_layer_prev(I,j) = A_layer + endif ; enddo ; enddo ; enddo + endif ! v-points do J=Jsq,Jeq ; do i=is,ie; do_I(i,J) = .False. ; enddo ; enddo @@ -123,16 +143,29 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) eta_v(i,J,nk+1), A_layer_prev(i,J), do_I(i,J)) endif ; enddo ; enddo - do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then - call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,K), A_layer, do_I(i,J)) - if (eta_v(i,J,K) - (eta_v(i,J,K+1)+h_min) > 0.0) then - pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1))) - else - pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice - endif - A_layer_prev(i,J) = A_layer - endif ; enddo ; enddo ; enddo + if (CS%answer_date < 20220806) then + do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), A_layer, do_I(i,J)) + if (eta_v(i,J,K) - eta_v(i,J,K+1) > 0.0) then + pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1)) + else + pbv%por_face_areaV(i,J,k) = 0.0 + endif + A_layer_prev(i,J) = A_layer + endif ; enddo ; enddo ; enddo + else + do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), A_layer, do_I(i,J)) + if (eta_v(i,J,K) - (eta_v(i,J,K+1)+h_min) > 0.0) then + pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1))) + else + pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice + endif + A_layer_prev(i,J) = A_layer + endif ; enddo ; enddo ; enddo + endif if (CS%debug) then call uvchksum("Interface height used by porous barrier for layer weights", & @@ -149,8 +182,7 @@ end subroutine porous_widths_layer !> subroutine to assign porous barrier widths at the layer interfaces subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) - !eta_bt, halo_size, eta_to_m not currently used - !variables needed to call find_eta + ! Note: eta_bt is not currently used 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 @@ -160,7 +192,7 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable !! used to dilate the layer thicknesses !! [H ~> m or kg m-2]. - type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type), intent(inout) :: pbv !< porous barrier fractional cell metrics type(porous_barrier_CS), intent(in) :: CS !< Control structure for porous barrier !local variables @@ -181,7 +213,11 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke Isq = G%IscB; Ieq = G%IecB; Jsq = G%JscB; Jeq = G%JecB - dmask = CS%mask_depth + if (CS%answer_date < 20220806) then + dmask = 0.0 + else + dmask = CS%mask_depth + endif call calc_eta_at_uv(eta_u, eta_v, CS%eta_interp, dmask, h, tv, G, GV, US) @@ -191,10 +227,17 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) if (G%porous_DavgU(I,j) < dmask) do_I(I,j) = .True. enddo ; enddo - do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then - call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) - endif ; enddo ; enddo ; enddo + if (CS%answer_date < 20220806) then + do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then + call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) + endif ; enddo ; enddo ; enddo + else + do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then + call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) + endif ; enddo ; enddo ; enddo + endif ! v-points do J=Jsq,Jeq ; do i=is,ie @@ -202,10 +245,17 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) if (G%porous_DavgV(i,J) < dmask) do_I(i,J) = .True. enddo ; enddo - do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then - call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) - endif ; enddo ; enddo ; enddo + if (CS%answer_date < 20220806) then + do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then + call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) + endif ; enddo ; enddo ; enddo + else + do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then + call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) + endif ; enddo ; enddo ; enddo + endif if (CS%debug) then call uvchksum("Interface height used by porous barrier for interface weights", & @@ -231,7 +281,7 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable !! used to dilate the layer thicknesses !! [H ~> m or kg m-2]. - real, intent(in) :: dmask !< The depth below which + real, intent(in) :: dmask !< The depth shaller than which !! porous barrier is not applied [Z ~> m] integer, intent(in) :: interp !< eta interpolation method real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta_u !< Layer interface heights at u points [Z ~> m] @@ -297,14 +347,14 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) end subroutine calc_eta_at_uv !> subroutine to calculate the profile fit (the three parameter fit from Adcroft 2013) -! for a single layer in a column +! of the open face area fraction below a certain depth (eta_layer) in a column subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, A_layer, do_next) - real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m] - real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m] - real, intent(in) :: D_avg !< mean topographic height [Z ~> m] - real, intent(in) :: eta_layer !< height of interface [Z ~> m] - real, intent(out) :: A_layer !< frac. open face area of below eta_layer [Z ~> m] - logical, intent(out) :: do_next !< False if eta_layer > D_max + real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m] + real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m] + real, intent(in) :: D_avg !< mean topographic height [Z ~> m] + real, intent(in) :: eta_layer !< height of interface [Z ~> m] + real, intent(out) :: A_layer !< frac. open face area of below eta_layer [Z ~> m] + logical, intent(out) :: do_next !< False if eta_layer>D_max ! local variables real :: m, & ! convenience constant for fit [nondim] @@ -329,13 +379,15 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, A_layer, do_next) endif end subroutine calc_por_layer +!> subroutine to calculate the profile fit (the three parameter fit from Adcroft 2013) +! of the open interface fraction at a certain depth (eta_layer) in a column subroutine calc_por_interface(D_min, D_max, D_avg, eta_layer, w_layer, do_next) - real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m] - real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m] - real, intent(in) :: D_avg !< mean topographic height [Z ~> m] - real, intent(in) :: eta_layer !< height of interface [Z ~> m] - real, intent(out) :: w_layer !< frac. open interface width at eta_layer [nondim] - logical, intent(out) :: do_next !< False if eta_layer > D_max + real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m] + real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m] + real, intent(in) :: D_avg !< mean topographic height [Z ~> m] + real, intent(in) :: eta_layer !< height of interface [Z ~> m] + real, intent(out) :: w_layer !< frac. open interface width at eta_layer [nondim] + logical, intent(out) :: do_next !< False if eta_layer>D_max ! local variables real :: m, a, & ! convenience constant for fit [nondim] @@ -362,29 +414,38 @@ subroutine calc_por_interface(D_min, D_max, D_avg, eta_layer, w_layer, do_next) end subroutine calc_por_interface subroutine porous_barriers_init(Time, US, param_file, diag, CS) - type(porous_barrier_CS), intent(inout) :: CS !< Module control structure - type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse - type(time_type), intent(in) :: Time !< Current model time - type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(porous_barrier_CS), intent(inout) :: CS !< Module control structure + type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse + type(time_type), intent(in) :: Time !< Current model time + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - !> This include declares and sets the variable "version". -# include "version_variable.h" ! local variables character(len=40) :: mdl = "MOM_porous_barriers" ! This module's name. character(len=20) :: interp_method ! String storing eta interpolation method + integer :: default_answer_date ! Global answer date + !> This include declares and sets the variable "version". +# include "version_variable.h" CS%initialized = .true. CS%diag => diag call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.false., & debugging=.false.) + call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & + "This sets the default value for the various _ANSWER_DATE parameters.", & + default=99991231) + call get_param(param_file, mdl, "PORBAR_ANSWER_DATE", CS%answer_date, & + "The vintage of the porous barrier weight function calculations. Values below "//& + "20220806 recover the old answers in which the layer averaged weights are not "//& + "strictly limited by an upper-bound of 1.0 .", default=default_answer_date) call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) call get_param(param_file, mdl, "PORBAR_MASKING_DEPTH", CS%mask_depth, & - "The depth below which porous barrier is not applied. "//& - "The effective average depths at the velocity cells are used "//& - "to test against this criterion.", units="m", default=0.0, & - scale=US%m_to_Z) + "If the effective average depth at the velocity cell is shallower than this "//& + "number, then porous barrier is not applied at that location. "//& + "PORBAR_MASKING_DEPTH is assumed to be positive below the sea surface.", & + units="m", default=0.0, scale=US%m_to_Z) + ! The sign needs to be inverted to be consistent with the sign convention of Davg_[UV] CS%mask_depth = -CS%mask_depth call get_param(param_file, mdl, "PORBAR_ETA_INTERP", interp_method, & "A string describing the method that decicdes how the "//& @@ -394,7 +455,7 @@ subroutine porous_barriers_init(Time, US, param_file, diag, CS) "\t MIN - minimum of the adjacent cells \n"//& "\t ARITHMETIC - arithmetic mean of the adjacent cells \n"//& "\t HARMOINIC - harmonic mean of the adjacent cells \n", & - default=ETA_INTERP_MAX_STRING) ! do_not_log=.not.CS%use_por_bar) + default=ETA_INTERP_MAX_STRING) select case (interp_method) case (ETA_INTERP_MAX_STRING) ; CS%eta_interp = ETA_INTERP_MAX case (ETA_INTERP_MIN_STRING) ; CS%eta_interp = ETA_INTERP_MIN diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index ef90b0b3c5..6fc81fa976 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -310,15 +310,15 @@ module MOM_variables end type BT_cont_type !> Container for grids modifying cell metric at porous barriers -! TODO: rename porous_barrier_ptrs to porous_barrier_type -type, public :: porous_barrier_ptrs +! TODO: rename porous_barrier_type to porous_barrier_type +type, public :: porous_barrier_type ! Each of the following fields has nz layers. real, allocatable :: por_face_areaU(:,:,:) !< fractional open area of U-faces [nondim] real, allocatable :: por_face_areaV(:,:,:) !< fractional open area of V-faces [nondim] ! Each of the following fields is found at nz+1 interfaces. real, allocatable :: por_layer_widthU(:,:,:) !< fractional open width of U-faces [nondim] real, allocatable :: por_layer_widthV(:,:,:) !< fractional open width of V-faces [nondim] -end type porous_barrier_ptrs +end type porous_barrier_type contains diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index bfc6f1b1a4..e97c8d981b 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -110,13 +110,13 @@ module MOM_dyn_horgrid areaCv !< The areas of the v-grid cells [L2 ~> m2]. real, allocatable, dimension(:,:) :: & - porous_DminU, & !< minimum topographic height of U-face [Z ~> m] - porous_DmaxU, & !< maximum topographic height of U-face [Z ~> m] + porous_DminU, & !< minimum topographic height (deepest) of U-face [Z ~> m] + porous_DmaxU, & !< maximum topographic height (shallowest) of U-face [Z ~> m] porous_DavgU !< average topographic height of U-face [Z ~> m] real, allocatable, dimension(:,:) :: & - porous_DminV, & !< minimum topographic height of V-face [Z ~> m] - porous_DmaxV, & !< maximum topographic height of V-face [Z ~> m] + porous_DminV, & !< minimum topographic height (deepest) of V-face [Z ~> m] + porous_DmaxV, & !< maximum topographic height (shallowest) of V-face [Z ~> m] porous_DavgV !< average topographic height of V-face [Z ~> m] real, allocatable, dimension(:,:) :: & diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 04eab60569..ff272e7fce 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1171,12 +1171,16 @@ end subroutine read_face_length_list !! for the use of porous barrier. !! Note that we assume the depth values in the sub-grid bathymetry file of the same !! convention as in-cell bathymetry file, i.e. positive below the sea surface and -!! increasing downward. This is the opposite of the convention in subroutine -!! porous_widths. Therefore, all signs of the variable are reverted here. +!! increasing downward; while in subroutine reset_face_lengths_list, it is implied +!! that read-in fields min_bathy, max_bathy and avg_bathy from the input file +!! CHANNEL_LIST_FILE all have negative values below the surface. Therefore, to ensure +!! backward compatibility, all signs of the variable are inverted here. +!! And porous_Dmax[UV] = shallowest point, porous_Dmin[UV] = deepest point subroutine set_subgrid_topo_at_vel_from_file(G, param_file, US) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + ! Local variables character(len=200) :: filename, topo_file, inputdir ! Strings for file/path character(len=200) :: varname_uhi, varname_ulo, varname_uav, & @@ -1225,8 +1229,8 @@ subroutine set_subgrid_topo_at_vel_from_file(G, param_file, US) call MOM_read_vector(filename, trim(varname_uav), trim(varname_vav), & G%porous_DavgU, G%porous_DavgV, G%Domain, stagger=CGRID_NE, scale=US%m_to_Z) - ! The signs of the depth parameters need to be reverted to comply with subroutine calc_por_layer, - ! which assumes depth is negative below the sea surface. + ! The signs of the depth parameters need to be inverted to be backward compatible with input files + ! used by subroutine reset_face_lengths_list, which assumes depth is negative below the sea surface. G%porous_DmaxU = -G%porous_DmaxU; G%porous_DminU = -G%porous_DminU; G%porous_DavgU = -G%porous_DavgU G%porous_DmaxV = -G%porous_DmaxV; G%porous_DminV = -G%porous_DminV; G%porous_DavgV = -G%porous_DavgV diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 448ae079f4..df6ed9063b 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -23,7 +23,7 @@ module MOM_set_visc use MOM_restart, only : register_restart_field_as_obsolete use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_ptrs +use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_type use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_DIRECTION_E @@ -139,7 +139,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) !! related fields. type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous !! call to set_visc_init. - type(porous_barrier_ptrs),intent(in) :: pbv !< porous barrier fractional cell metrics + type(porous_barrier_type),intent(in) :: pbv !< porous barrier fractional cell metrics ! Local variables real, dimension(SZIB_(G)) :: & From b8bf337894c49d72d243751aaf8a68743bb700cf Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 8 Aug 2022 00:42:01 -0400 Subject: [PATCH 14/15] Bugfix for uninitialized eta_[uv] * Fixed a bug that eta_[uv] was not properly initialized, which caused issues with global chksums with DEBUG=True. * Documents added to comply with Doxygen tests --- src/core/MOM.F90 | 3 ++- src/core/MOM_porous_barriers.F90 | 8 ++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9f1fd99a26..1cd4ea578e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -413,7 +413,8 @@ module MOM !! ensemble model state vectors and data assimilation !! increments and priors type(dbcomms_CS_type) :: dbcomms_CS !< Control structure for database client used for online ML/AI - logical :: use_porbar + logical :: use_porbar !< If true, use porous barrier to constrain the widths and face areas + !! at the edges of the grid cells. type(porous_barrier_type) :: pbv !< porous barrier fractional cell metrics type(particles), pointer :: particles => NULL() ! NULL() !< a pointer to the stochastics control structure diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 39a2de5b7f..cac8a44bf0 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -23,6 +23,7 @@ module MOM_porous_barriers #include +!> The control structure for the MOM_porous_barriers module type, public :: porous_barrier_CS; private logical :: initialized = .false. !< True if this control structure has been initialized. type(diag_ctrl), pointer :: & @@ -35,8 +36,10 @@ module MOM_porous_barriers integer :: answer_date !< The vintage of the porous barrier weight function calculations. !! Values below 20220806 recover the old answers in which the layer !! averaged weights are not strictly limited by an upper-bound of 1.0 . + !>@{ Diagnostic IDs integer :: id_por_layer_widthU = -1, id_por_layer_widthV = -1, & id_por_face_areaU = -1, id_por_face_areaV = -1 + !>@} end type porous_barrier_CS integer :: id_clock_porous_barrier !< CPU clock for porous barrier @@ -303,6 +306,11 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) H_to_eta = GV%H_to_m * US%m_to_Z * Z_to_eta h_neglect = GV%H_subroundoff * H_to_eta + do K=1,nk+1 + do j=js,je ; do I=Isq,Ieq ; eta_u(I,j,K) = dmask ; enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie ; eta_v(i,J,K) = dmask ; enddo ; enddo + enddo + select case (interp) case (ETA_INTERP_MAX) ! The shallower interface height do K=1,nk+1 From 34247a9b8276a6eea06b9f71abdb5b0dbb96d266 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 22 Aug 2022 14:25:00 -0400 Subject: [PATCH 15/15] Bugfix for allocating porous barrier variables --- src/core/MOM.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 1cd4ea578e..7444597bbb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2485,10 +2485,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0 !allocate porous topography variables - ALLOC_(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0 - ALLOC_(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%pbv%por_face_areaV(:,:,:) = 1.0 - ALLOC_(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%pbv%por_layer_widthU(:,:,:) = 1.0 - ALLOC_(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%pbv%por_layer_widthV(:,:,:) = 1.0 + allocate(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0 + allocate(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%pbv%por_face_areaV(:,:,:) = 1.0 + allocate(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%pbv%por_layer_widthU(:,:,:) = 1.0 + allocate(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%pbv%por_layer_widthV(:,:,:) = 1.0 ! Use the Wright equation of state by default, unless otherwise specified ! Note: this line and the following block ought to be in a separate @@ -3785,8 +3785,8 @@ subroutine MOM_end(CS) if (CS%use_ALE_algorithm) call ALE_end(CS%ALE_CSp) !deallocate porous topography variables - DEALLOC_(CS%pbv%por_face_areaU) ; DEALLOC_(CS%pbv%por_face_areaV) - DEALLOC_(CS%pbv%por_layer_widthU) ; DEALLOC_(CS%pbv%por_layer_widthV) + deallocate(CS%pbv%por_face_areaU) ; deallocate(CS%pbv%por_face_areaV) + deallocate(CS%pbv%por_layer_widthU) ; deallocate(CS%pbv%por_layer_widthV) ! NOTE: Allocated in PressureForce_FV_Bouss if (associated(CS%tv%varT)) deallocate(CS%tv%varT)