diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 74fc6f7ed0..256f203a65 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -118,7 +118,7 @@ module MOM use MOM_open_boundary, only : ocean_OBC_type, open_boundary_end use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : read_OBC_dynamics_data, read_OBC_tracer_data -use MOM_open_boundary, only : initialize_OBC_segment_reservoirs +use MOM_open_boundary, only : initialize_OBC_tracer_reservoirs use MOM_open_boundary, only : setup_OBC_tracer_reservoirs use MOM_open_boundary, only : setup_OBC_thickness_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields @@ -127,6 +127,8 @@ module MOM use MOM_open_boundary, only : rotate_OBC_config use MOM_open_boundary, only : open_boundary_halo_update, write_OBC_info, chksum_OBC_segments use MOM_open_boundary, only : segment_thickness_reservoir_init +use MOM_open_boundary, only : fill_temp_salt_segments, fill_thickness_segments +use MOM_open_boundary, only : open_boundary_test_extern_h use MOM_open_boundary, only : copy_OBC_radiation_coefs use MOM_open_boundary, only : copy_OBC_tracer_reservoirs, copy_OBC_thickness_reservoirs use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init @@ -136,7 +138,7 @@ module MOM use MOM_set_visc, only : set_visc_init, set_visc_end use MOM_shared_initialization, only : write_ocean_geometry_file use MOM_sponge, only : init_sponge_diags, sponge_CS -use MOM_state_initialization, only : MOM_initialize_state, MOM_initialize_OBCs +use MOM_state_initialization, only : MOM_initialize_state, initialize_user_OBCs use MOM_stoch_eos, only : MOM_stoch_eos_init, MOM_stoch_eos_run, MOM_stoch_eos_CS use MOM_stoch_eos, only : stoch_EOS_register_restarts, post_stoch_EOS_diags, mom_calc_varT use MOM_sum_output, only : write_energy, accumulate_net_input @@ -2412,6 +2414,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & logical :: calc_dtbt ! Indicates whether the dynamically adjusted barotropic ! time step needs to be updated before it is used. logical :: debug_truncations ! If true, turn on diagnostics useful for debugging truncations. + logical :: obc_debug_test ! If true, call open_boundary_test_extern_uv. integer :: first_direction ! An integer that indicates which direction is to be ! updated first in directionally split parts of the ! calculation. @@ -2424,6 +2427,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & logical :: semi_Boussinesq ! If true, this run is partially non-Boussinesq logical :: use_KPP ! If true, diabatic is using KPP vertical mixing logical :: MLE_use_PBL_MLD ! If true, use stored boundary layer depths for submesoscale restratification. + logical :: OBC_TS_reservoir_init_bug logical :: OBC_reservoir_init_bug logical :: OBC_bgc_time_ref_bug ! If true, use the start of the current run (not the overall ! start time) as the reference for OBC BGC tracer update schedule. @@ -2622,6 +2626,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call get_param(param_file, "MOM", "DEBUG_OBCS", CS%debug_OBCs, & "If true, write out verbose debugging data about OBCs.", & default=.false., debuggingParam=.true., do_not_log=(number_of_OBC_segments<=0)) + call get_param(param_file, "MOM", "OBC_DEBUGGING_TESTS", obc_debug_test, & + "If true, do additional calls resetting values to help verify the correctness "//& + "of the open boundary condition code.", default=.false., do_not_log=.true., & + old_name="DEBUG_OBC", debuggingParam=.true.) call get_param(param_file, "MOM", "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & "If true, the defaults for certain recently added bug-fix flags are set to "//& "recreate the bugs so that the code can be moved forward without changing "//& @@ -2631,7 +2639,15 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & "means that bugs are only used if they are actively selected, but it also "//& "means that answers may change when code is updated due to newly found bugs.", & default=.true.) - + ! Log this parameter in MOM_initialize_state + call get_param(param_file, "MOM", "OBC_TS_RESERVOIR_INIT_BUG", OBC_TS_reservoir_init_bug, & + "If true, set the OBC temperature and salinity reservoirs at the startup of a "//& + "new run from initial values that are set before remapping.", default=enable_bugs, & + do_not_log=.true.) + call get_param(param_file, "MOM", "OBC_RESERVOIR_INIT_BUG", OBC_reservoir_init_bug, & + "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& + "interior tracer concentrations regardless of properties that may be explicitly "//& + "specified for the reservoir concentrations.", default=enable_bugs) call get_param(param_file, "MOM", "DT", CS%dt, & "The (baroclinic) dynamics time step. The time-step that "//& "is actually used will be an integer fraction of the "//& @@ -3311,11 +3327,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%tv%S => S_in if (associated(CS%OBC)) then - ! Log this parameter in MOM_initialize_state - call get_param(param_file, "MOM", "OBC_RESERVOIR_INIT_BUG", OBC_reservoir_init_bug, & - "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& - "interior tracer concentrations regardless of properties that may be explicitly "//& - "specified for the reservoir concentrations.", default=enable_bugs, do_not_log=.true.) if (OBC_reservoir_init_bug .and. (allocated(CS%OBC%tres_x) .or. allocated(CS%OBC%tres_y))) & call MOM_error(FATAL, "OBC_RESERVOIR_INIT_BUG can not be set to true with grid rotation.") endif @@ -3414,21 +3425,46 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif if (associated(CS%OBC)) then - call MOM_initialize_OBCs(CS%h, CS%tv, CS%OBC, Time, G, GV, US, param_file, restart_CSp, CS%tracer_Reg) - - if (use_temperature) then - call pass_var(CS%tv%T, G%Domain, complete=.false.) - call pass_var(CS%tv%S, G%Domain, complete=.true.) - endif + ! Read and calculate segment dynamic fields. Needed for + ! (a) continuity calls in initialize_dyn_* + ! (b) segment%dz used for subsequent read_OBC_tracer_data call + ! (c) One-time dynamic update when OBC_VALUE_UPDATE_BUG=True + ! Calcualte tv%SpV_avg for dz used in read_OBC_segment_data call calc_derived_thermo(CS%tv, CS%h, G, GV, US) - - ! Call this during initialization to fill boundary arrays from fixed values + ! OBC file -> segment%field()%buffer_dst call read_OBC_dynamics_data(G, GV, US, CS%OBC, CS%tv, CS%h, Time) + ! segment%field()%buffer_dst -> segment%normal_trans etc call update_OBC_dynamics_data(G, GV, US, CS%OBC, CS%h, Time) - ! BGC data is not read/updated at initialization since OBC%update_OBC_seg_data is false. - call read_OBC_tracer_data(G, GV, US, CS%OBC, Time, include_bgc=.false.) - call update_OBC_tracer_data(CS%OBC, include_bgc=.false.) - call initialize_OBC_segment_reservoirs(GV, CS%OBC) + ! h -> %h and %h_res + if (CS%OBC%use_h_res) & + call fill_thickness_segments(G, GV, US, CS%OBC, CS%h) + + ! Initialize OBC segment T/S reservoir [segment%tr_Reg%Tr()%tres] + if (is_new_run(restart_CSp) .and. use_temperature) then + ! T/S reservoirs are initialized by interior data. + ! They may be overriden by external data again a few lines later. The only reason that this + ! call is made unconditional to OBC_reservoir_init_bug or diabatic_first is to provide a + ! fallback for buggy cases where use_temperature is True but external T/S are not given. + if (.not. OBC_TS_reservoir_init_bug) & + ! Interior T/S -> segment%tr_Reg%Tr()%t and segment%tr_Reg%Tr()%tres + call fill_temp_salt_segments(G, GV, US, CS%OBC, CS%tv) + + ! T/S reservoirs are initialized by external data. + ! OBC file -> segment%field()%buffer_dst -> segment%tr_Reg%Tr()%t -> segment%tr_Reg%Tr()%tres. + ! Note BGC data is deliberately not read/updated at initialization since OBC%update_OBC_seg_data is false. + if ((.not. OBC_reservoir_init_bug) .or. CS%diabatic_first) then + ! OBC file -> segment%field()%buffer_dst + call read_OBC_tracer_data(G, GV, US, CS%OBC, Time, include_bgc=.false.) + ! [T/S] segment%field()%buffer_dst -> segment%tr_Reg%Tr()%t + call update_OBC_tracer_data(CS%OBC, include_bgc=.false.) + ! [T/S] segment%tr_Reg%Tr()%t -> segment%tr_Reg%Tr()%tres + call initialize_OBC_tracer_reservoirs(CS%OBC) + endif + endif + + call initialize_user_OBCs(CS%tv, CS%OBC, G, GV, US, param_file, CS%tracer_Reg) + + if (obc_debug_test) call open_boundary_test_extern_h(G, GV, CS%OBC, CS%h) endif if (use_ice_shelf .and. CS%debug) then @@ -3840,19 +3876,20 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call register_diags_offline_transport(Time, CS%diag, CS%offline_CSp, GV, US) endif + ! [if new_sim] Save %tres to OBC%tres_[xy] (needs to be after tracer_flow_control_init to init BGC) + ! [if restart] copy OBC%tres_[xy] to %tres if (associated(CS%OBC)) then - ! At this point any information related to the tracer reservoirs has either been read from - ! the restart file or has been specified in the segments. Initialize the tracer reservoir - ! values from the segments if they have not been set via the restart file. - call setup_OBC_tracer_reservoirs(G, GV, CS%OBC, restart_CSp) - call setup_OBC_thickness_reservoirs(G, GV, CS%OBC, restart_CSp) - call open_boundary_halo_update(G, CS%OBC) - call copy_OBC_radiation_coefs(CS%OBC) - if (.not. (CS%OBC%reservoir_init_bug .and. new_sim .and. CS%diabatic_first)) & - ! The if-guard is needed to preserve old answers with OBC_RESERVOIR_INIT_BUG=True, in which case - ! segment T/S reservoir %tres and global restart arrays OBC%tres_x/y have diverged at this point. + if (is_new_run(restart_CSp)) then + if (save_IC) then + call setup_OBC_tracer_reservoirs(G, GV, CS%OBC, restart_CSp) + call setup_OBC_thickness_reservoirs(G, GV, CS%OBC, restart_CSp) + endif + else + call open_boundary_halo_update(G, CS%OBC) + call copy_OBC_radiation_coefs(CS%OBC) call copy_OBC_tracer_reservoirs(CS%OBC) - call copy_OBC_thickness_reservoirs(CS%OBC, G, GV) + call copy_OBC_thickness_reservoirs(CS%OBC, G, GV) + endif endif call register_obsolete_diagnostics(param_file, CS%diag) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index ecd5e51106..71dfbddf57 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -53,7 +53,7 @@ module MOM_open_boundary public read_OBC_tracer_data public update_OBC_dynamics_data public update_OBC_tracer_data -public initialize_OBC_segment_reservoirs +public initialize_OBC_tracer_reservoirs public open_boundary_test_extern_uv public open_boundary_test_extern_h public open_boundary_zero_normal_flow @@ -76,7 +76,6 @@ module MOM_open_boundary public copy_OBC_thickness_reservoirs public update_segment_tracer_reservoirs public update_segment_thickness_reservoirs -public set_initialized_OBC_tracer_reservoirs public update_OBC_ramp public remap_OBC_fields public rotate_OBC_config @@ -155,7 +154,6 @@ module MOM_open_boundary !> Tracer on OBC segment data structure, for putting into a segment tracer registry. type, public :: OBC_segment_tracer_type - logical :: is_initialized !< Reservoir values have been set when True character(len=32) :: name !< Tracer name used for error messages integer :: ntr_index = -1 !< Index of segment tracer in the global tracer registry real, allocatable :: t(:,:,:) !< External tracer concentration array in rescaled @@ -955,6 +953,8 @@ subroutine open_boundary_config(G, US, param_file, OBC) call open_boundary_dealloc(OBC) endif + id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) + end subroutine open_boundary_config !> Setup vertical remapping for open boundaries @@ -2294,13 +2294,12 @@ subroutine open_boundary_halo_update(G, OBC) if (.not.associated(OBC)) return - id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) if (OBC%radiation_BCs_exist_globally) call pass_vector(OBC%rx_normal, OBC%ry_normal, G%Domain, & To_All+Scalar_Pair) if (OBC%oblique_BCs_exist_globally) then -! call pass_vector(OBC%rx_oblique_u, OBC%ry_oblique_v, G%Domain, To_All+Scalar_Pair) -! call pass_vector(OBC%ry_oblique_u, OBC%rx_oblique_v, G%Domain, To_All+Scalar_Pair) -! call pass_vector(OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) + ! call pass_vector(OBC%rx_oblique_u, OBC%ry_oblique_v, G%Domain, To_All+Scalar_Pair) + ! call pass_vector(OBC%ry_oblique_u, OBC%rx_oblique_v, G%Domain, To_All+Scalar_Pair) + ! call pass_vector(OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) call create_group_pass(OBC%pass_oblique, OBC%rx_oblique_u, OBC%ry_oblique_v, G%Domain, To_All+Scalar_Pair) call create_group_pass(OBC%pass_oblique, OBC%ry_oblique_u, OBC%rx_oblique_v, G%Domain, To_All+Scalar_Pair) call create_group_pass(OBC%pass_oblique, OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) @@ -2604,26 +2603,14 @@ subroutine setup_OBC_tracer_reservoirs(G, GV, OBC, restart_CS) if (segment%is_E_or_W .and. set_tres_x) then I = segment%HI%IsdB - if (segment%tr_Reg%Tr(m)%is_initialized) then - do k=1,GV%ke ; do j=segment%HI%jsd,segment%HI%jed - OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,j,k) - enddo ; enddo - else - do k=1,GV%ke ; do j=segment%HI%jsd,segment%HI%jed - OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%t(i,j,k) - enddo ; enddo - endif + do k=1,GV%ke ; do j=segment%HI%jsd,segment%HI%jed + OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,j,k) + enddo ; enddo elseif (segment%is_N_or_S .and. set_tres_y) then J = segment%HI%JsdB - if (segment%tr_Reg%Tr(m)%is_initialized) then - do k=1,GV%ke ; do i=segment%HI%isd,segment%HI%ied - OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) - enddo ; enddo - else - do k=1,GV%ke ; do i=segment%HI%isd,segment%HI%ied - OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%t(i,J,k) - enddo ; enddo - endif + do k=1,GV%ke ; do i=segment%HI%isd,segment%HI%ied + OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) + enddo ; enddo endif endif ; endif enddo @@ -2696,30 +2683,6 @@ subroutine setup_OBC_thickness_reservoirs(G, GV, OBC, restart_CS) end subroutine setup_OBC_thickness_reservoirs -!> Record that the tracer reservoirs have been initialized so that their values are not reset later. -subroutine set_initialized_OBC_tracer_reservoirs(G, OBC, restart_CS) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(ocean_OBC_type), intent(in) :: OBC !< Open boundary control structure - type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure - character(len=12) :: x_var_name, y_var_name - integer :: m - - do m=1,OBC%ntr - ! Set the names of the reservoirs for this tracer in the restart file - if (modulo(G%HI%turns, 2) == 0) then - write(x_var_name,'("tres_x_",I3.3)') m - write(y_var_name,'("tres_y_",I3.3)') m - else - write(x_var_name,'("tres_y_",I3.3)') m - write(y_var_name,'("tres_x_",I3.3)') m - endif - - if (OBC%tracer_x_reservoirs_used(m)) call set_initialized(OBC%tres_x, x_var_name, restart_CS) - if (OBC%tracer_y_reservoirs_used(m)) call set_initialized(OBC%tres_y, y_var_name, restart_CS) - enddo - -end subroutine set_initialized_OBC_tracer_reservoirs - !> Copy radiation and oblique boundary condition coefficients (phase speeds and normalizing !! denominator) from the global restart arrays into the per-segment arrays. subroutine copy_OBC_radiation_coefs(OBC) @@ -4999,25 +4962,22 @@ subroutine update_OBC_tracer_data(OBC, include_bgc) enddo ! end segment loop end subroutine update_OBC_tracer_data -!> Initialize thickness and tracer reservoirs to external value. -subroutine initialize_OBC_segment_reservoirs(GV, OBC) - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - type(ocean_OBC_type), pointer :: OBC !< Open boundary structure +!> Initialize tracer reservoirs to the current external tracer value. +subroutine initialize_OBC_tracer_reservoirs(OBC) + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure ! Local variables type(OBC_segment_type), pointer :: segment => NULL() integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: is_seg, ie_seg, js_seg, je_seg, nz integer :: n, m, nt, i, j, k - character(len=256) :: msg ! Error message if (.not. associated(OBC)) return - nz = GV%ke + nz = OBC%ke do n=1,OBC%number_of_segments segment => OBC%segment(n) - if (.not. segment%on_pe) cycle isd = segment%HI%isd ; ied = segment%HI%ied ; IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB @@ -5031,32 +4991,17 @@ subroutine initialize_OBC_segment_reservoirs(GV, OBC) js_seg = JsdB ; je_seg = JedB ! = js_seg endif - ! Thickness - ! If the thickness reservoir has not yet been initialized, then set to external value. - if (OBC%thickness_x_reservoirs_used .or. OBC%thickness_y_reservoirs_used) then - if (.not. segment%h_Reg%is_initialized) then ! h_Reg may be initialized by fill_thickness_segments - do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg - segment%h_Reg%h_res(i,j,k) = segment%h_Reg%h(i,j,k) - enddo ; enddo ; enddo - segment%h_Reg%is_initialized = .true. - endif - endif - - ! Tracers ! If the tracer reservoir has not yet been initialized, then set to external value. do m=NUM_PHYS_FIELDS-1, segment%num_fields ! F_T = NUM_PHYS_FIELDS-1 and F_S = NUM_PHYS_FIELDS if ((.not. allocated(segment%field(m)%buffer_dst)) .or. & (segment%field(m)%bgc_tracer .and. (.not. OBC%update_OBC_seg_data))) cycle nt = segment%field(m)%tr_index - if (.not. segment%tr_Reg%Tr(nt)%is_initialized) then ! T/S may be initialized by fill_temp_salt_segments - do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg - 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 + do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg + segment%tr_Reg%Tr(nt)%tres(i,j,k) = segment%tr_Reg%Tr(nt)%t(i,j,k) + enddo ; enddo ; enddo enddo ! end tracer field loop enddo ! end segment loop -end subroutine initialize_OBC_segment_reservoirs +end subroutine initialize_OBC_tracer_reservoirs !> Update the OBC ramp value as a function of time. !! If called with the optional argument activate=.true., record the @@ -5270,12 +5215,10 @@ subroutine register_segment_tracer(tr_ptr, ntr_index, param_file, GV, segment, O if (present(OBC_scalar)) then init_value = OBC_scalar - segment%tr_Reg%Tr(ntseg)%is_initialized = .true. segment%tr_Reg%Tr(ntseg)%resrv_lfac_in = 0.0 segment%tr_Reg%Tr(ntseg)%resrv_lfac_out = 0.0 else init_value = 0.0 - segment%tr_Reg%Tr(ntseg)%is_initialized = .false. ! Currently, resrv_lfac_in/out are for BGC tracers only. if (present(resrv_lfac_in)) segment%tr_Reg%Tr(ntseg)%resrv_lfac_in = resrv_lfac_in if (present(resrv_lfac_out)) segment%tr_Reg%Tr(ntseg)%resrv_lfac_out = resrv_lfac_out @@ -5499,31 +5442,7 @@ subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) enddo ; enddo endif - if (.not.segment%tr_Reg%Tr(nt)%is_initialized) & - segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) - - if (OBC%reservoir_init_bug) then - ! OBC%tres_x and OBC%tres_y should not be set here, but in a subsequent call to setup_OBC_tracer_reservoirs. - ! Note that fill_obgc_segments is not called for runs that start from a restart file. - I_scale = 1.0 - if (segment%tr_Reg%Tr(nt)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(nt)%scale - if (segment%is_E_or_W) then - if (allocated(OBC%tres_x)) then - I = segment%HI%IsdB - do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed - OBC%tres_x(I,j,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%tres(I,j,k) - enddo ; enddo - endif - else ! segment%is_N_or_S - if (allocated(OBC%tres_y)) then - J = segment%HI%JsdB - do k=1,nz ; do i=segment%HI%isd,segment%HI%ied - OBC%tres_y(i,J,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%tres(i,J,k) - enddo ; enddo - endif - endif - endif - + segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) enddo ! End of loop over segments. end subroutine fill_obgc_segments @@ -5579,10 +5498,9 @@ subroutine fill_temp_salt_segments(G, GV, US, OBC, tv) endif enddo ; enddo endif - if (.not.segment%tr_Reg%Tr(1)%is_initialized) & - segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:) - if (.not.segment%tr_Reg%Tr(2)%is_initialized) & - segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) + + segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:) + segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) enddo end subroutine fill_temp_salt_segments diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 95e1260459..b3d6421352 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -22,10 +22,8 @@ module MOM_state_initialization use MOM_interface_heights, only : find_eta, dz_to_thickness, dz_to_thickness_simple use MOM_interface_heights, only : calc_derived_thermo use MOM_io, only : file_exists, field_size, MOM_read_data, MOM_read_vector, slasher -use MOM_open_boundary, only : ocean_OBC_type, open_boundary_test_extern_h -use MOM_open_boundary, only : fill_temp_salt_segments, setup_OBC_tracer_reservoirs -use MOM_open_boundary, only : fill_thickness_segments -use MOM_open_boundary, only : set_initialized_OBC_tracer_reservoirs +use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : fill_temp_salt_segments use MOM_restart, only : restore_state, is_new_run, copy_restart_var, copy_restart_vector use MOM_restart, only : restart_registry_lock, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density @@ -100,7 +98,7 @@ module MOM_state_initialization #include -public MOM_initialize_state, MOM_initialize_OBCs +public MOM_initialize_state, initialize_user_OBCs ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -450,7 +448,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call get_param(PF, mdl, "OBC_TS_RESERVOIR_INIT_BUG", OBC_TS_reservoir_init_bug, & "If true, set the OBC temperature and salinity reservoirs at the startup of a "//& "new run from initial values that are set before remapping.", & - default=enable_bugs, do_not_log=.true.) + default=enable_bugs) if (OBC_TS_reservoir_init_bug) then ! These calls should be moved down to join the OBC code, but doing so changes answers because ! the temperatures and salinities can change due to the remapping and reading from the restarts. @@ -653,112 +651,60 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & end subroutine MOM_initialize_state -subroutine MOM_initialize_OBCs(h, tv, OBC, Time, G, GV, US, PF, restart_CS, tracer_Reg) +subroutine initialize_user_OBCs(tv, OBC, G, GV, US, PF, tracer_Reg) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic !! variables type(ocean_OBC_type), pointer :: OBC !< The open boundary condition control structure. - type(time_type), intent(in) :: Time !< Time at the start of the run segment. type(param_file_type), intent(in) :: PF !< A structure indicating the open file to parse !! for model parameter values. - type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure type(tracer_registry_type), pointer :: tracer_Reg !< A pointer to the tracer registry ! Local variables character(len=200) :: config - logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to - ! recreate the bugs, or if false bugs are only used if actively selected. - logical :: debug ! If true, write debugging output. - logical :: debug_obc ! If true, do additional calls resetting values to help debug the correctness - ! of the open boundary condition code. - logical :: OBC_TS_reservoir_init_bug ! If true, set the OBC temperature and salinity reservoirs - ! at the startup of a new run from initial values that are set before remapping. - logical :: OBC_reservoir_init_bug ! If true, set the OBC tracer reservoirs at the startup of a new - ! run from the interior tracer concentrations regardless of properties that - ! may be explicitly specified for the reservoir concentrations. - - call callTree_enter('MOM_initialize_OBCs()') - if (associated(OBC)) then - call get_param(PF, mdl, "DEBUG", debug, default=.false.) - call get_param(PF, mdl, "OBC_DEBUGGING_TESTS", debug_obc, & - "If true, do additional calls resetting values to help verify the correctness "//& - "of the open boundary condition code.", default=.false., & - do_not_log=.true., old_name="DEBUG_OBC", debuggingParam=.true.) - call get_param(PF, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & - default=.true., do_not_log=.true.) ! This is logged from MOM.F90. - call get_param(PF, mdl, "OBC_TS_RESERVOIR_INIT_BUG", OBC_TS_reservoir_init_bug, & - "If true, set the OBC temperature and salinity reservoirs at the startup of a "//& - "new run from initial values that are set before remapping.", default=enable_bugs) - if (associated(tv%T) .and. (.not.OBC_TS_reservoir_init_bug)) then - ! Store the updated temperatures and salinities at the open boundaries, noting that they may - ! still be updated by the calls in the next 50 lines, so the code setting the tracer - ! reservoir values will come later in the calling routine. - call fill_temp_salt_segments(G, GV, US, OBC, tv) - endif - call get_param(PF, mdl, "OBC_RESERVOIR_INIT_BUG", OBC_reservoir_init_bug, & - "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& - "interior tracer concentrations regardless of properties that may be explicitly "//& - "specified for the reservoir concentrations.", default=enable_bugs) - if (OBC_reservoir_init_bug .and. associated(tv%T) .and. is_new_run(restart_CS)) then - ! Set up OBC%trex_x and OBC%tres_y as they have not been read from a restart file. - ! When OBC_RESERVOIR_INIT_BUG is false, setup_OBC_tracer_reservoirs() is called from initialize_MOM - ! after all tracer package initialization is finished and grid rotation has been dealt with. - call setup_OBC_tracer_reservoirs(G, GV, OBC) - ! Ensure that the values of the tracer reservoirs that have just been set will not be revised. - call set_initialized_OBC_tracer_reservoirs(G, OBC, restart_CS) - endif - ! This controls user code for setting open boundary data - call get_param(PF, mdl, "OBC_USER_CONFIG", config, & - "A string that sets how the user code is invoked to set open boundary data: \n"//& - " DOME - specified inflow on northern boundary\n"//& - " dyed_channel - supercritical with dye on the inflow boundary\n"//& - " dyed_obcs - circle_obcs with dyes on the open boundaries\n"//& - " Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//& - " shelfwave - Flather with shelf wave forcing on western boundary\n"//& - " supercritical - now only needed here for the allocations\n"//& - " tidal_bay - Flather with tidal forcing on eastern boundary\n"//& - " USER - user specified", default="none") - if (trim(config) == "DOME") then - call DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tracer_Reg) - elseif (trim(config) == "dyed_channel") then - call dyed_channel_set_OBC_tracer_data(OBC, G, GV, PF, tracer_Reg) - OBC%update_OBC = .true. - elseif (trim(config) == "dyed_obcs") then - call dyed_obcs_set_OBC_data(OBC, G, GV, PF, tracer_Reg) - elseif (trim(config) == "Kelvin") then - OBC%update_OBC = .true. - elseif (trim(config) == "shelfwave") then - OBC%update_OBC = .true. - elseif (lowercase(trim(config)) == "supercritical") then - call supercritical_set_OBC_data(OBC, G, GV, US, PF) - elseif (trim(config) == "tidal_bay") then - OBC%update_OBC = .true. - elseif (trim(config) == "USER") then - call user_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) - elseif (.not. trim(config) == "none") then - call MOM_error(FATAL, "The open boundary conditions specified by "//& - "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.") - endif + if (.not. associated(OBC)) return - if (debug) then - call hchksum(G%mask2dT, 'MOM_initialize_OBCs: mask2dT ', G%HI) - call uvchksum('MOM_initialize_OBCs: mask2dC[uv]', G%mask2dCu, G%mask2dCv, G%HI) - call qchksum(G%mask2dBu, 'MOM_initialize_OBCs: mask2dBu ', G%HI) - endif - if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h) + call callTree_enter('initialize_user_OBCs()') - if (OBC%use_h_res) & - call fill_thickness_segments(G, GV, US, OBC, h) + ! This controls user code for setting open boundary data + call get_param(PF, mdl, "OBC_USER_CONFIG", config, & + "A string that sets how the user code is invoked to set open boundary data: \n"//& + " \t DOME - specified inflow on northern boundary\n"//& + " \t dyed_channel - supercritical with dye on the inflow boundary\n"//& + " \t dyed_obcs - circle_obcs with dyes on the open boundaries\n"//& + " \t Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//& + " \t shelfwave - Flather with shelf wave forcing on western boundary\n"//& + " \t supercritical - now only needed here for the allocations\n"//& + " \t tidal_bay - Flather with tidal forcing on eastern boundary\n"//& + " \t USER - user specified", default="none") + if (trim(config) == "DOME") then + call DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tracer_Reg) + elseif (trim(config) == "dyed_channel") then + call dyed_channel_set_OBC_tracer_data(OBC, G, GV, PF, tracer_Reg) + OBC%update_OBC = .true. + elseif (trim(config) == "dyed_obcs") then + call dyed_obcs_set_OBC_data(OBC, G, GV, PF, tracer_Reg) + elseif (trim(config) == "Kelvin") then + OBC%update_OBC = .true. + elseif (trim(config) == "shelfwave") then + OBC%update_OBC = .true. + elseif (lowercase(trim(config)) == "supercritical") then + call supercritical_set_OBC_data(OBC, G, GV, US, PF) + elseif (trim(config) == "tidal_bay") then + OBC%update_OBC = .true. + elseif (trim(config) == "USER") then + call user_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) + elseif (.not. trim(config) == "none") then + call MOM_error(FATAL, "The open boundary conditions specified by "//& + "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.") endif - call callTree_leave('MOM_initialize_OBCs()') + call callTree_leave('initialize_user_OBCs()') -end subroutine MOM_initialize_OBCs +end subroutine initialize_user_OBCs !> Reads the layer thicknesses or interface heights from a file. subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, file_has_thickness, &