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 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ba784f4a40..e5fde8082d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -278,6 +278,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 :: 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 :: & @@ -1032,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 @@ -1078,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, @@ -1996,6 +2010,14 @@ 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=CS%dt) + + ! 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 @@ -2400,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, & @@ -2433,7 +2450,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(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(HI, GV, CS%OBC, CS%tracer_Reg, & @@ -2728,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 2c3f016005..6399757089 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -54,13 +54,17 @@ 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 open_boundary_register_restarts public update_segment_tracer_reservoirs public update_OBC_ramp 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 @@ -77,14 +81,19 @@ 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, 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 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. @@ -202,13 +211,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 :: 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 @@ -237,6 +246,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 @@ -331,6 +342,16 @@ 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 + 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 integer :: id_clock_pass !< A CPU time clock character(len=40) :: mdl = "MOM_open_boundary" !< This module's name. @@ -647,7 +668,7 @@ subroutine initialize_segment_data(G, 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 @@ -664,6 +685,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 @@ -715,8 +737,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. @@ -729,9 +752,24 @@ 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), m, trim(fields(m)), & - value, filename, 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), m, trim(fields(m)), & + value, filename, fieldname) + else + segment%field(m)%genre='obgc' + 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()) & + 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 @@ -1676,8 +1714,8 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) 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 @@ -1689,6 +1727,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 @@ -2067,20 +2122,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, GV, OBC) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure +subroutine setup_OBC_tracer_reservoirs(GV, OBC, skip_ts) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(ocean_OBC_type), target, intent(inout) :: 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 (allocated(segment%tr_Reg%Tr(m)%tres)) then do k=1,GV%ke do j=segment%HI%jsd,segment%HI%jed @@ -2091,7 +2153,7 @@ subroutine setup_OBC_tracer_reservoirs(G, GV, OBC) enddo else J = segment%HI%JsdB - do m=1,OBC%ntr + do m=m0,OBC%ntr if (allocated(segment%tr_Reg%Tr(m)%tres)) then do k=1,GV%ke do i=segment%HI%isd,segment%HI%ied @@ -3481,6 +3543,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) @@ -3708,7 +3786,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, 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() @@ -3794,6 +3872,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) @@ -4138,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 @@ -4325,6 +4409,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 @@ -4591,9 +4694,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, & @@ -4606,6 +4706,115 @@ 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,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, 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 + 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 + +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 + + 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 + ! 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(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) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -4662,7 +4871,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(G, 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 @@ -5027,12 +5236,29 @@ 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 + if (segment%is_E_or_W) then I = segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed @@ -5046,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 / & - ((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 / & - ((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 @@ -5071,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 / & - ((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 / & - ((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 @@ -5401,7 +5628,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/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) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index bf9f01e266..d87b43f7cd 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -27,7 +27,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_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS use MOM_coms, only : max_across_PEs, min_across_PEs, PE_here use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -38,7 +39,8 @@ module MOM_generic_tracer use MOM_grid, only : ocean_grid_type 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 : ocean_OBC_type, register_obgc_segments, fill_obgc_segments + use MOM_open_boundary, only : set_obgc_segments_props use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_spatial_means, only : global_area_mean use MOM_sponge, only : set_up_sponge_field, sponge_CS @@ -79,7 +81,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() @@ -95,7 +97,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 @@ -103,10 +105,11 @@ 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 - + 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 - + 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? @@ -114,6 +117,8 @@ 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 :: 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 @@ -158,7 +163,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 @@ -204,6 +209,13 @@ 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 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(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 call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & restart_CS, longname=longname, units=units) @@ -228,7 +240,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. @@ -239,15 +251,13 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), 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 !! 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 @@ -349,6 +359,8 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, endif 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, 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 2ae72a3270..323cdf0716 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -147,7 +147,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 @@ -160,7 +160,10 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) !! advection and diffusion module. type(MOM_restart_CS), intent(inout) :: 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" @@ -257,7 +260,7 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) tr_Reg, restart_CS) 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) if (CS%use_pseudo_salt_tracer) CS%use_pseudo_salt_tracer = & register_pseudo_salt_tracer(HI, GV, param_file, CS%pseudo_salt_tracer_CSp, & tr_Reg, restart_CS) @@ -340,7 +343,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag call initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS%CFC_cap_CSp) 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) if (CS%use_pseudo_salt_tracer) & call initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS%pseudo_salt_tracer_CSp, & diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 2c77df3e74..d5c68061fd 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -852,6 +852,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