From 79e579713b9a326f6f5c1594024ee624b530159c Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 4 Jun 2020 13:06:26 -0400 Subject: [PATCH 1/5] 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 2/5] 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 3/5] 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 d06e1444984e6be7bda5c5efe5a00519b0f9185a Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 9 Oct 2020 15:45:41 -0400 Subject: [PATCH 4/5] Fix to fill OBC segments only for prognostic tracers --- 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 8acd7d04915c5211281941d005c808526a921470 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 27 Jul 2021 19:47:37 -0400 Subject: [PATCH 5/5] Fixes and set OBC_update period - These are some fixes towards understanding the restart issue for OBC models. - Also, introduce a time period for OBC_segment update --- src/core/MOM.F90 | 38 ++++++++++- src/core/MOM_open_boundary.F90 | 65 +++++++++++-------- .../MOM_state_initialization.F90 | 7 +- 3 files changed, 77 insertions(+), 33 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 69c759e4e2..95a6c090ce 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -259,6 +259,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 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 :: & @@ -1021,7 +1024,6 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call disable_averaging(CS%diag) endif - if (CS%do_dynamics .and. CS%split) then !--------------------------- start SPLIT ! This section uses a split time stepping scheme for the dynamic equations, ! basically the stacked shallow water equations with viscosity. @@ -1035,6 +1037,17 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif endif + !OBC hack + if(associated(CS%OBC)) then + CS%OBC%update_OBC = .false. + if (CS%update_OBC_period == 0.0) CS%OBC%update_OBC = .true. + if (CS%update_OBC_period > 0.0) then + if (Time_local >= CS%update_OBC_time) then !### Change >= to > here. + CS%OBC%update_OBC = .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, & CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, & @@ -1882,6 +1895,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)) + !OBC hack + CS%update_OBC_period = -1.0 + call get_param(param_file, "MOM", "UPDATE_OBC_PERIOD", CS%update_OBC_period, & + "The period between recalculations OBC updates. "//& + "If DTBT_RESET_PERIOD 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. @@ -2541,7 +2563,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call VarMix_init(Time, G, GV, US, param_file, diag, CS%VarMix) call set_visc_init(Time, G, GV, US, param_file, diag, CS%visc, CS%set_visc_CSp, restart_CSp, CS%OBC) call thickness_diffuse_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%thickness_diffuse_CSp) - if (CS%split) then allocate(eta(SZI_(G),SZJ_(G))) ; eta(:,:) = 0.0 call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, & @@ -2562,6 +2583,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 + !OBC hack + 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 c41068e1fa..ab4a59db89 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -66,6 +66,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 @@ -82,7 +83,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 @@ -329,7 +330,6 @@ module MOM_open_boundary 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. ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1911,8 +1911,8 @@ 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) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure +subroutine setup_OBC_tracer_reservoirs(Gke, OBC) + integer, intent(in) :: Gke !< Ocean grid ke type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables type(OBC_segment_type), pointer :: segment => NULL() @@ -1925,7 +1925,7 @@ subroutine setup_OBC_tracer_reservoirs(G, OBC) I = segment%HI%IsdB do m=1,OBC%ntr if (associated(segment%tr_Reg%Tr(m)%tres)) then - do k=1,G%ke + do k=1,Gke do j=segment%HI%jsd,segment%HI%jed OBC%tres_x(I,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k) enddo @@ -1936,7 +1936,7 @@ subroutine setup_OBC_tracer_reservoirs(G, OBC) J = segment%HI%JsdB do m=1,OBC%ntr if (associated(segment%tr_Reg%Tr(m)%tres)) then - do k=1,G%ke + do k=1,Gke do i=segment%HI%isd,segment%HI%ied OBC%tres_y(i,J,k,m) = segment%tr_Reg%Tr(m)%t(i,J,k) enddo @@ -3618,6 +3618,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (.not. associated(OBC)) return + call MOM_error(NOTE,"update_OBC_segment_data: called! ") do n = 1, OBC%number_of_segments segment => OBC%segment(n) @@ -4537,7 +4538,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%ke, OBC) !This will redo the T&S end subroutine fill_obgc_segments subroutine fill_temp_salt_segments(G, OBC, tv) @@ -4596,7 +4597,7 @@ subroutine fill_temp_salt_segments(G, OBC, tv) segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) enddo - call setup_OBC_tracer_reservoirs(G, OBC) + !call setup_OBC_tracer_reservoirs(G%ke, OBC) end subroutine fill_temp_salt_segments !> Find the region outside of all open boundary segments and @@ -4856,7 +4857,7 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart ! Local variables type(vardesc) :: vd(2) integer :: m, n - character(len=100) :: mesg + character(len=256) :: mesg,longname type(OBC_segment_type), pointer :: segment=>NULL() if (.not. associated(OBC)) & @@ -4924,37 +4925,41 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart ! Still painfully inefficient, now in four dimensions. ! Allocating both for now so that the pass_vector works. - if (any(OBC%tracer_x_reservoirs_used) .or. any(OBC%tracer_y_reservoirs_used)) then + if (any(OBC%tracer_x_reservoirs_used)) then allocate(OBC%tres_x(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke,OBC%ntr)) OBC%tres_x(:,:,:,:) = 0.0 do m=1,OBC%ntr if (OBC%tracer_x_reservoirs_used(m)) then - if (modulo(HI%turns, 2) /= 0) then - write(mesg,'("tres_y_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') - call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CSp) - else +! if (modulo(HI%turns, 2) /= 0) then +! write(mesg,'("tres_y_",I3.3)') m +! longname="Tracer concentration for NS OBCs for "//trim(Reg%Tr(m)%name) +! vd(1) = var_desc(mesg,"Conc", longname,'v','L') +! call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CSp) +! else write(mesg,'("tres_x_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') + longname="Tracer concentration for EW OBCs for "//trim(Reg%Tr(m)%name) + vd(1) = var_desc(mesg,"Conc", longname,'u','L') call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CSp) - endif +! endif endif enddo -! endif -! if (any(OBC%tracer_y_reservoirs_used)) then + endif + if (any(OBC%tracer_y_reservoirs_used)) then allocate(OBC%tres_y(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke,OBC%ntr)) OBC%tres_y(:,:,:,:) = 0.0 do m=1,OBC%ntr if (OBC%tracer_y_reservoirs_used(m)) then - if (modulo(HI%turns, 2) /= 0) then - write(mesg,'("tres_x_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') - call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CSp) - else +! if (modulo(HI%turns, 2) /= 0) then +! write(mesg,'("tres_x_",I3.3)') m +! longname="Tracer concentration for EW OBCs for "//trim(Reg%Tr(m)%name) +! vd(1) = var_desc(mesg,"Conc", longname,'u','L') +! call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CSp) +! else write(mesg,'("tres_y_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') + longname="Tracer concentration for NS OBCs for "//trim(Reg%Tr(m)%name) + vd(1) = var_desc(mesg,"Conc", longname,'v','L') call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CSp) - endif +! endif endif enddo endif @@ -5005,7 +5010,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) u_L_in = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_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) + & + segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + & !NOT reproduce across restart + !segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(OBC%tres_x(I,j,k,m) + & !NOT reproduce + !segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%t(I,j,k) + & !Reproduce (u_L_out*Reg%Tr(m)%t(I+ishift,j,k) - & u_L_in*segment%tr_Reg%Tr(m)%t(I,j,k))) if (associated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = segment%tr_Reg%Tr(m)%tres(I,j,k) @@ -5030,7 +5037,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) v_L_in = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_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) + & + segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + & !NOT reproduce across restart + !segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(OBC%tres_y(i,J,k,m) + & !NOT reproduce + !segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%t(i,J,k) + & !Reproduce (v_L_out*Reg%Tr(m)%t(i,J+jshift,k) - & v_L_in*segment%tr_Reg%Tr(m)%t(i,J,k))) if (associated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = segment%tr_Reg%Tr(m)%tres(i,J,k) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index b613648c7c..a08d4883c8 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 @@ -380,9 +380,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, OBC, tv) - + call setup_OBC_tracer_reservoirs(G%ke, OBC) + endif ! The thicknesses in halo points might be needed to initialize the velocities. if (new_sim) call pass_var(h, G%Domain)