From 79e579713b9a326f6f5c1594024ee624b530159c Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 4 Jun 2020 13:06:26 -0400 Subject: [PATCH 01/21] First attempt to add a generic tracer to OBC scheme - This commit is the first attempt to add a (generic) tracer to the OBC. - The tracer "gtr1" is added from the ocean_bgc module generic_CFC.F90 - It is to imitate "salt" as if it was a passive tracer. --- src/core/MOM.F90 | 2 +- src/core/MOM_open_boundary.F90 | 131 ++++++++++++++++++++++++- src/tracer/MOM_generic_tracer.F90 | 17 ++-- src/tracer/MOM_tracer_flow_control.F90 | 11 ++- src/tracer/MOM_tracer_registry.F90 | 1 + 5 files changed, 149 insertions(+), 13 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 78d53e0b76..be2274182b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2133,7 +2133,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! This subroutine calls user-specified tracer registration routines. ! Additional calls can be added to MOM_tracer_flow_control.F90. call call_tracer_register(dG%HI, GV, US, param_file, CS%tracer_flow_CSp, & - CS%tracer_Reg, restart_CSp) + CS%tracer_Reg, restart_CSp, CS%OBC) call MEKE_alloc_register_restart(dG%HI, param_file, CS%MEKE, restart_CSp) call set_visc_register_restarts(dG%HI, GV, param_file, CS%visc, restart_CSp) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3b1559ab81..c6a2bb3cbc 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -53,7 +53,9 @@ module MOM_open_boundary public segment_tracer_registry_end public register_segment_tracer public register_temp_salt_segments +public register_obgc_segments public fill_temp_salt_segments +public fill_obgc_segments public open_boundary_register_restarts public update_segment_tracer_reservoirs public update_OBC_ramp @@ -1414,6 +1416,13 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) OBC%tracer_y_reservoirs_used(2) = .true. endif endif + if (fields(m) == 'gtr1') then + if (segment%is_E_or_W_2) then + OBC%tracer_x_reservoirs_used(2) = .true. + else + OBC%tracer_y_reservoirs_used(2) = .true. + endif + endif endif enddo ! Alternately, set first two to true if use_temperature is true @@ -3206,6 +3215,22 @@ function lookup_seg_field(OBC_seg,field) end function lookup_seg_field +!> Return the tracer index from its name +function get_tracer_index(OBC_seg,tr_name) + type(OBC_segment_type), pointer :: OBC_seg !< OBC segment + character(len=*), intent(in) :: tr_name !< The field name + integer :: get_tracer_index, it + get_tracer_index=-1 + it=1 + do while(allocated(OBC_seg%tr_Reg%Tr(it)%t)) + if (trim(OBC_seg%tr_Reg%Tr(it)%name) == trim(tr_name)) then + get_tracer_index=it + exit + endif + it=it+1 + enddo + return +end function get_tracer_index !> Allocate segment data fields subroutine allocate_OBC_segment_data(OBC, segment) @@ -3448,7 +3473,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) type(time_type), intent(in) :: Time !< Model time ! Local variables integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed - integer :: IsdB, IedB, JsdB, JedB, n, m, nz + integer :: IsdB, IedB, JsdB, JedB, n, m, nz, nt, it character(len=40) :: mdl = "set_OBC_segment_data" ! This subroutine's name. character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path type(OBC_segment_type), pointer :: segment => NULL() @@ -3937,6 +3962,25 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) else segment%tr_Reg%Tr(2)%OBC_inflow_conc = segment%field(m)%value endif + elseif (trim(segment%field(m)%name) == 'gtr1') then + nt=get_tracer_index(segment,'gtr1') + if(nt .lt. 0) then + call MOM_error(FATAL,"update_OBC_segment_data: Did not find tracer gtr1!") + endif + if (associated(segment%field(m)%buffer_dst)) then + do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc + segment%tr_Reg%Tr(nt)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k) + enddo ; enddo ; enddo + if (.not. segment%tr_Reg%Tr(nt)%is_initialized) then + !if the tracer reservoir has not yet been initialized, then set to external value. + do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc + segment%tr_Reg%Tr(nt)%tres(i,j,k) = segment%tr_Reg%Tr(nt)%t(i,j,k) + enddo ; enddo ; enddo + segment%tr_Reg%Tr(nt)%is_initialized=.true. + endif + else + segment%tr_Reg%Tr(nt)%OBC_inflow_conc = segment%field(m)%value + endif endif enddo ! end field loop @@ -4213,6 +4257,91 @@ subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file) end subroutine register_temp_salt_segments +subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name) + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry + type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values + character(len=*), intent(in) :: tr_name!< Tracer name +! Local variables + integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, nz, nf + integer :: i, j, k, n + type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list + type(tracer_type), pointer :: tr_ptr => NULL() + + if (.not. associated(OBC)) return + + do n=1, OBC%number_of_segments + segment=>OBC%segment(n) + if (.not. segment%on_pe) cycle + !For testing activate only one particular tracer for OBC + !This could be later generalized to all or a list of tracers + if(trim(tr_name) == 'gtr1') then + call tracer_name_lookup(tr_Reg, tr_ptr, tr_name) + call register_segment_tracer(tr_ptr, param_file, GV, segment, OBC_array=.True.) + endif + enddo + +end subroutine register_obgc_segments + +subroutine fill_obgc_segments(G, OBC, tr_ptr, tr_name) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + real, dimension(:,:,:), pointer :: tr_ptr !< Pointer to tracer field + character(len=*), intent(in) :: tr_name!< Tracer name + +! Local variables + integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n, nz, nt + integer :: i, j, k + type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list + + if (.not. associated(OBC)) return + + if(trim(tr_name) /= 'gtr1') return !Test for one particular tracer + + call pass_var(tr_ptr, G%Domain) + + nz = G%ke + + do n=1, OBC%number_of_segments + segment => OBC%segment(n) + if (.not. segment%on_pe) cycle + + nt=get_tracer_index(segment,tr_name) + if(nt .lt. 0) then + call MOM_error(FATAL,"fill_obgc_segments: Did not find tracer "// tr_name) + endif + + isd = segment%HI%isd ; ied = segment%HI%ied + jsd = segment%HI%jsd ; jed = segment%HI%jed + IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB + JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB + + ! Fill with Tracer values + if (segment%is_E_or_W) then + I=segment%HI%IsdB + do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed + if (segment%direction == OBC_DIRECTION_W) then + segment%tr_Reg%Tr(nt)%t(I,j,k) = tr_ptr(i+1,j,k) + else + segment%tr_Reg%Tr(nt)%t(I,j,k) = tr_ptr(i,j,k) + endif + enddo ; enddo + else + J=segment%HI%JsdB + do k=1,nz ; do i=segment%HI%isd,segment%HI%ied + if (segment%direction == OBC_DIRECTION_S) then + segment%tr_Reg%Tr(nt)%t(i,J,k) = tr_ptr(i,j+1,k) + else + segment%tr_Reg%Tr(nt)%t(i,J,k) = tr_ptr(i,j,k) + endif + enddo ; enddo + endif + segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) + enddo + call setup_OBC_tracer_reservoirs(G, OBC) !This will redo the T&S +end subroutine fill_obgc_segments + subroutine fill_temp_salt_segments(G, OBC, tv) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary structure diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 83c2c9a8e7..3bdb837b09 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -41,7 +41,7 @@ module MOM_generic_tracer use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs - use MOM_open_boundary, only : ocean_OBC_type + use MOM_open_boundary, only : ocean_OBC_type, register_obgc_segments, fill_obgc_segments use MOM_verticalGrid, only : verticalGrid_type @@ -69,7 +69,7 @@ module MOM_generic_tracer type(diag_ctrl), pointer :: diag => NULL() ! A structure that is used to ! regulate the timing of diagnostic output. type(MOM_restart_CS), pointer :: restart_CSp => NULL() - + type(ocean_OBC_type), pointer :: OBC => NULL() ! The following pointer will be directed to the first element of the ! linked list of generic tracers. type(g_tracer_type), pointer :: g_tracer_list => NULL() @@ -86,7 +86,7 @@ module MOM_generic_tracer !> Initializes the generic tracer packages and adds their tracers to the list !! Adds the tracers in the list of generic tracers to the set of MOM tracers (i.e., MOM-register them) !! Register these tracers for restart - function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) + function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, OBC) type(hor_index_type), intent(in) :: HI !< Horizontal index ranges type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters @@ -94,6 +94,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer !! advection and diffusion module. type(MOM_restart_CS), pointer :: restart_CS !< Pointer to the restart control structure. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. ! Local variables logical :: register_MOM_generic_tracer @@ -149,7 +151,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) "restart files of a restarted run.", default=.false.) CS%restart_CSp => restart_CS - + CS%OBC => OBC ntau=1 ! MOM needs the fields at only one time step @@ -195,6 +197,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) name=g_tracer_name, longname=longname, units=units, & registry_diags=.false., & !### CHANGE TO TRUE? restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + if (associated(CS%OBC)) & + call register_obgc_segments(GV, CS%OBC, tr_Reg, param_file, g_tracer_name) else call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & restart_CS, longname=longname, units=units) @@ -219,7 +223,7 @@ end function register_MOM_generic_tracer !! !! This subroutine initializes the NTR tracer fields in tr(:,:,:,:) !! and it sets up the tracer output. - subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, CS, & + subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, CS, & sponge_CSp, ALE_sponge_CSp) logical, intent(in) :: restart !< .true. if the fields have already been !! read from a restart file. @@ -230,8 +234,6 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(diag_ctrl), target, intent(in) :: diag !< Regulates diagnostic output. - type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, - !! where, and what open boundary conditions are used. type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. type(sponge_CS), pointer :: sponge_CSp !< Pointer to the control structure for the sponges. type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< Pointer to the control structure for the @@ -340,6 +342,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, endif endif + call fill_obgc_segments(G, CS%OBC, tr_ptr, g_tracer_name) !traverse the linked list till hit NULL call g_tracer_get_next(g_tracer, g_tracer_next) if (.NOT. associated(g_tracer_next)) exit diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 6e28477d26..d4c8df96d8 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -148,7 +148,7 @@ end subroutine call_tracer_flux_init !> This subroutine determines which tracer packages are to be used and does the calls to !! register their tracers to be advected, diffused, and read from restarts. -subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) +subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS, OBC) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -161,7 +161,10 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) !! advection and diffusion module. type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control !! structure. - + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition + !! type specifies whether, where, + !! and what open boundary + !! conditions are used. ! This include declares and sets the variable "version". # include "version_variable.h" @@ -256,7 +259,7 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) #ifdef _USE_GENERIC_TRACER if (CS%use_MOM_generic_tracer) CS%use_MOM_generic_tracer = & register_MOM_generic_tracer(HI, GV, param_file, CS%MOM_generic_tracer_CSp, & - tr_Reg, restart_CS) + tr_Reg, restart_CS, OBC) #endif if (CS%use_pseudo_salt_tracer) CS%use_pseudo_salt_tracer = & register_pseudo_salt_tracer(HI, GV, param_file, CS%pseudo_salt_tracer_CSp, & @@ -336,7 +339,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag sponge_CSp) #ifdef _USE_GENERIC_TRACER if (CS%use_MOM_generic_tracer) & - call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, & + call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, & CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp) #endif if (CS%use_pseudo_salt_tracer) & diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 07ca30dec8..9cadfa2159 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -847,6 +847,7 @@ subroutine tracer_name_lookup(Reg, tr_ptr, name) character(len=32), intent(in) :: name !< tracer name integer n + tr_ptr => null() do n=1,Reg%ntr if (lowercase(Reg%Tr(n)%name) == lowercase(name)) tr_ptr => Reg%Tr(n) enddo From 1638c0bf01e358d043475f4149dfa06cac2d8668 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 14 Aug 2020 18:15:20 -0400 Subject: [PATCH 02/21] Initialize OBC segments for OBGC tracers - This update queries the obgc modules (generic_tracers) for the OBC source files and var names for each generic tracers and initializes the segment fields accoringly. - With this update the obgc tracers should NOT appear in OBC_SEGMENT_XXX_DATA in MOM parameter files (MOM_override). - The default source file name is obgc_obc.nc - The default source file var name is the generic_tracer name - These can be overriden by field_table mechanism. E.g., "namelists","ocean_mod","generic_CFC" cfc11_obc_src_field_name = salt cfc12_obc_src_field_name = salt / --- src/core/MOM_open_boundary.F90 | 69 ++++++++++++++++++++++++++++--- src/tracer/MOM_generic_tracer.F90 | 13 +++--- 2 files changed, 71 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c786c58116..32c0f9a6cf 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -59,6 +59,7 @@ module MOM_open_boundary public register_obgc_segments public fill_temp_salt_segments public fill_obgc_segments +public set_obgc_segments_props public open_boundary_register_restarts public update_segment_tracer_reservoirs public update_OBC_ramp @@ -317,6 +318,15 @@ module MOM_open_boundary !! When locked=.true.,no more boundaries can be registered. end type OBC_registry_type +!> Type to carry OBC information needed for setting segments for OBGC tracers +type, private :: external_tracers_segments_props + type(external_tracers_segments_props), pointer :: next => NULL() + character(len=128) :: tracer_name + character(len=128) :: tracer_src_file + character(len=128) :: tracer_src_field +end type external_tracers_segments_props +type(external_tracers_segments_props), pointer, save :: obgc_segments_props => NULL() !< Linked-list of obgc tracers properties +integer, save :: num_obgc_tracers = 0 !< Keeps the total number of obgc tracers integer :: id_clock_pass !< A CPU time clock character(len=40) :: mdl = "MOM_open_boundary" !< This module's name. @@ -609,7 +619,7 @@ subroutine initialize_segment_data(G, OBC, PF) type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure type(param_file_type), intent(in) :: PF !< Parameter file handle - integer :: n,m,num_fields + integer :: n,m,num_fields,mm character(len=256) :: segstr, filename character(len=20) :: segnam, suffix character(len=32) :: varnam, fieldname @@ -625,6 +635,7 @@ subroutine initialize_segment_data(G, OBC, PF) integer, dimension(:), allocatable :: saved_pelist integer :: current_pe integer, dimension(1) :: single_pelist + type(external_tracers_segments_props), pointer :: obgc_segments_props_list =>NULL() !will be able to dynamically switch between sub-sampling refined grid data or model grid is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -677,8 +688,9 @@ subroutine initialize_segment_data(G, OBC, PF) cycle ! cycle to next segment endif - allocate(segment%field(num_fields)) - segment%num_fields = num_fields + !There are num_obgc_tracers obgc tracers are there that are not listed in param file + segment%num_fields = num_fields +num_obgc_tracers + allocate(segment%field(segment%num_fields)) segment%temp_segment_data_exists=.false. segment%salt_segment_data_exists=.false. @@ -691,8 +703,21 @@ subroutine initialize_segment_data(G, OBC, PF) IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB - do m=1,num_fields - call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname) + obgc_segments_props_list => obgc_segments_props !Get a pointer to the saved type + do m=1,segment%num_fields + if(m .le. num_fields) then + call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname) + else + call get_obgc_segments_props(obgc_segments_props_list,fields(m),filename,fieldname) + do mm=1,num_fields + if(trim(fields(m))==trim(fields(mm))) then + if (is_root_pe()) & + call MOM_error(FATAL,"MOM_open_boundary:initialize_segment_data(): obgc tracer " //trim(fields(m))// & +" appears in OBC_SEGMENT_XXX_DATA string in MOM6 param file. This is not supported!") + endif + enddo + endif + if (trim(filename) /= 'none') then OBC%update_OBC = .true. ! Data is assumed to be time-dependent if we are reading from file OBC%needs_IO_for_data = .true. ! At least one segment is using I/O for OBC data @@ -3297,7 +3322,7 @@ function get_tracer_index(OBC_seg,tr_name) integer :: get_tracer_index, it get_tracer_index=-1 it=1 - do while(allocated(OBC_seg%tr_Reg%Tr(it)%t)) + do while(associated(OBC_seg%tr_Reg%Tr(it)%t)) if (trim(OBC_seg%tr_Reg%Tr(it)%name) == trim(tr_name)) then get_tracer_index=it exit @@ -4395,6 +4420,38 @@ subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file) end subroutine register_temp_salt_segments +!> Sets the OBC properties of external obgc tracers, such as their source file and field name +subroutine set_obgc_segments_props(tr_name,obc_src_file_name,obc_src_field_name) + character(len=*), intent(in) :: tr_name !< Tracer name + character(len=*), intent(in) :: obc_src_file_name !< OBC source file name + character(len=*), intent(in) :: obc_src_field_name !< name of the field in the source file + + type(external_tracers_segments_props),pointer :: node_ptr => NULL() + allocate(node_ptr) + node_ptr%tracer_name = trim(tr_name) + node_ptr%tracer_src_file = trim(obc_src_file_name) + node_ptr%tracer_src_field = trim(obc_src_field_name) + !Reversed Linked List implementation! Make this new node to be the head of the list. + node_ptr%next => obgc_segments_props + obgc_segments_props => node_ptr + num_obgc_tracers = num_obgc_tracers+1 +end subroutine set_obgc_segments_props + +!> Get the OBC properties of external obgc tracers, such as their source file and field name +subroutine get_obgc_segments_props(node, tr_name,obc_src_file_name,obc_src_field_name) + type(external_tracers_segments_props),pointer :: node + character(len=*), intent(out) :: tr_name !< Tracer name + character(len=*), intent(out) :: obc_src_file_name !< OBC source file name + character(len=*), intent(out) :: obc_src_field_name !< name of the field in the source file + + tr_name=trim(node%tracer_name) + obc_src_file_name=trim(node%tracer_src_file) + obc_src_field_name=trim(node%tracer_src_field) + !pop the head. + node => node%next + +end subroutine get_obgc_segments_props + subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name) type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary structure diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 224a56f63c..c948e545fe 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -13,9 +13,6 @@ module MOM_generic_tracer #define _ALLOCATED allocated #endif - ! ### These imports should not reach into FMS directly ### - use field_manager_mod, only: fm_string_len - use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set @@ -27,7 +24,8 @@ module MOM_generic_tracer use g_tracer_utils, only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init use g_tracer_utils, only: g_tracer_send_diag,g_tracer_get_values use g_tracer_utils, only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag - + use g_tracer_utils, only: g_tracer_get_obc_segment_props + use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe @@ -47,7 +45,8 @@ module MOM_generic_tracer use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs - use MOM_open_boundary, only : ocean_OBC_type, register_obgc_segments, fill_obgc_segments + use MOM_open_boundary, only : ocean_OBC_type, register_obgc_segments, fill_obgc_segments + use MOM_open_boundary, only : set_obgc_segments_props use MOM_verticalGrid, only : verticalGrid_type @@ -56,6 +55,7 @@ module MOM_generic_tracer !> An state hidden in module data that is very much not allowed in MOM6 ! ### This needs to be fixed logical :: g_registered = .false. + integer, parameter :: fm_string_len=128 !>string lengths used in obgc packages public register_MOM_generic_tracer, initialize_MOM_generic_tracer public MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state @@ -118,6 +118,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, integer :: ntau, k,i,j,axes(3) type(g_tracer_type), pointer :: g_tracer,g_tracer_next character(len=fm_string_len) :: g_tracer_name,longname,units + character(len=fm_string_len) :: obc_src_file_name,obc_src_field_name real, dimension(:,:,:,:), pointer :: tr_field real, dimension(:,:,:), pointer :: tr_ptr real, dimension(HI%isd:HI%ied, HI%jsd:HI%jed,GV%ke) :: grid_tmask @@ -209,6 +210,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, registry_diags=.false., & !### CHANGE TO TRUE? restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) if (associated(CS%OBC)) & + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_src_file_name,obc_src_field_name) + call set_obgc_segments_props(g_tracer_name,obc_src_file_name,obc_src_field_name) call register_obgc_segments(GV, CS%OBC, tr_Reg, param_file, g_tracer_name) else call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & From d0c431b603141b3797b59659564e6c8e7b65a7bb Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 1 Sep 2020 10:42:55 -0400 Subject: [PATCH 03/21] Add OBC reservoirs for ocean_BGC tracers - This update adds OBC reservoirs for all obgc tracers that are registered to have OBC. - This fixes restart issue with obgc tracers in test runs --- src/core/MOM_open_boundary.F90 | 63 ++++++++++++++++--------------- src/tracer/MOM_generic_tracer.F90 | 16 +++++--- 2 files changed, 43 insertions(+), 36 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 32c0f9a6cf..c41068e1fa 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -83,6 +83,7 @@ module MOM_open_boundary integer :: fid !< handle from FMS associated with segment data on disk integer :: fid_dz !< handle from FMS associated with segment thicknesses on disk character(len=8) :: name !< a name identifier for the segment data + character(len=8) :: genre !< a family identifier for the segment data real, dimension(:,:,:), allocatable :: buffer_src !< buffer for segment data located at cell faces !! and on the original vertical grid integer :: nk_src !< Number of vertical levels in the source data @@ -708,6 +709,7 @@ subroutine initialize_segment_data(G, OBC, PF) if(m .le. num_fields) then call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname) else + segment%field(m)%genre='obgc' call get_obgc_segments_props(obgc_segments_props_list,fields(m),filename,fieldname) do mm=1,num_fields if(trim(fields(m))==trim(fields(mm))) then @@ -1372,7 +1374,7 @@ end function interpret_int_expr end subroutine parse_segment_str !> Parse an OBC_SEGMENT_%%%_DATA string - subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fields, num_fields, debug ) + subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fields, num_fields, debug, has_var ) character(len=*), intent(in) :: segment_str !< A string in form of !! "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..." character(len=*), optional, intent(in) :: var !< The name of the variable for which parameters are needed @@ -1384,14 +1386,16 @@ subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fi optional, intent(out) :: fields !< List of fieldnames for each segment integer, optional, intent(out) :: num_fields !< The number of fields in the segment data logical, optional, intent(in) :: debug !< If present and true, write verbose debugging messages + logical, optional, intent(out) :: has_var !< .true. if var is found in segment_str ! Local variables character(len=128) :: word1, word2, word3, method integer :: lword, nfields, n, m - logical :: continue,dbg + logical :: continue,dbg,has character(len=32), dimension(MAX_OBC_FIELDS) :: flds nfields=0 continue=.true. + has=.false. dbg=.false. if (PRESENT(debug)) dbg=debug @@ -1419,11 +1423,13 @@ subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fi do n=1,nfields if (trim(var)==trim(flds(n))) then m=n + has=.true. exit endif enddo + if (PRESENT(has_var)) has_var=has if (m==0) then - call abort() + return ! Why call abort() ? endif ! Process first word which will start with the fieldname @@ -1505,15 +1511,8 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) else OBC%tracer_y_reservoirs_used(2) = .true. endif - endif - if (fields(m) == 'gtr1') then - if (segment%is_E_or_W_2) then - OBC%tracer_x_reservoirs_used(2) = .true. - else - OBC%tracer_y_reservoirs_used(2) = .true. - endif - endif - endif + endif + endif enddo ! Alternately, set first two to true if use_temperature is true if (use_temperature) then @@ -1525,6 +1524,23 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) OBC%tracer_y_reservoirs_used(2) = .true. endif endif + !Add reservoirs for external/obgc tracers + !There is a diconnect in the above logic between tracer index and reservoir index. + !It arbitarily assigns reservoir indexes 1&2 to tracers T&S, + !so we need to start from 3 for the rest of tracers, hence the m+2 in the following. + !num_fields is the number of vars in segstr (6 of them now, U,V,SSH,TEMP,SALT,dye) + !but OBC%tracer_x_reservoirs_used is allocated to size Reg%ntr, which is the total number of tracers + !(t,s,dye,obgc1,obcg2,obgc3,... 6 of them by chance) + do m=1,num_obgc_tracers + !This logic assumes all external tarcers need a reservoir + !The segments for tracers are not initialized yet (that happens later in initialize_segment_data()) + !so we cannot query to determine if this tracer needs a reservoir. + if (segment%is_E_or_W_2) then + OBC%tracer_x_reservoirs_used(m+2) = .true. + else + OBC%tracer_y_reservoirs_used(m+2) = .true. + endif + enddo enddo return @@ -4125,10 +4141,10 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) else segment%tr_Reg%Tr(2)%OBC_inflow_conc = segment%field(m)%value endif - elseif (trim(segment%field(m)%name) == 'gtr1') then - nt=get_tracer_index(segment,'gtr1') + elseif (trim(segment%field(m)%genre) == 'obgc') then + nt=get_tracer_index(segment,trim(segment%field(m)%name)) if(nt .lt. 0) then - call MOM_error(FATAL,"update_OBC_segment_data: Did not find tracer gtr1!") + call MOM_error(FATAL,"update_OBC_segment_data: Did not find tracer "//trim(segment%field(m)%name)) endif if (associated(segment%field(m)%buffer_dst)) then do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc @@ -4469,12 +4485,8 @@ subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name) do n=1, OBC%number_of_segments segment=>OBC%segment(n) if (.not. segment%on_pe) cycle - !For testing activate only one particular tracer for OBC - !This could be later generalized to all or a list of tracers - if(trim(tr_name) == 'gtr1') then - call tracer_name_lookup(tr_Reg, tr_ptr, tr_name) - call register_segment_tracer(tr_ptr, param_file, GV, segment, OBC_array=.True.) - endif + call tracer_name_lookup(tr_Reg, tr_ptr, tr_name) + call register_segment_tracer(tr_ptr, param_file, GV, segment, OBC_array=.True.) enddo end subroutine register_obgc_segments @@ -4484,34 +4496,25 @@ subroutine fill_obgc_segments(G, OBC, tr_ptr, tr_name) type(ocean_OBC_type), pointer :: OBC !< Open boundary structure real, dimension(:,:,:), pointer :: tr_ptr !< Pointer to tracer field character(len=*), intent(in) :: tr_name!< Tracer name - ! Local variables integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n, nz, nt integer :: i, j, k type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list if (.not. associated(OBC)) return - - if(trim(tr_name) /= 'gtr1') return !Test for one particular tracer - call pass_var(tr_ptr, G%Domain) - nz = G%ke - do n=1, OBC%number_of_segments segment => OBC%segment(n) if (.not. segment%on_pe) cycle - nt=get_tracer_index(segment,tr_name) if(nt .lt. 0) then call MOM_error(FATAL,"fill_obgc_segments: Did not find tracer "// tr_name) endif - isd = segment%HI%isd ; ied = segment%HI%ied jsd = segment%HI%jsd ; jed = segment%HI%jed IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB - ! Fill with Tracer values if (segment%is_E_or_W) then I=segment%HI%IsdB diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index c948e545fe..0de68fa6ed 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -110,7 +110,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, ! Local variables logical :: register_MOM_generic_tracer - + logical :: obc_has character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer' character(len=200) :: inputdir ! The directory where NetCDF input files are. ! These can be overridden later in via the field manager? @@ -210,9 +210,12 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, registry_diags=.false., & !### CHANGE TO TRUE? restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) if (associated(CS%OBC)) & - call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_src_file_name,obc_src_field_name) - call set_obgc_segments_props(g_tracer_name,obc_src_file_name,obc_src_field_name) - call register_obgc_segments(GV, CS%OBC, tr_Reg, param_file, g_tracer_name) + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,& + obc_src_file_name,obc_src_field_name ) + if(obc_has) then + call set_obgc_segments_props(g_tracer_name,obc_src_file_name,obc_src_field_name) + call register_obgc_segments(GV, CS%OBC, tr_Reg, param_file, g_tracer_name) + endif else call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & restart_CS, longname=longname, units=units) @@ -254,7 +257,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, !! ALE sponges. character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer' - logical :: OK + logical :: OK,obc_has integer :: i, j, k, isc, iec, jsc, jec, nk type(g_tracer_type), pointer :: g_tracer,g_tracer_next character(len=fm_string_len) :: g_tracer_name @@ -356,7 +359,8 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, endif endif - call fill_obgc_segments(G, CS%OBC, tr_ptr, g_tracer_name) + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ) + if(obc_has) call fill_obgc_segments(G, CS%OBC, tr_ptr, g_tracer_name) !traverse the linked list till hit NULL call g_tracer_get_next(g_tracer, g_tracer_next) if (.NOT. associated(g_tracer_next)) exit From cbf2d52699332dafe1fd27a8e3843d0773a42ea3 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Thu, 28 Jan 2021 16:27:57 -0500 Subject: [PATCH 04/21] Increase length of segment data name --- src/core/MOM_open_boundary.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c41068e1fa..6f5b80bc6c 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -82,7 +82,7 @@ module MOM_open_boundary type, public :: OBC_segment_data_type integer :: fid !< handle from FMS associated with segment data on disk integer :: fid_dz !< handle from FMS associated with segment thicknesses on disk - character(len=8) :: name !< a name identifier for the segment data + character(len=32) :: name !< a name identifier for the segment data character(len=8) :: genre !< a family identifier for the segment data real, dimension(:,:,:), allocatable :: buffer_src !< buffer for segment data located at cell faces !! and on the original vertical grid From 038781176c458d4349ef63365bd7a863a5b86d71 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Thu, 28 Jan 2021 16:30:50 -0500 Subject: [PATCH 05/21] Check if tracer is prog --- src/tracer/MOM_generic_tracer.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 0de68fa6ed..466890349b 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -360,7 +360,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, endif call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ) - if(obc_has) call fill_obgc_segments(G, CS%OBC, tr_ptr, g_tracer_name) + if(obc_has .and. g_tracer_is_prog(g_tracer)) call fill_obgc_segments(G, CS%OBC, tr_ptr, g_tracer_name) !traverse the linked list till hit NULL call g_tracer_get_next(g_tracer, g_tracer_next) if (.NOT. associated(g_tracer_next)) exit From 8080fb2380edc3c4580d894ed0b31837ec2d6fa9 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Thu, 4 Mar 2021 12:03:46 -0500 Subject: [PATCH 06/21] Fix restart bug due to setting up T&S reservoirs twice --- src/core/MOM_open_boundary.F90 | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 6f5b80bc6c..1390d587c3 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1911,19 +1911,27 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) end subroutine open_boundary_impose_land_mask !> Make sure the OBC tracer reservoirs are initialized. -subroutine setup_OBC_tracer_reservoirs(G, OBC) +subroutine setup_OBC_tracer_reservoirs(G, OBC, skip_ts) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + logical, optional :: skip_ts !< If true, skip setting up first 2 tracers (T and S) (default: False) ! Local variables type(OBC_segment_type), pointer :: segment => NULL() - integer :: i, j, k, m, n + integer :: i, j, k, m, m0, n + + ! By default, start tracer loop at 1. + ! But, if skip_ts is true, skip T and S and start loop at 3. + m0 = 1 + if(present(skip_ts)) then + if(skip_ts) m0 = 3 + endif do n=1,OBC%number_of_segments segment=>OBC%segment(n) if (associated(segment%tr_Reg)) then if (segment%is_E_or_W) then I = segment%HI%IsdB - do m=1,OBC%ntr + do m=m0,OBC%ntr if (associated(segment%tr_Reg%Tr(m)%tres)) then do k=1,G%ke do j=segment%HI%jsd,segment%HI%jed @@ -1934,7 +1942,7 @@ subroutine setup_OBC_tracer_reservoirs(G, OBC) enddo else J = segment%HI%JsdB - do m=1,OBC%ntr + do m=m0,OBC%ntr if (associated(segment%tr_Reg%Tr(m)%tres)) then do k=1,G%ke do i=segment%HI%isd,segment%HI%ied @@ -4537,7 +4545,7 @@ subroutine fill_obgc_segments(G, OBC, tr_ptr, tr_name) endif segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) enddo - call setup_OBC_tracer_reservoirs(G, OBC) !This will redo the T&S + call setup_OBC_tracer_reservoirs(G, OBC, .true.) ! Skip T and S to avoid restart bug end subroutine fill_obgc_segments subroutine fill_temp_salt_segments(G, OBC, tv) From 104201f7e2034cc8ebd6919c35de56488ddafbaa Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Thu, 4 Mar 2021 15:11:52 -0500 Subject: [PATCH 07/21] Cleanup missed merge --- src/core/MOM_open_boundary.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 41fc3d13dc..6597e29d05 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -3808,11 +3808,9 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Thickness [H ~> m or kg m-2] type(time_type), intent(in) :: Time !< Model time ! Local variables -==== BASE ==== - integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed - integer :: IsdB, IedB, JsdB, JedB, n, m, nz - character(len=40) :: mdl = "set_OBC_segment_data" ! This subroutine's name. -==== BASE ==== + integer :: c, i, j, k, is, ie, js, je, isd, ied, jsd, jed + integer :: IsdB, IedB, JsdB, JedB, n, m, nz, nt, it + character(len=40) :: mdl = "update_OBC_segment_data" ! This subroutine's name. character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path type(OBC_segment_type), pointer :: segment => NULL() integer, dimension(4) :: siz From f57c8c787df9b8d9efb3c37a93f5a3d1ecac9fd6 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Thu, 4 Mar 2021 15:29:57 -0500 Subject: [PATCH 08/21] Fix merge of MOM_generic_tracer.F90 --- src/tracer/MOM_generic_tracer.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 96a7cb2255..644fff9efd 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -59,7 +59,6 @@ module MOM_generic_tracer !> An state hidden in module data that is very much not allowed in MOM6 ! ### This needs to be fixed logical :: g_registered = .false. - integer, parameter :: fm_string_len=128 !>string lengths used in obgc packages public register_MOM_generic_tracer, initialize_MOM_generic_tracer public MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state @@ -361,7 +360,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, endif call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ) - if(obc_has .and. g_tracer_is_prog(g_tracer)) call fill_obgc_segments(G, CS%OBC, tr_ptr, g_tracer_name) + if(obc_has .and. g_tracer_is_prog(g_tracer)) call fill_obgc_segments(G, GV, CS%OBC, tr_ptr, g_tracer_name) !traverse the linked list till hit NULL call g_tracer_get_next(g_tracer, g_tracer_next) if (.NOT. associated(g_tracer_next)) exit From 528ccb6794fc1d882773f86fa140153eabc75587 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Thu, 4 Mar 2021 15:50:38 -0500 Subject: [PATCH 09/21] Add null subroutine g_tracer_get_obc_segment_props --- .../external/GFDL_ocean_BGC/generic_tracer_utils.F90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 b/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 index de513a7f11..f685a5f114 100644 --- a/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 +++ b/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 @@ -25,6 +25,8 @@ module g_tracer_utils character(len=fm_string_len) :: src_var_name !< Tracer source variable name character(len=fm_string_len) :: src_var_unit !< Tracer source variable units character(len=fm_string_len) :: src_var_gridspec !< Tracer source grid file name + character(len=fm_string_len) :: obc_src_file_name !< Boundary condition tracer source filename + character(len=fm_string_len) :: obc_src_field_name !< Boundary condition tracer source fieldname integer :: src_var_record !< Unknown logical :: requires_src_info = .false. !< Unknown real :: src_var_unit_conversion = 1.0 !< This factor depends on the tracer. Ask Jasmin @@ -61,6 +63,7 @@ module g_tracer_utils public :: g_tracer_get_next public :: g_tracer_is_prog public :: g_diag_type + public :: g_tracer_get_obc_segment_props !> Set the values of various (array) members of the tracer node g_tracer_type !! @@ -286,6 +289,14 @@ subroutine g_tracer_get_next(g_tracer,g_tracer_next) type(g_tracer_type), pointer :: g_tracer_next !< Pointer to the next tracer node in the list end subroutine g_tracer_get_next + subroutine g_tracer_get_obc_segment_props(g_tracer_list, name, obc_has, src_file, src_var_name) + type(g_tracer_type), pointer :: g_tracer_list !< pointer to the head of the generic tracer list + type(g_tracer_type), pointer :: g_tracer !< Pointer to tracer node + character(len=*), intent(in) :: name + logical, intent(out):: obc_has !<.true. if This tracer has OBC + character(len=*),optional,intent(out):: src_file, src_var_name !Vertical Diffusion of a tracer node !! !! This subroutine solves a tridiagonal equation to find and set values of vertically diffused field From dfcd0d4a33ebc6790b856458cc8dd440a71ca123 Mon Sep 17 00:00:00 2001 From: andrew-c-ross Date: Wed, 4 Aug 2021 09:47:07 -0400 Subject: [PATCH 10/21] Fix for array already allocated --- src/core/MOM_open_boundary.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index e9a1777061..b90eeaadc6 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4714,9 +4714,6 @@ subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file) segment=>OBC%segment(n) if (.not. segment%on_pe) cycle - if (associated(segment%tr_Reg)) & - call MOM_error(FATAL,"register_temp_salt_segments: tracer array was previously allocated") - name = 'temp' call tracer_name_lookup(tr_Reg, tr_ptr, name) call register_segment_tracer(tr_ptr, param_file, GV, segment, & From b33c000052bd2a2e95ef763c5c08d479c5381fd1 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 4 Aug 2021 10:10:37 -0400 Subject: [PATCH 11/21] Allow OBC segment data to update with a different period for OBGC tracers - This introduces a time scale for updating OBC segment data for OBGC tracers by a parameter override, e.g., 1 day instead of DT via #override UPDATE_OBC_PERIOD_MAX = 86400 - Won't affect T&S or other fields , but it's easy to include any field by adding their names to the if() cycle statement --- src/core/MOM.F90 | 38 ++++++++++++++++++++++++++++++++++ src/core/MOM_open_boundary.F90 | 6 ++++++ 2 files changed, 44 insertions(+) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a57cb01f42..2a7f12edcb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -262,6 +262,10 @@ module MOM !! calculated, and if it is 0, dtbt is calculated every step. type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period. type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated. + real :: update_OBC_period!< The maximum allowed time interval between OBC updates + type(time_type) :: update_OBC_interval !< A time_time representation of update_OBC_period. + type(time_type) :: update_OBC_time !< The next time OBC is applied. + real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied !! by ice shelf [nondim] real, dimension(:,:,:), pointer :: & @@ -1054,6 +1058,17 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & CS%dtbt_reset_time = CS%dtbt_reset_time + CS%dtbt_reset_interval endif endif + !OBC segment data update for some fields can be less frequent than others + if(associated(CS%OBC)) then + CS%OBC%update_OBC_seg_data = .false. + if (CS%update_OBC_period == 0.0) CS%OBC%update_OBC_seg_data = .true. + if (CS%update_OBC_period > 0.0) then + if (Time_local >= CS%update_OBC_time) then !### Change >= to > here. + CS%OBC%update_OBC_seg_data = .true. + CS%update_OBC_time = CS%update_OBC_time + CS%update_OBC_interval + endif + endif + endif call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & @@ -1924,6 +1939,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "DTBT will be set every dynamics time step. The default "//& "is set by DT_THERM. This is only used if SPLIT is true.", & units="s", default=default_val, do_not_read=(dtbt > 0.0)) + + CS%update_OBC_period = -1.0 + call get_param(param_file, "MOM", "UPDATE_OBC_PERIOD_MAX", CS%update_OBC_period, & + "The maximum period between OBC segment data updates. "//& + "If UPDATE_OBC_PERIOD_MAX is negative, DTBT is set based "//& + "only on information available at initialization. If 0, "//& + "OBC updates every dynamics time step. The default "//& + "is set by DT_THERM. This is only used if SPLIT is true.", & + units="s", default=default_val, do_not_read=(dtbt > 0.0)) + endif ! This is here in case these values are used inappropriately. @@ -2622,6 +2647,19 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%dtbt_reset_time = CS%dtbt_reset_time - CS%dtbt_reset_interval endif endif + !Set OBC segment data update period + if (associated(CS%OBC) .and. CS%update_OBC_period > 0.0) then + CS%update_OBC_interval = real_to_time(CS%update_OBC_period) + ! Set update_OBC_time to be the next even multiple of update_OBC_interval. + CS%update_OBC_time = Time_init + CS%update_OBC_interval * & + ((Time - Time_init) / CS%update_OBC_interval) + if ((CS%update_OBC_time > Time) .and. CS%OBC%update_OBC) then + ! Back up dtbt_reset_time one interval to force dtbt to be calculated, + ! because the restart was not aligned with the interval to recalculate + ! dtbt, and dtbt was not read from a restart file. + CS%update_OBC_time = CS%update_OBC_time - CS%update_OBC_interval + endif + endif elseif (CS%use_RK2) then call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, Time, G, GV, US, & param_file, diag, CS%dyn_unsplit_RK2_CSp, restart_CSp, & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 6597e29d05..62c979849b 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -241,6 +241,8 @@ module MOM_open_boundary logical :: user_BCs_set_globally = .false. !< True if any OBC_USER_CONFIG is set !! for input from user directory. logical :: update_OBC = .false. !< Is OBC data time-dependent + logical :: update_OBC_seg_data = .false. !< Is it the time for OBC segment data update for fields that + !! require less frequent update logical :: needs_IO_for_data = .false. !< Is any i/o needed for OBCs logical :: zero_vorticity = .false. !< If True, sets relative vorticity to zero on open boundaries. logical :: freeslip_vorticity = .false. !< If True, sets normal gradient of tangential velocity to zero @@ -3897,6 +3899,10 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) allocate(h_stack(GV%ke)) h_stack(:) = 0.0 do m = 1,segment%num_fields + !This field may not require a high frequency OBC segment update and might be allowed + !a less frequent update as set by the parameter update_OBC_period_max in MOM.F90. + !Cycle if it is not the time to update OBC segment data for this field. + if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle if (segment%field(m)%fid > 0) then siz(1)=size(segment%field(m)%buffer_src,1) siz(2)=size(segment%field(m)%buffer_src,2) From e4bd8ebdf1ff5f1d963fd18384aafe7149a7fe15 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 4 Aug 2021 13:24:12 -0400 Subject: [PATCH 12/21] Allow OBC segment data to update with a different period for OBGC tracers - This introduces a time scale for updating OBC segment data for OBGC tracers by a parameter override, e.g., 1 day instead of DT via #override UPDATE_OBC_PERIOD_MAX = 86400 - Won't affect T&S or other fields , but it's easy to include any field by adding their names to the if() cycle statement --- src/core/MOM.F90 | 41 +++++++++++++++++++---------------------- 1 file changed, 19 insertions(+), 22 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2a7f12edcb..d73108b740 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -262,9 +262,9 @@ module MOM !! calculated, and if it is 0, dtbt is calculated every step. type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period. type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated. - real :: update_OBC_period!< The maximum allowed time interval between OBC updates - type(time_type) :: update_OBC_interval !< A time_time representation of update_OBC_period. - type(time_type) :: update_OBC_time !< The next time OBC is applied. + real :: dt_obc_seg_period!< The time interval between OBC segment updates for OBGC tracers + type(time_type) :: dt_obc_seg_interval !< A time_time representation of dt_obc_seg_period. + type(time_type) :: dt_obc_seg_time !< The next time OBC segment update is applied to OBGC tracers. real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied !! by ice shelf [nondim] @@ -1061,11 +1061,11 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & !OBC segment data update for some fields can be less frequent than others if(associated(CS%OBC)) then CS%OBC%update_OBC_seg_data = .false. - if (CS%update_OBC_period == 0.0) CS%OBC%update_OBC_seg_data = .true. - if (CS%update_OBC_period > 0.0) then - if (Time_local >= CS%update_OBC_time) then !### Change >= to > here. + if (CS%dt_obc_seg_period == 0.0) CS%OBC%update_OBC_seg_data = .true. + if (CS%dt_obc_seg_period > 0.0) then + if (Time_local >= CS%dt_obc_seg_time) then !### Change >= to > here. CS%OBC%update_OBC_seg_data = .true. - CS%update_OBC_time = CS%update_OBC_time + CS%update_OBC_interval + CS%dt_obc_seg_time = CS%dt_obc_seg_time + CS%dt_obc_seg_interval endif endif endif @@ -1940,14 +1940,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "is set by DT_THERM. This is only used if SPLIT is true.", & units="s", default=default_val, do_not_read=(dtbt > 0.0)) - CS%update_OBC_period = -1.0 - call get_param(param_file, "MOM", "UPDATE_OBC_PERIOD_MAX", CS%update_OBC_period, & - "The maximum period between OBC segment data updates. "//& - "If UPDATE_OBC_PERIOD_MAX is negative, DTBT is set based "//& - "only on information available at initialization. If 0, "//& - "OBC updates every dynamics time step. The default "//& - "is set by DT_THERM. This is only used if SPLIT is true.", & - units="s", default=default_val, do_not_read=(dtbt > 0.0)) + CS%dt_obc_seg_period = -1.0 + call get_param(param_file, "MOM", "DT_OBC_SEG_UPDATE_OBGC", CS%dt_obc_seg_period, & + "The time between OBC segment data updates for OBGC tracers. "//& + "The default is set to DT. This is only used if SPLIT is true.", & + units="s", default=CS%dt) endif @@ -2648,16 +2645,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif endif !Set OBC segment data update period - if (associated(CS%OBC) .and. CS%update_OBC_period > 0.0) then - CS%update_OBC_interval = real_to_time(CS%update_OBC_period) - ! Set update_OBC_time to be the next even multiple of update_OBC_interval. - CS%update_OBC_time = Time_init + CS%update_OBC_interval * & - ((Time - Time_init) / CS%update_OBC_interval) - if ((CS%update_OBC_time > Time) .and. CS%OBC%update_OBC) then + if (associated(CS%OBC) .and. CS%dt_obc_seg_period > 0.0) then + CS%dt_obc_seg_interval = real_to_time(CS%dt_obc_seg_period) + ! Set dt_obc_seg_time to be the next even multiple of dt_obc_seg_interval. + CS%dt_obc_seg_time = Time_init + CS%dt_obc_seg_interval * & + ((Time - Time_init) / CS%dt_obc_seg_interval) + if ((CS%dt_obc_seg_time > Time) .and. CS%OBC%update_OBC_seg_data) then ! Back up dtbt_reset_time one interval to force dtbt to be calculated, ! because the restart was not aligned with the interval to recalculate ! dtbt, and dtbt was not read from a restart file. - CS%update_OBC_time = CS%update_OBC_time - CS%update_OBC_interval + CS%dt_obc_seg_time = CS%dt_obc_seg_time - CS%dt_obc_seg_interval endif endif elseif (CS%use_RK2) then From 201288ecdc2a070ca6e971e1ef4f546fdd7a619b Mon Sep 17 00:00:00 2001 From: andrew-c-ross Date: Wed, 4 Aug 2021 16:47:17 -0400 Subject: [PATCH 13/21] Place call call_tracer_register in the correct order to avoid having extra tracers registered first on boundary --- src/core/MOM.F90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index b965c431c6..6f27e577ef 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2303,11 +2303,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%dyn_unsplit_CSp, restart_CSp) endif - ! This subroutine calls user-specified tracer registration routines. - ! Additional calls can be added to MOM_tracer_flow_control.F90. - call call_tracer_register(dG%HI, GV, US, param_file, CS%tracer_flow_CSp, & - CS%tracer_Reg, restart_CSp, CS%OBC) - call MEKE_alloc_register_restart(dG%HI, param_file, CS%MEKE, restart_CSp) call set_visc_register_restarts(dG%HI, GV, param_file, CS%visc, restart_CSp) call mixedlayer_restrat_register_restarts(dG%HI, param_file, & @@ -2336,7 +2331,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! could occur with the call to update_OBC_data or after the main initialization. if (use_temperature) & call register_temp_salt_segments(GV, CS%OBC, CS%tracer_Reg, param_file) + endif + ! This subroutine calls user-specified tracer registration routines. + ! Additional calls can be added to MOM_tracer_flow_control.F90. + ! Needs to be after registering temperature and salinity OBCs above, + ! or else the user-specified tracers will be first. + call call_tracer_register(dG%HI, GV, US, param_file, CS%tracer_flow_CSp, & + CS%tracer_Reg, restart_CSp, CS%OBC) + + if (associated(CS%OBC)) then ! This needs the number of tracers and to have called any code that sets whether ! reservoirs are used. call open_boundary_register_restarts(dg%HI, GV, CS%OBC, CS%tracer_Reg, & From d9fffc2609781fa35bd5e21b6783ee85bf6a1f90 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Mon, 27 Dec 2021 11:27:42 -0500 Subject: [PATCH 14/21] Attempt to fix errors from merge --- src/core/MOM_open_boundary.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 7990c8c084..02d55ae86e 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -80,8 +80,8 @@ module MOM_open_boundary type, public :: OBC_segment_data_type integer :: fid !< handle from FMS associated with segment data on disk integer :: fid_dz !< handle from FMS associated with segment thicknesses on disk - character(len=32) :: name !< a name identifier for the segment data - character(len=8) :: name !< a name identifier for the segment data + character(len=32) :: name !< a name identifier for the segment data + character(len=8) :: genre !< a family identifier for the segment data real, allocatable :: buffer_src(:,:,:) !< buffer for segment data located at cell faces !! and on the original vertical grid integer :: nk_src !< Number of vertical levels in the source data @@ -3543,7 +3543,7 @@ function get_tracer_index(OBC_seg,tr_name) integer :: get_tracer_index, it get_tracer_index=-1 it=1 - do while(associated(OBC_seg%tr_Reg%Tr(it)%t)) + do while(allocated(OBC_seg%tr_Reg%Tr(it)%t)) if (trim(OBC_seg%tr_Reg%Tr(it)%name) == trim(tr_name)) then get_tracer_index=it exit @@ -4405,7 +4405,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if(nt .lt. 0) then call MOM_error(FATAL,"update_OBC_segment_data: Did not find tracer "//trim(segment%field(m)%name)) endif - if (associated(segment%field(m)%buffer_dst)) then + if (allocatedZ(segment%field(m)%buffer_dst)) then do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc segment%tr_Reg%Tr(nt)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k) enddo ; enddo ; enddo From 1ec41d51b8e2ebe531f455e5363d08babf8f80c1 Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Mon, 27 Dec 2021 11:40:25 -0500 Subject: [PATCH 15/21] Stray typo --- src/core/MOM_open_boundary.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 02d55ae86e..f7feda4aaf 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4405,7 +4405,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if(nt .lt. 0) then call MOM_error(FATAL,"update_OBC_segment_data: Did not find tracer "//trim(segment%field(m)%name)) endif - if (allocatedZ(segment%field(m)%buffer_dst)) then + if (allocated(segment%field(m)%buffer_dst)) then do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc segment%tr_Reg%Tr(nt)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k) enddo ; enddo ; enddo From 4ca2bcfdd60dbd396f6e16d979bcbad8909c714d Mon Sep 17 00:00:00 2001 From: Andrew C Ross Date: Mon, 27 Dec 2021 11:57:19 -0500 Subject: [PATCH 16/21] Remove duplicated tracer register from merge --- src/core/MOM.F90 | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c87f57253e..aca2da3cec 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2422,11 +2422,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%dyn_unsplit_CSp) endif - ! This subroutine calls user-specified tracer registration routines. - ! Additional calls can be added to MOM_tracer_flow_control.F90. - call call_tracer_register(HI, GV, US, param_file, CS%tracer_flow_CSp, & - CS%tracer_Reg, restart_CSp) - call MEKE_alloc_register_restart(HI, param_file, CS%MEKE, restart_CSp) call set_visc_register_restarts(HI, GV, param_file, CS%visc, restart_CSp) call mixedlayer_restrat_register_restarts(HI, param_file, & @@ -2461,7 +2456,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Additional calls can be added to MOM_tracer_flow_control.F90. ! Needs to be after registering temperature and salinity OBCs above, ! or else the user-specified tracers will be first. - call call_tracer_register(dG%HI, GV, US, param_file, CS%tracer_flow_CSp, & + call call_tracer_register(HI, GV, US, param_file, CS%tracer_flow_CSp, & CS%tracer_Reg, restart_CSp, CS%OBC) if (associated(CS%OBC)) then From 11bfe1401b86b8a7f0e5d58ef33d5345b585bcd1 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 7 Jan 2022 16:26:51 -0500 Subject: [PATCH 17/21] Allow obgc tracers have different reservoir length scales - This update along with an update to generic_tracers_utils.F90 brings in the capability to set tracer reservoir (inverse) length scales for individual obgc tracers via ocean_bgc field_table mechanism. - Tracer resevoir inverse length scales are set as before to 1.0/Lscale_in and 1.0/Lscale_out - Users should be able to specify a multiplicative factor for each obgc tracer for each IN or OUT direction in the field_table, e.g., alk_obc_lfac_in = 0.5 alk_obc_lfac_out = 0.75 - The factors are multiplied by the inverse length scales before reservoir calculations. - The default factors for all tracers are 1. - No answer changes unless non-one factors are specified in field_table --- src/core/MOM_open_boundary.F90 | 60 ++++++++++++++++++------------- src/tracer/MOM_generic_tracer.F90 | 5 +-- 2 files changed, 38 insertions(+), 27 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f7feda4aaf..5dd23e6b0f 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -88,7 +88,11 @@ module MOM_open_boundary real, allocatable :: dz_src(:,:,:) !< vertical grid cell spacing of the incoming segment !! data, set in [Z ~> m] then scaled to [H ~> m or kg m-2] real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid - real :: value !< constant value if fid is equal to -1 + real :: value !< constant value if fid is equal to -1 + real :: resrv_lfac_in = 1. !< reservoir inverse length scale factor for IN direction per field + !< the general 1/Lscale_IN is multiplied by this factor for each tracer + real :: resrv_lfac_out= 1. !< reservoir inverse length scale factor for OUT direction per field + !< the general 1/Lscale_OUT is multiplied by this factor for each tracer end type OBC_segment_data_type !> Tracer on OBC segment data structure, for putting into a segment tracer registry. @@ -206,13 +210,13 @@ module MOM_open_boundary !! can occur [T-1 ~> s-1]. type(segment_tracer_registry_type), pointer :: tr_Reg=> NULL()!< A pointer to the tracer registry for the segment. type(hor_index_type) :: HI !< Horizontal index ranges - real :: Tr_InvLscale_out !< An effective inverse length scale for restoring - !! the tracer concentration in a fictitious - !! reservoir towards interior values when flow - !! is exiting the domain [L-1 ~> m-1] - real :: Tr_InvLscale_in !< An effective inverse length scale for restoring - !! the tracer concentration towards an externally - !! imposed value when flow is entering [L-1 ~> m-1] + real, allocatable :: Tr_InvLscale_out !< An effective inverse length scale for restoring + !! the tracer concentration in a fictitious + !! reservoir towards interior values when flow + !! is exiting the domain [L-1 ~> m-1] + real, allocatable :: Tr_InvLscale_in !< An effective inverse length scale for restoring + !! the tracer concentration towards an externally + !! imposed value when flow is entering [L-1 ~> m-1] end type OBC_segment_type !> Open-boundary data @@ -343,6 +347,7 @@ module MOM_open_boundary character(len=128) :: tracer_name character(len=128) :: tracer_src_file character(len=128) :: tracer_src_field + real :: lfac_in,lfac_out end type external_tracers_segments_props type(external_tracers_segments_props), pointer, save :: obgc_segments_props => NULL() !< Linked-list of obgc tracers properties integer, save :: num_obgc_tracers = 0 !< Keeps the total number of obgc tracers @@ -753,7 +758,8 @@ subroutine initialize_segment_data(G, OBC, PF) value, filename, fieldname) else segment%field(m)%genre='obgc' - call get_obgc_segments_props(obgc_segments_props_list,fields(m),filename,fieldname) + call get_obgc_segments_props(obgc_segments_props_list,fields(m),filename,fieldname,& + segment%field(m)%resrv_lfac_in,segment%field(m)%resrv_lfac_out) do mm=1,num_fields if(trim(fields(m))==trim(fields(mm))) then if (is_root_pe()) & @@ -4698,32 +4704,36 @@ subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file) end subroutine register_temp_salt_segments !> Sets the OBC properties of external obgc tracers, such as their source file and field name -subroutine set_obgc_segments_props(tr_name,obc_src_file_name,obc_src_field_name) - character(len=*), intent(in) :: tr_name !< Tracer name - character(len=*), intent(in) :: obc_src_file_name !< OBC source file name - character(len=*), intent(in) :: obc_src_field_name !< name of the field in the source file - +subroutine set_obgc_segments_props(tr_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + character(len=*), intent(in) :: tr_name !< Tracer name + character(len=*), intent(in) :: obc_src_file_name !< OBC source file name + character(len=*), intent(in) :: obc_src_field_name !< name of the field in the source file + real, intent(in) :: lfac_in,lfac_out type(external_tracers_segments_props),pointer :: node_ptr => NULL() allocate(node_ptr) node_ptr%tracer_name = trim(tr_name) node_ptr%tracer_src_file = trim(obc_src_file_name) node_ptr%tracer_src_field = trim(obc_src_field_name) + node_ptr%lfac_in = lfac_in + node_ptr%lfac_out = lfac_out !Reversed Linked List implementation! Make this new node to be the head of the list. node_ptr%next => obgc_segments_props obgc_segments_props => node_ptr num_obgc_tracers = num_obgc_tracers+1 end subroutine set_obgc_segments_props -!> Get the OBC properties of external obgc tracers, such as their source file and field name -subroutine get_obgc_segments_props(node, tr_name,obc_src_file_name,obc_src_field_name) +!> Get the OBC properties of external obgc tracers, such as their source file, field name, reservoir length scale factors +subroutine get_obgc_segments_props(node, tr_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) type(external_tracers_segments_props),pointer :: node - character(len=*), intent(out) :: tr_name !< Tracer name - character(len=*), intent(out) :: obc_src_file_name !< OBC source file name - character(len=*), intent(out) :: obc_src_field_name !< name of the field in the source file - + character(len=*), intent(out) :: tr_name !< Tracer name + character(len=*), intent(out) :: obc_src_file_name !< OBC source file name + character(len=*), intent(out) :: obc_src_field_name !< name of the field in the source file + real, intent(out) :: lfac_in,lfac_out !< multiplicative factors for inverse reservoir length scale tr_name=trim(node%tracer_name) obc_src_file_name=trim(node%tracer_src_file) obc_src_field_name=trim(node%tracer_src_field) + lfac_in =node%lfac_in + lfac_out=node%lfac_out !pop the head. node => node%next @@ -5243,9 +5253,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(I+ishift,j) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz - u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out / & + u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / & ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j))) - u_L_in = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in / & + u_L_in = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / & ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j))) fac1 = 1.0 + (u_L_out-u_L_in) segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + & @@ -5268,9 +5278,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(i,j+jshift) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz - v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out / & + v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / & ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J))) - v_L_in = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in / & + v_L_in = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / & ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J))) fac1 = 1.0 + (v_L_out-v_L_in) segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + & @@ -5597,7 +5607,7 @@ subroutine rotate_OBC_segment_config(segment_in, G_in, segment, G, turns) end select ! These are conditionally set if Lscale_{in,out} are present - segment%Tr_InvLscale_in = segment_in%Tr_InvLscale_in + segment%Tr_InvLscale_in = segment_in%Tr_InvLscale_in segment%Tr_InvLscale_out = segment_in%Tr_InvLscale_out end subroutine rotate_OBC_segment_config diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 684f38ee19..d87b43f7cd 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -118,6 +118,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, type(g_tracer_type), pointer :: g_tracer,g_tracer_next character(len=fm_string_len) :: g_tracer_name,longname,units character(len=fm_string_len) :: obc_src_file_name,obc_src_field_name + real :: lfac_in,lfac_out real, dimension(:,:,:,:), pointer :: tr_field real, dimension(:,:,:), pointer :: tr_ptr real, dimension(HI%isd:HI%ied, HI%jsd:HI%jed,GV%ke) :: grid_tmask @@ -210,9 +211,9 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS, restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) if (associated(CS%OBC)) & call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,& - obc_src_file_name,obc_src_field_name ) + obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) if(obc_has) then - call set_obgc_segments_props(g_tracer_name,obc_src_file_name,obc_src_field_name) + call set_obgc_segments_props(g_tracer_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) call register_obgc_segments(GV, CS%OBC, tr_Reg, param_file, g_tracer_name) endif else From 7bac290168df138242e84302d46ceae6c87b9cff Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Mon, 10 Jan 2022 12:16:06 -0500 Subject: [PATCH 18/21] Allow obgc tracers have different reservoir length scales - This update along with an update to generic_tracers_utils.F90 brings in the capability to set tracer reservoir (inverse) length scales for individual obgc tracers via ocean_bgc field_table mechanism. - Tracer resevoir inverse length scales are set as before to 1.0/Lscale_in and 1.0/Lscale_out - Users should be able to specify a multiplicative factor for each obgc tracer for each IN or OUT direction in the field_table, e.g., alk_obc_lfac_in = 0.5 alk_obc_lfac_out = 0.75 - The factors are multiplied by the inverse length scales before reservoir calculations. - The default factors for all tracers are 1. - No answer changes unless non-one factors are specified in field_table --- src/core/MOM_open_boundary.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 5dd23e6b0f..7f9f1962f7 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -210,13 +210,13 @@ module MOM_open_boundary !! can occur [T-1 ~> s-1]. type(segment_tracer_registry_type), pointer :: tr_Reg=> NULL()!< A pointer to the tracer registry for the segment. type(hor_index_type) :: HI !< Horizontal index ranges - real, allocatable :: Tr_InvLscale_out !< An effective inverse length scale for restoring - !! the tracer concentration in a fictitious - !! reservoir towards interior values when flow - !! is exiting the domain [L-1 ~> m-1] - real, allocatable :: Tr_InvLscale_in !< An effective inverse length scale for restoring - !! the tracer concentration towards an externally - !! imposed value when flow is entering [L-1 ~> m-1] + real :: Tr_InvLscale_out !< An effective inverse length scale for restoring + !! the tracer concentration in a fictitious + !! reservoir towards interior values when flow + !! is exiting the domain [L-1 ~> m-1] + real :: Tr_InvLscale_in !< An effective inverse length scale for restoring + !! the tracer concentration towards an externally + !! imposed value when flow is entering [L-1 ~> m-1] end type OBC_segment_type !> Open-boundary data From 2dcc7613301f3bed41127e81de619b1b93ac2404 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 23 Dec 2021 10:11:17 -0500 Subject: [PATCH 19/21] Fix restart issue for OBC, reservoir length scale zero - This update fixes the restart issues for OBC runs if reservoir length scale is set to zero - There is still a restart issue when the reservoir length scale is non-zero even for the physical (noCOBALT) models that has to be fixed. --- src/core/MOM_open_boundary.F90 | 5 +++-- src/initialization/MOM_state_initialization.F90 | 7 ++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 7f9f1962f7..48044a2223 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -64,6 +64,7 @@ module MOM_open_boundary public rotate_OBC_config public rotate_OBC_init public initialize_segment_data +public setup_OBC_tracer_reservoirs integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary integer, parameter, public :: OBC_SIMPLE = 1 !< Indicates the use of a simple inflow open boundary @@ -4809,7 +4810,7 @@ subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) endif segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) enddo - call setup_OBC_tracer_reservoirs(GV, OBC, .true.) ! Skip T and S to avoid restart bug +! call setup_OBC_tracer_reservoirs(GV, OBC, .true.) ! Skip T and S to avoid restart bug end subroutine fill_obgc_segments subroutine fill_temp_salt_segments(G, GV, OBC, tv) @@ -4868,7 +4869,7 @@ subroutine fill_temp_salt_segments(G, GV, OBC, tv) segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) enddo - call setup_OBC_tracer_reservoirs(GV, OBC) +! call setup_OBC_tracer_reservoirs(GV, OBC) end subroutine fill_temp_salt_segments !> Find the region outside of all open boundary segments and diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 8055440cce..86963d2968 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -24,7 +24,7 @@ module MOM_state_initialization use MOM_open_boundary, only : open_boundary_query use MOM_open_boundary, only : set_tracer_data, initialize_segment_data use MOM_open_boundary, only : open_boundary_test_extern_h -use MOM_open_boundary, only : fill_temp_salt_segments +use MOM_open_boundary, only : fill_temp_salt_segments,setup_OBC_tracer_reservoirs use MOM_open_boundary, only : update_OBC_segment_data !use MOM_open_boundary, only : set_3D_OBC_data use MOM_grid_initialize, only : initialize_masks, set_grid_metrics @@ -401,9 +401,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & end select endif endif ! not from_Z_file. - if (use_temperature .and. use_OBC) & + if (use_temperature .and. use_OBC) then call fill_temp_salt_segments(G, GV, OBC, tv) - + call setup_OBC_tracer_reservoirs(GV, OBC) + endif ! The thicknesses in halo points might be needed to initialize the velocities. if (new_sim) call pass_var(h, G%Domain) From 159ba650cf4fe6d0e40f5f20aa8b6c56bf570ebd Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 6 Apr 2022 19:15:54 -0400 Subject: [PATCH 20/21] Fix restart issue step 1 - Apply halo updates to h,uhr,vhr before use in tracer reservoir update fixes the restart issue for the physical model T&S fixes the restart issue for the BGC tracers when segment update period = DT - There is still a restart issue for BGC tracers when segment update period > DT - However the later restart issue is absent when tracer reservoir INPUT length scale is 0 which points to an issue in segment%tr_Reg%Tr(m)%t --- src/core/MOM_open_boundary.F90 | 52 +++++++++++++++++++++++----------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 48044a2223..7ea9a181ba 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -5234,12 +5234,31 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) real :: fac1 ! The denominator of the expression for tracer updates [nondim] integer :: i, j, k, m, n, ntr, nz integer :: ishift, idir, jshift, jdir + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h1 + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhr1 + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhr1 + if(.NOT. associated(OBC)) return nz = GV%ke ntr = Reg%ntr - if (associated(OBC)) then ; if (OBC%OBC_pe) then ; do n=1,OBC%number_of_segments - segment=>OBC%segment(n) - if (.not. associated(segment%tr_Reg)) cycle + !Input arrays h,uhr,vhr have not had halo updates since they were last updated outside this routine. + !This was causing restart issues for experiments that include tracers with reservoirs. + !We need to perform halo updates for input arrays, since they are intent(in) they cannot be modified. + !Hence we make copies of them. + h1=h + uhr1=uhr + vhr1=vhr + call pass_var(h1, G%Domain) + call pass_vector(uhr1,vhr1, G%Domain) + + if (OBC%OBC_pe) then ; do n=1,OBC%number_of_segments + segment=>OBC%segment(n) + if (.not. associated(segment%tr_Reg)) cycle + do m=1,ntr + if (.not. allocated(segment%tr_Reg%Tr(m)%tres)) cycle + !OBGC tracers may have less frequent segment updates + if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle + if (segment%is_E_or_W) then I = segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed @@ -5253,17 +5272,17 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ! Can keep this or take it out, either way if (G%mask2dT(I+ishift,j) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep - do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz - u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / & - ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j))) - u_L_in = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / & - ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j))) + do k=1,nz + u_L_out = max(0.0, (idir*uhr1(I,j,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / & + ((h1(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j))) + u_L_in = min(0.0, (idir*uhr1(I,j,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / & + ((h1(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j))) fac1 = 1.0 + (u_L_out-u_L_in) segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + & (u_L_out*Reg%Tr(m)%t(I+ishift,j,k) - & u_L_in*segment%tr_Reg%Tr(m)%t(I,j,k))) if (allocated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = segment%tr_Reg%Tr(m)%tres(I,j,k) - enddo ; endif ; enddo + enddo enddo elseif (segment%is_N_or_S) then J = segment%HI%JsdB @@ -5278,20 +5297,21 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ! Can keep this or take it out, either way if (G%mask2dT(i,j+jshift) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep - do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz - v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / & - ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J))) - v_L_in = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / & - ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J))) + do k=1,nz + v_L_out = max(0.0, (jdir*vhr1(i,J,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / & + ((h1(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J))) + v_L_in = min(0.0, (jdir*vhr1(i,J,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / & + ((h1(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J))) fac1 = 1.0 + (v_L_out-v_L_in) segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + & (v_L_out*Reg%Tr(m)%t(i,J+jshift,k) - & v_L_in*segment%tr_Reg%Tr(m)%t(i,J,k))) if (allocated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = segment%tr_Reg%Tr(m)%tres(i,J,k) - enddo ; endif ; enddo + enddo enddo endif - enddo ; endif ; endif + enddo !m=1,ntr + enddo ; endif end subroutine update_segment_tracer_reservoirs From 5ee87fc74260f72b7300ea6564f3032a6c36a204 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 12 Apr 2022 13:55:14 -0400 Subject: [PATCH 21/21] Fix restart issue step 2 - Apply halo updates to h,uhr,vhr before use in tracer reservoir update fixes the restart issue for the physical model T&S fixes the restart issue for the BGC tracers when segment update period = DT - Fix the restart issue for BGC tracers when segment update period > DT fixes the logic for setting the update time (it was not correct) Note that DT_OBC_SEG_UPDATE_OBGC has to be an integer multiple of DT and DT_THERM ! --- src/core/MOM.F90 | 55 +++++++++++++++------------------- src/core/MOM_open_boundary.F90 | 4 +-- 2 files changed, 26 insertions(+), 33 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index aca2da3cec..e5fde8082d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1036,7 +1036,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 @@ -1082,6 +1081,17 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call disable_averaging(CS%diag) endif + !OBC segment data update for some fields can be less frequent than others + if(associated(CS%OBC)) then + CS%OBC%update_OBC_seg_data = .false. + if (CS%dt_obc_seg_period == 0.0) CS%OBC%update_OBC_seg_data = .true. + if (CS%dt_obc_seg_period > 0.0) then + if (Time_local >= CS%dt_obc_seg_time) then !### Change >= to > here. + CS%OBC%update_OBC_seg_data = .true. + CS%dt_obc_seg_time = CS%dt_obc_seg_time + CS%dt_obc_seg_interval + endif + endif + endif if (CS%do_dynamics .and. CS%split) then !--------------------------- start SPLIT ! This section uses a split time stepping scheme for the dynamic equations, @@ -1095,17 +1105,6 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & CS%dtbt_reset_time = CS%dtbt_reset_time + CS%dtbt_reset_interval endif endif - !OBC segment data update for some fields can be less frequent than others - if(associated(CS%OBC)) then - CS%OBC%update_OBC_seg_data = .false. - if (CS%dt_obc_seg_period == 0.0) CS%OBC%update_OBC_seg_data = .true. - if (CS%dt_obc_seg_period > 0.0) then - if (Time_local >= CS%dt_obc_seg_time) then !### Change >= to > here. - CS%OBC%update_OBC_seg_data = .true. - CS%dt_obc_seg_time = CS%dt_obc_seg_time + CS%dt_obc_seg_interval - endif - endif - endif call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & @@ -2009,14 +2008,15 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "DTBT will be set every dynamics time step. The default "//& "is set by DT_THERM. This is only used if SPLIT is true.", & units="s", default=default_val, do_not_read=(dtbt > 0.0)) + endif - CS%dt_obc_seg_period = -1.0 - call get_param(param_file, "MOM", "DT_OBC_SEG_UPDATE_OBGC", CS%dt_obc_seg_period, & - "The time between OBC segment data updates for OBGC tracers. "//& - "The default is set to DT. This is only used if SPLIT is true.", & - units="s", default=CS%dt) + CS%dt_obc_seg_period = -1.0 + call get_param(param_file, "MOM", "DT_OBC_SEG_UPDATE_OBGC", CS%dt_obc_seg_period, & + "The time between OBC segment data updates for OBGC tracers. "//& + "This must be an integer multiple of DT and DT_THERM. "//& + "The default is set to DT.", & + units="s", default=CS%dt) - endif ! This is here in case these values are used inappropriately. use_frazil = .false. ; bound_salinity = .false. @@ -2740,19 +2740,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%dtbt_reset_time = CS%dtbt_reset_time - CS%dtbt_reset_interval endif endif - !Set OBC segment data update period - if (associated(CS%OBC) .and. CS%dt_obc_seg_period > 0.0) then - CS%dt_obc_seg_interval = real_to_time(CS%dt_obc_seg_period) - ! Set dt_obc_seg_time to be the next even multiple of dt_obc_seg_interval. - CS%dt_obc_seg_time = Time_init + CS%dt_obc_seg_interval * & - ((Time - Time_init) / CS%dt_obc_seg_interval) - if ((CS%dt_obc_seg_time > Time) .and. CS%OBC%update_OBC_seg_data) then - ! Back up dtbt_reset_time one interval to force dtbt to be calculated, - ! because the restart was not aligned with the interval to recalculate - ! dtbt, and dtbt was not read from a restart file. - CS%dt_obc_seg_time = CS%dt_obc_seg_time - CS%dt_obc_seg_interval - endif - endif elseif (CS%use_RK2) then call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, Time, G, GV, US, & param_file, diag, CS%dyn_unsplit_RK2_CSp, & @@ -2767,6 +2754,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%ntrunc, cont_stencil=CS%cont_stencil) endif + !Set OBC segment data update period + if (associated(CS%OBC) .and. CS%dt_obc_seg_period > 0.0) then + CS%dt_obc_seg_interval = real_to_time(CS%dt_obc_seg_period) + CS%dt_obc_seg_time = Time + CS%dt_obc_seg_interval + endif + call callTree_waypoint("dynamics initialized (initialize_MOM)") CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 7ea9a181ba..6399757089 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4220,6 +4220,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! Start second loop to update all fields now that data for all fields are available. ! (split because tides depend on multiple variables). do m = 1,segment%num_fields + !cycle if it is not the time to update OBGC tracers from source + if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle ! if (segment%field(m)%fid>0) then ! calculate external BT velocity and transport if needed if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then @@ -5256,8 +5258,6 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (.not. associated(segment%tr_Reg)) cycle do m=1,ntr if (.not. allocated(segment%tr_Reg%Tr(m)%tres)) cycle - !OBGC tracers may have less frequent segment updates - if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle if (segment%is_E_or_W) then I = segment%HI%IsdB