From 9d5a320db71309220061cf2a5dfbfa0c9e7e81c1 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 10 Nov 2022 04:33:34 -0500 Subject: [PATCH] Obc setup plus segment update period (#198) * Setup OBC segments for COBALT/OBGC tracers - These are updates required to setup OBC segments for OBGC tracers. - Since COBALT package has more than 50 tracers using the MOM6 table mechanism for setting up OBC segments is not feasible. Rather, this update delegates such setup to mechanims used in ocean_BGS tracers leaving MOM6 mechanism for native tracers intact. - Fixed issues caught by MOM6 githubCI * Add capability to change obc segment update period - COBALT tracers do not need as frequent segment bc updates and can use a larger update period to speed up the model. This commit introduces a new parameter DT_OBC_SEG_UPDATE_OBGC that can be adjusted for obc segment update period. - This commit applies the change only to BGC tracers but can easily be changed to apply for all. * Insert missing US%T_to_sec - The unit conversion factor was missing causing a crash in a newer test. * Updates from Andrew Ross - Avoid low initial values in the tracer reservoirs * Per Andrew Ross review * corrected indentation per review * Avoid using module vars per review request - Reviewer asked to avoid using module variables with "save" attributes. - This commit hides the module variables inside the existing OBC type. * Coding style corrections per review * Modification per review: do_not_log if .not.associated(CS%OBC) Co-authored-by: Robert Hallberg --- .../GFDL_ocean_BGC/generic_tracer_utils.F90 | 14 + src/core/MOM.F90 | 34 ++- src/core/MOM_open_boundary.F90 | 246 +++++++++++++++++- src/tracer/MOM_generic_tracer.F90 | 62 ++++- src/tracer/MOM_tracer_flow_control.F90 | 24 +- 5 files changed, 358 insertions(+), 22 deletions(-) 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 cbc310eb7d..e0d3b1d6a9 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 !! @@ -284,6 +287,17 @@ 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 + !> get obc segment properties for each tracer + subroutine g_tracer_get_obc_segment_props(g_tracer_list, name, obc_has, src_file, src_var_name,lfac_in,lfac_out) + type(g_tracer_type), pointer :: g_tracer_list !< pointer to the head of the generic tracer list + character(len=*), intent(in) :: name !< tracer name + logical, intent(out):: obc_has !< .true. if This tracer has OBC + real, optional,intent(out):: lfac_in !< OBC reservoir inverse lengthscale factor + real, optional,intent(out):: lfac_out !< OBC reservoir inverse lengthscale factor + character(len=*),optional,intent(out):: src_file !< OBC source file + character(len=*),optional,intent(out):: src_var_name !< OBC source variable in file + end subroutine g_tracer_get_obc_segment_props + !>Vertical Diffusion of a tracer node !! !! This subroutine solves a tridiagonal equation to find and set values of vertically diffused field diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c61f130ef7..fef71ab4d5 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -133,7 +133,7 @@ module MOM use MOM_tracer_registry, only : lock_tracer_registry, tracer_registry_end use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_CS use MOM_tracer_flow_control, only : tracer_flow_control_init, call_tracer_surface_state -use MOM_tracer_flow_control, only : tracer_flow_control_end +use MOM_tracer_flow_control, only : tracer_flow_control_end, call_tracer_register_obc_segments use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling @@ -294,7 +294,11 @@ module MOM !! barotropic time step [s]. If this is negative dtbt is never !! 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. + type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated. + 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] real, dimension(:,:), pointer :: mass_shelf => NULL() !< Mass of ice shelf [R Z ~> kg m-2] @@ -1132,6 +1136,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 + 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, @@ -2152,6 +2167,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & 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. "//& + "This must be an integer multiple of DT and DT_THERM. "//& + "The default is set to DT.", & + units="s", default=US%T_to_s*CS%dt, do_not_log=.not.associated(CS%OBC)) + ! This is here in case these values are used inappropriately. use_frazil = .false. ; bound_salinity = .false. CS%tv%P_Ref = 2.0e7*US%kg_m3_to_R*US%m_s_to_L_T**2 @@ -2627,6 +2649,8 @@ 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, US, CS%OBC, CS%tracer_Reg, param_file) + !This is the equivalent call to register_temp_salt_segments for external tracers with OBC + call call_tracer_register_obc_segments(GV, param_file, CS%tracer_flow_CSp, CS%tracer_Reg, CS%OBC) ! This needs the number of tracers and to have called any code that sets whether ! reservoirs are used. @@ -2962,6 +2986,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(US%T_to_s*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 1cc8505d17..7d89d97c0f 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -57,7 +57,11 @@ 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 set_obgc_segments_props +public setup_OBC_tracer_reservoirs public open_boundary_register_restarts public update_segment_tracer_reservoirs public update_OBC_ramp @@ -78,7 +82,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=8) :: name !< a name identifier for the segment data + character(len=32) :: name !< a name identifier for the segment data + character(len=8) :: genre !< an identifier for the segment data real :: scale !< A scaling factor for converting input data to !! the internal units of this field real, allocatable :: buffer_src(:,:,:) !< buffer for segment data located at cell faces @@ -91,6 +96,10 @@ module MOM_open_boundary !! The values for tracers should have the same units as the field !! they are being applied to? 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. @@ -262,6 +271,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 @@ -304,6 +315,9 @@ module MOM_open_boundary ! Which segment object describes the current point. integer, allocatable :: segnum_u(:,:) !< Segment number of u-points. integer, allocatable :: segnum_v(:,:) !< Segment number of v-points. + ! Keep the OBC segment properties for external BGC tracers + type(external_tracers_segments_props), pointer :: obgc_segments_props => NULL() !< obgc segment properties + integer :: num_obgc_tracers = 0 !< The total number of obgc tracers ! The following parameters are used in the baroclinic radiation code: real :: gamma_uv !< The relative weighting for the baroclinic radiation @@ -370,6 +384,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() !< pointer to the next node + character(len=128) :: tracer_name !< tracer name + character(len=128) :: tracer_src_file !< tracer source file for BC + character(len=128) :: tracer_src_field !< name of the field in source file to extract BC + real :: lfac_in !< multiplicative factor for inbound tracer reservoir length scale + real :: lfac_out !< multiplicative factor for outbound tracer reservoir length scale +end type external_tracers_segments_props integer :: id_clock_pass !< A CPU time clock character(len=40) :: mdl = "MOM_open_boundary" !< This module's name. @@ -704,7 +727,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) type(ocean_OBC_type), target, 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=1024) :: segstr character(len=256) :: filename character(len=20) :: segnam, suffix @@ -721,6 +744,7 @@ subroutine initialize_segment_data(G, GV, US, 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 @@ -772,8 +796,9 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) cycle ! cycle to next segment endif - allocate(segment%field(num_fields)) - segment%num_fields = num_fields + !There are OBC%num_obgc_tracers obgc tracers are there that are not listed in param file + segment%num_fields = num_fields + OBC%num_obgc_tracers + allocate(segment%field(segment%num_fields)) segment%temp_segment_data_exists = .false. segment%salt_segment_data_exists = .false. @@ -786,9 +811,28 @@ subroutine initialize_segment_data(G, GV, US, 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), m, trim(fields(m)), & - value, filename, fieldname) + obgc_segments_props_list => OBC%obgc_segments_props !pointer to the head node + do m=1,segment%num_fields + if (m .le. num_fields) then + !These are tracers with segments specified in MOM6 style override files + call parse_segment_data_str(trim(segstr), m, trim(fields(m)), value, filename, fieldname) + else + !These are obgc tracers with segments specified by external modules. + !Set a flag so that these can be distinguished from native tracers as they may need + !extra steps for preparation and handling. + segment%field(m)%genre = 'obgc' + !Query the obgc segment properties by traversing the linkedlist + 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) + !Make sure the obgc tracer is not specified in the MOM6 param file too. + 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 @@ -1737,7 +1781,7 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) logical, intent(in) :: use_temperature !< If true, T and S are used ! Local variables - integer :: n,m,num_fields + integer :: n,m,num_fields,na character(len=1024) :: segstr character(len=256) :: filename character(len=20) :: segnam, suffix @@ -1789,6 +1833,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 reservoir index for non-native tracers from 3, hence na=2 below. + !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 + na=2 !number of native MOM6 tracers (T&S) with reservoirs + do m=1,OBC%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+na) = .true. + else + OBC%tracer_y_reservoirs_used(m+na) = .true. + endif + enddo enddo return @@ -3491,6 +3552,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) @@ -3715,7 +3792,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) type(time_type), intent(in) :: Time !< Model time ! Local variables integer :: c, 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 type(OBC_segment_type), pointer :: segment => NULL() integer, dimension(4) :: siz real, dimension(:,:,:), pointer :: tmp_buffer_in => NULL() ! Unrotated input [various units] @@ -3810,6 +3887,10 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) allocate(h_stack(GV%ke), source=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) @@ -4173,6 +4254,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 @@ -4359,6 +4442,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)%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 "//trim(segment%field(m)%name)) + endif + 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 + 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 @@ -4660,6 +4762,123 @@ subroutine register_temp_salt_segments(GV, US, 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(OBC,tr_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + type(ocean_OBC_type),pointer :: OBC !< Open boundary structure + 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 !< factors for tracer reservoir length scales + real, intent(in) :: lfac_out !< factors for tracer reservoir length scales + + type(external_tracers_segments_props),pointer :: node_ptr => NULL() !pointer to type that keeps + ! the tracer segment properties + 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 => OBC%obgc_segments_props + OBC%obgc_segments_props => node_ptr + OBC%num_obgc_tracers = OBC%num_obgc_tracers+1 +end subroutine set_obgc_segments_props + +!> 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 !< pointer to tracer segment properties + 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 !< multiplicative factor for inbound reservoir length scale + real, intent(out) :: lfac_out !< multiplicative factor for outbound 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 + 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 + 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 + 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 + +subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical 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 + real :: I_scale + + if (.not. associated(OBC)) return + 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 + I_scale = 1.0 + if (segment%tr_Reg%Tr(nt)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(nt)%scale + ! 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 + OBC%tres_x(I,j,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%t(I,j,k) + 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 + OBC%tres_y(i,J,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%t(i,J,k) + enddo ; enddo + endif + segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) + enddo +end subroutine fill_obgc_segments + subroutine fill_temp_salt_segments(G, GV, US, OBC, tv) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -5109,7 +5328,6 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ! e.g. a_in is -1 only if b_in ==1 and uhr or vhr is inward ! e.g. a_out is 1 only if b_out==1 and uhr or vhr is outward ! It's clear that a_in and a_out cannot be both non-zero [nodim] - nz = GV%ke ntr = Reg%ntr @@ -5140,9 +5358,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ! When InvLscale_in is 0 and inflow, only nudged data is applied to reservoirs a_out = b_out * max(0.0, sign(1.0, idir*uhr(I,j,k))) a_in = b_in * min(0.0, sign(1.0, idir*uhr(I,j,k))) - 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 - (a_out - a_in)) + ((u_L_out + a_out) - (u_L_in + a_in)) segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1) * & @@ -5171,9 +5389,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz a_out = b_out * max(0.0, sign(1.0, jdir*vhr(i,J,k))) a_in = b_in * min(0.0, sign(1.0, jdir*vhr(i,J,k))) - 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) fac1 = (1.0 - (a_out - a_in)) + ((v_L_out + a_out) - (v_L_in + a_in)) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 3cbed68467..4e88944958 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -27,6 +27,7 @@ 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_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here @@ -39,6 +40,8 @@ module MOM_generic_tracer use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, MOM_read_data, slasher use MOM_open_boundary, only : ocean_OBC_type + use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments + use MOM_open_boundary, only : set_obgc_segments_props use MOM_restart, only : register_restart_field, query_initialized, set_initialized, MOM_restart_CS use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP use MOM_sponge, only : set_up_sponge_field, sponge_CS @@ -65,6 +68,7 @@ module MOM_generic_tracer public MOM_generic_flux_init public MOM_generic_tracer_min_max public MOM_generic_tracer_fluxes_accumulate + public register_MOM_generic_tracer_segments !> Control structure for generic tracers type, public :: MOM_generic_tracer_CS ; private @@ -79,7 +83,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() !< Restart control structure - + type(ocean_OBC_type), pointer :: OBC => NULL() ! Pointer to the first element of the linked list of generic tracers. type(g_tracer_type), pointer :: g_tracer_list => NULL() @@ -98,10 +102,9 @@ 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), target, intent(inout) :: restart_CS !< MOM restart control struct - ! Local variables logical :: register_MOM_generic_tracer - + logical :: obc_has ! This include declares and sets the variable "version". # include "version_variable.h" @@ -112,6 +115,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) integer :: ntau, 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 :: 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 @@ -156,7 +161,6 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) CS%restart_CSp => restart_CS - ntau=1 ! MOM needs the fields at only one time step @@ -216,6 +220,52 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) register_MOM_generic_tracer = .true. end function register_MOM_generic_tracer + !> Register OBC segments for generic tracers + subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file) + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. + type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer + !! advection and diffusion module. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + ! Local variables + logical :: obc_has + ! This include declares and sets the variable "version". +# include "version_variable.h" + + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer_segments' + type(g_tracer_type), pointer :: g_tracer,g_tracer_next + character(len=fm_string_len) :: g_tracer_name + character(len=fm_string_len) :: obc_src_file_name,obc_src_field_name + real :: lfac_in,lfac_out + + if (.NOT. associated(OBC)) return + !Get the tracer list + call generic_tracer_get_list(CS%g_tracer_list) + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + + g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,g_tracer_name) + if (g_tracer_is_prog(g_tracer)) then + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,& + obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + if (obc_has) then + call set_obgc_segments_props(OBC,g_tracer_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + call register_obgc_segments(GV, OBC, tr_Reg, param_file, g_tracer_name) + endif + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + + end subroutine register_MOM_generic_tracer_segments !> Initialize phase II: Initialize required variables for generic tracers !! There are some steps of initialization that cannot be done in register_MOM_generic_tracer !! This is the place and time to do them: @@ -244,7 +294,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 @@ -348,6 +398,8 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, call set_initialized(tr_ptr, g_tracer_name, CS%restart_CSp) 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, GV, 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 1345126d73..af0dded244 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -60,6 +60,7 @@ module MOM_tracer_flow_control use MOM_generic_tracer, only : MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state use MOM_generic_tracer, only : end_MOM_generic_tracer, MOM_generic_tracer_get, MOM_generic_flux_init use MOM_generic_tracer, only : MOM_generic_tracer_stock, MOM_generic_tracer_min_max, MOM_generic_tracer_CS +use MOM_generic_tracer, only : register_MOM_generic_tracer_segments use pseudo_salt_tracer, only : register_pseudo_salt_tracer, initialize_pseudo_salt_tracer use pseudo_salt_tracer, only : pseudo_salt_tracer_column_physics, pseudo_salt_tracer_surface_state use pseudo_salt_tracer, only : pseudo_salt_stock, pseudo_salt_tracer_end, pseudo_salt_tracer_CS @@ -75,6 +76,7 @@ module MOM_tracer_flow_control public call_tracer_register, tracer_flow_control_init, call_tracer_set_forcing public call_tracer_column_fns, call_tracer_surface_state, call_tracer_stocks public call_tracer_flux_init, get_chl_from_model, tracer_flow_control_end +public call_tracer_register_obc_segments !> The control structure for orchestrating the calling of tracer packages type, public :: tracer_flow_control_CS ; private @@ -114,6 +116,7 @@ module MOM_tracer_flow_control contains + !> This subroutine carries out a series of calls to initialize the air-sea !! tracer fluxes, but it does not record the generated indicies, and it may !! be called _before_ the ocean model has been initialized and may be called @@ -163,7 +166,6 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) type(MOM_restart_CS), intent(inout) :: restart_CS !< A pointer to the restart control !! structure. - ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_tracer_flow_control" ! This module's name. @@ -357,6 +359,26 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag end subroutine tracer_flow_control_init +!> This subroutine calls all registered tracers to register their OBC segments +!! similar to register_temp_salt_segments for T&S +subroutine call_tracer_register_obc_segments(GV, param_file, CS, tr_Reg, OBC) + 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. + type(tracer_flow_control_CS), pointer :: CS !< A pointer that is set to point to the + !! control structure for this module. + type(tracer_registry_type), pointer :: tr_Reg !< A pointer that is set to point to the + !! control structure for the tracer + !! advection and diffusion module. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition + !! type specifies whether, where, + !! and what open boundary + !! conditions are used. + + if (CS%use_MOM_generic_tracer) & + call register_MOM_generic_tracer_segments(CS%MOM_generic_tracer_CSp, GV, OBC, tr_Reg, param_file) + +end subroutine call_tracer_register_obc_segments !> This subroutine extracts the chlorophyll concentrations from the model state, if possible subroutine get_chl_from_model(Chl_array, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.