diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 31b7b1e889..4d24b9108e 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -742,7 +742,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%OBC)) then call radiation_open_bdry_conds(CS%OBC, u_av, u_old_rad_OBC, v_av, & - v_old_rad_OBC, hp, h_old_rad_OBC, G) + v_old_rad_OBC, hp, h_old_rad_OBC, G, dt) endif ! h_av = (h + hp)/2 @@ -1005,7 +1005,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%OBC)) then call radiation_open_bdry_conds(CS%OBC, u, u_old_rad_OBC, v, & - v_old_rad_OBC, h, h_old_rad_OBC, G) + v_old_rad_OBC, h, h_old_rad_OBC, G, dt) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index a46eb074d8..1a357ce1ea 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -639,7 +639,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (associated(CS%OBC)) then call radiation_open_bdry_conds(CS%OBC, u_av, u_old_rad_OBC, v_av, & - v_old_rad_OBC, hp, h_old_rad_OBC, G) + v_old_rad_OBC, hp, h_old_rad_OBC, G, dt) endif ! h_av = (h + hp)/2 @@ -855,7 +855,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (associated(CS%OBC)) then call radiation_open_bdry_conds(CS%OBC, u, u_old_rad_OBC, v, & - v_old_rad_OBC, h, h_old_rad_OBC, G) + v_old_rad_OBC, h, h_old_rad_OBC, G, dt) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3d2a687dd0..77ccf509b6 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -9,14 +9,19 @@ module MOM_open_boundary use MOM_domains, only : To_All, SCALAR_PAIR, CGRID_NE use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type, log_param -use MOM_grid, only : ocean_grid_type +use MOM_grid, only : ocean_grid_type, hor_index_type use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_io, only : EAST_FACE, NORTH_FACE -use MOM_io, only : slasher, read_data +use MOM_io, only : slasher, read_data, field_size use MOM_obsolete_params, only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char use MOM_string_functions, only : extract_word, remove_spaces use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type use MOM_variables, only : thermo_var_ptrs +use time_interp_external_mod, only : init_external_field, time_interp_external +use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS, initialize_remapping +use MOM_remapping, only : remapping_core_h, end_remapping +use MOM_regridding, only : regridding_CS +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -30,6 +35,8 @@ module MOM_open_boundary public open_boundary_impose_land_mask public radiation_open_bdry_conds public set_Flather_data +public update_obc_segment_data +public fill_OBC_halos integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 integer, parameter, public :: OBC_FLATHER = 3 @@ -38,9 +45,21 @@ module MOM_open_boundary integer, parameter, public :: OBC_DIRECTION_S = 200 !< Indicates the boundary is an effective southern boundary integer, parameter, public :: OBC_DIRECTION_E = 300 !< Indicates the boundary is an effective eastern boundary integer, parameter, public :: OBC_DIRECTION_W = 400 !< Indicates the boundary is an effective western boundary +integer, parameter :: MAX_OBC_FIELDS = 100 !< Maximum number of data fields needed for OBC segments + +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 + real, pointer, dimension(:,:,:) :: buffer_src=>NULL() ! 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, dimension(:,:,:), pointer :: dz_src=>NULL() ! vertical grid cell spacing of the incoming segment data (m) + real, dimension(:,:,:), pointer :: buffer_dst=>NULL() ! buffer src data remapped to the target vertical grid + real, dimension(:,:), pointer :: bt_vel=>NULL() ! barotropic velocity (m s-1) + real :: value ! constant value if fid is equal to -1 +end type OBC_segment_data_type -!> Open boundary segment type - we'll have one for each open segment -!! to describe that segment. type, public :: OBC_segment_type logical :: Flather !< If true, applies Flather + Chapman radiation of barotropic gravity waves. logical :: radiation !< If true, 1D Orlanksi radiation boundary conditions are applied. @@ -49,20 +68,38 @@ module MOM_open_boundary logical :: nudged !< Optional supplement to radiation boundary. logical :: specified !< Boundary fixed to external value. logical :: gradient !< Zero gradient at boundary. - logical :: values_needed !< Whether or not baroclinic OBC fields are needed. + logical :: values_needed !< Whether or not external OBC fields are needed. logical :: legacy !< Old code for tangential BT velocities. integer :: direction !< Boundary faces one of the four directions. + type(OBC_segment_data_type), pointer, dimension(:) :: field=>NULL() !< OBC data + integer :: num_fields !< number of OBC data fields (e.g. u_normal,u_parallel and eta for Flather) + character(len=32), pointer, dimension(:) :: field_names=>NULL() !< field names for this segment integer :: Is_obc !< i-indices of boundary segment. integer :: Ie_obc !< i-indices of boundary segment. integer :: Js_obc !< j-indices of boundary segment. integer :: Je_obc !< j-indices of boundary segment. real :: Tnudge_in !< Inverse nudging timescale on inflow (1/s). real :: Tnudge_out !< Inverse nudging timescale on outflow (1/s). + logical :: on_pe !< true if segment is located in the computational domain + real, pointer, dimension(:,:) :: Cg=>NULL() !NULL() !NULL() !NULL() !NULL() !NULL() !NULL() !NULL() !NULL() !NULL()!NULL() ! Open-boundary data type, public :: ocean_OBC_type - integer :: number_of_segments = 0 !< The number of open-boundary segments. + integer :: number_of_segments = 0 !< The number of open-boundary segments. + integer :: ke = 0 !< The number of model layers logical :: Flather_u_BCs_exist_globally = .false. !< True if any zonal velocity points !! in the global domain use Flather BCs. logical :: Flather_v_BCs_exist_globally = .false. !< True if any meridional velocity points @@ -75,6 +112,9 @@ module MOM_open_boundary !! in the global domain use specified BCs. logical :: specified_v_BCs_exist_globally = .false. !< True if any meridional velocity points !! in the global domain use specified BCs. + logical :: user_BCs_set_globally = .false. !< True if any OBC_CONFIG or OBC_VALUES_CONFIG + !! set for input from user directory. + logical :: update_OBC = .false. !< Is the open boundary info going to get updated? real :: g_Earth ! Properties of the segments used. type(OBC_segment_type), pointer, dimension(:) :: & @@ -121,8 +161,8 @@ module MOM_open_boundary !! velocity (or speed of characteristics), in m s-1. The !! default value is 10 m s-1. logical :: OBC_pe !< Is there an open boundary on this tile? - logical :: update_OBC = .false. !< Is the open boundary info going to get updated? character(len=200) :: OBC_values_config + type(remapping_CS), pointer :: remap_CS ! ALE remapping control structure for segments only end type ocean_OBC_type integer :: id_clock_pass @@ -134,6 +174,12 @@ module MOM_open_boundary contains !> Enables OBC module and reads configuration parameters +!> This routine is called from MOM_initialize_fixed which +!> occurs before the initialization of the vertical coordinate +!> and ALE_init. Therefore segment data are not fully initialized +!> here. The remainder of the segment data are initialized in a +!> later call to update_open_boundary_data + subroutine open_boundary_config(G, param_file, OBC) type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle @@ -141,7 +187,9 @@ subroutine open_boundary_config(G, param_file, OBC) ! Local variables integer :: l ! For looping over segments character(len=15) :: segment_param_str ! The run-time parameter name for each segment - character(len=100) :: segment_str ! The contents (rhs) for parameter "segment_param_str" + character(len=100) :: segment_str ! The contents (rhs) for parameter "segment_param_str" + character(len=200) :: config1 ! String for OBC_CONFIG + character(len=200) :: config2 ! String for OBC_VALUES_CONFIG allocate(OBC) @@ -153,6 +201,17 @@ subroutine open_boundary_config(G, param_file, OBC) call get_param(param_file, mod, "G_EARTH", OBC%g_Earth, & "The gravitational acceleration of the Earth.", & units="m s-2", default = 9.80) + call get_param(param_file, mod, "OBC_CONFIG", config1, & + "A string that sets how the open boundary conditions are \n"//& + " configured: \n", default="None", do_not_log=.true.) + call get_param(param_file, mod, "OBC_VALUES_CONFIG", config2, & + "A string that sets how the open boundary values are \n"//& + " configured: \n", default="None", do_not_log=.true.) + call get_param(param_file, mod, "NZ", OBC%ke, & + "The number of model layers", default=0, do_not_log=.true.) + + if (config1 .ne. "None" .or. config2 .ne. "None") OBC%user_BCs_set_globally = .true. + if (OBC%number_of_segments > 0) then ! Allocate everything allocate(OBC%OBC_segment_number(0:OBC%number_of_segments)) @@ -187,6 +246,10 @@ subroutine open_boundary_config(G, param_file, OBC) "Unable to interpret "//segment_param_str//" = "//trim(segment_str)) endif enddo + +! if (open_boundary_query(OBC, needs_ext_seg_data=.true.) .and. .not. OBC%user_BCs_set_globally) & +! if (open_boundary_query(OBC, needs_ext_seg_data=.true.)) & +! call initialize_segment_data(G, OBC, param_file) endif ! Safety check @@ -203,6 +266,200 @@ subroutine open_boundary_config(G, param_file, OBC) end subroutine open_boundary_config +subroutine initialize_segment_data(G, OBC, PF) + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure + type(param_file_type), intent(in) :: PF !< Parameter file handle + + integer :: num_segs,n,m,num_fields + character(len=256) :: segstr, filename + character(len=20) :: segnam,suffix + character(len=32) :: varnam, fieldname + real :: value + integer :: orient + character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names + character(len=128) :: inputdir + type(OBC_segment_type), pointer, dimension(:) :: OBC_segments ! pointer to segment type list + character(len=32) :: remappingScheme + logical :: check_reconstruction, check_remapping, force_bounds_in_subcell + integer, dimension(4) :: siz,siz2 + integer :: is_obc,ie_obc,js_obc,je_obc + + num_segs = OBC%number_of_segments + + ! There is a problem with the order of the OBC initialization + ! with respect to ALE_init. Currently handling this by copying the + ! param file so that I can use it later in step_MOM in order to finish + ! initializing segments on the first step. + + call get_param(PF, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + + call get_param(PF, mod, "REMAPPING_SCHEME", remappingScheme, & + "This sets the reconstruction scheme used\n"//& + "for vertical remapping for all variables.\n"//& + "It can be one of the following schemes:\n"//& + trim(remappingSchemesDoc), default=remappingDefaultScheme) + call get_param(PF, mod, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & + "If true, cell-by-cell reconstructions are checked for\n"//& + "consistency and if non-monotonicty or an inconsistency is\n"//& + "detected then a FATAL error is issued.", default=.false.) + call get_param(PF, mod, "FATAL_CHECK_REMAPPING", check_remapping, & + "If true, the results of remapping are checked for\n"//& + "conservation and new extrema and if an inconsistency is\n"//& + "detected then a FATAL error is issued.", default=.false.) + call get_param(PF, mod, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, & + "If true, the values on the intermediate grid used for remapping\n"//& + "are forced to be bounded, which might not be the case due to\n"//& + "round off.", default=.false.) + + allocate(OBC%remap_CS) + call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & + check_reconstruction=check_reconstruction, & + check_remapping=check_remapping, force_bounds_in_subcell=force_bounds_in_subcell) + + OBC_segments => OBC%OBC_segment_number(1:num_segs) + + do n=1, num_segs + write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n + write(suffix,"('_segment_',i3.3)") n + call get_param(PF, mod, segnam, segstr) + + call parse_segment_data_str(trim(segstr),fields=fields, num_fields=num_fields) + if (num_fields == 0) cycle ! cycle to next segment + allocate(OBC_segments(n)%field(num_fields)) + + if (OBC_segments(n)%Flather) then + if (num_fields /= 3) call MOM_error(FATAL, & + "MOM_open_boundary, initialize_segment_data: "//& + "Need three inputs for Flather") + OBC_segments(n)%num_fields = 3 ! these are the input fields required for the Flather option + ! note that this is assuming that the inputs are coming in this order + ! and independent of the input param string . Needs cleanup - mjh + allocate(OBC_segments(n)%field_names(OBC_segments(n)%num_fields)) + OBC_segments(n)%field_names(:)='None' + OBC_segments(n)%field_names(1)='UO' + OBC_segments(n)%field_names(2)='VO' + OBC_segments(n)%field_names(3)='ZOS' + endif + +!! +! CODE HERE FOR OTHER OPTIONS (CLAMPED, NUDGED,..) +!! + + do m=1,num_fields + call parse_segment_data_str(segstr,var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname) + OBC_segments(n)%field(m)%name = trim(fields(m)) + if (trim(filename) /= 'none') then + filename = trim(inputdir)//trim(filename) + fieldname = trim(fieldname)//trim(suffix) + call field_size(filename,fieldname,siz,no_domain=.true.) + if (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0) then + call MOM_error(FATAL,'segment data are not on the supergrid') + endif + + siz2(1)=1 + if (siz(1)>1) then + siz2(1)=(siz(1)-1)/2 + endif + siz2(2)=1 + if (siz(2)>1) then + siz2(2)=(siz(2)-1)/2 + endif + siz2(3)=siz(3) + + allocate(OBC_segments(n)%field(m)%buffer_src(is_obc:ie_obc,js_obc:je_obc,siz2(3))) + OBC_segments(n)%field(m)%buffer_src(:,:,:)=0.0 + OBC_segments(n)%field(m)%fid = init_external_field(trim(filename),trim(fieldname)) + if (siz(3) > 1) then + fieldname = 'dz_'//trim(fieldname) + call field_size(filename,fieldname,siz,no_domain=.true.) + if (OBC_segments(n)%direction == OBC_DIRECTION_E .or. OBC_segments(n)%direction == OBC_DIRECTION_W) then + allocate(OBC_segments(n)%field(m)%dz_src(is_obc:ie_obc,js_obc+1:je_obc,siz(3))) + else + allocate(OBC_segments(n)%field(m)%dz_src(is_obc+1:ie_obc,js_obc:je_obc,siz(3))) + endif + OBC_segments(n)%field(m)%dz_src(:,:,:)=0.0 + OBC_segments(n)%field(m)%nk_src=siz(3) + OBC_segments(n)%field(m)%fid_dz = init_external_field(trim(filename),trim(fieldname)) + else + OBC_segments(n)%field(m)%nk_src=1 + endif + else + OBC_segments(n)%field(m)%fid = -1 + OBC_segments(n)%field(m)%value = value + endif + enddo + + enddo + + return + +end subroutine initialize_segment_data + +!< Define indices for segment and store in hor_index_type +!< using global segment bounds corresponding to q-points +subroutine setup_segment_indices(G, seg, Is_obc, Ie_obc, Js_obc, Je_obc) + type(dyn_horgrid_type), intent(in) :: G !< grid type + type(OBC_segment_type), intent(inout) :: seg !< Open boundary segment + integer, intent(in) :: Is_obc !< Q-point global i-index of start of segment + integer, intent(in) :: Ie_obc !< Q-point global i-index of end of segment + integer, intent(in) :: Js_obc !< Q-point global j-index of start of segment + integer, intent(in) :: Je_obc !< Q-point global j-index of end of segment + ! Local variables + integer :: Isg,Ieg,Jsg,Jeg + +! if (.not. G%Domain%symmetric) call MOM_error(FATAL, "MOM_open_boundary.F90, setup_segment_indices: "//& +! "Need to compile in symmetric mode") + + ! Isg, Ieg will be I*_obc in global space + if (Ie_obc Parse an OBC_SEGMENT_%%% string starting with "I=" and configure placement and type of OBC accordingly subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure @@ -216,9 +473,14 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) ! This returns the global indices for the segment call parse_segment_str(G%ieg, G%jeg, segment_str, I_obc, Js_obc, Je_obc, action_str ) + + call setup_segment_indices(G, OBC%OBC_segment_number(l_seg),I_obc,I_obc,Js_obc,Je_obc) + call allocate_OBC_segment_data(OBC, OBC%OBC_segment_number(l_seg)) + I_obc = I_obc - G%idg_offset ! Convert to local tile indices on this tile Js_obc = Js_obc - G%jdg_offset ! Convert to local tile indices on this tile Je_obc = Je_obc - G%jdg_offset ! Convert to local tile indices on this tile + this_kind = OBC_NONE ! Hack to extend segment by one point @@ -228,14 +490,23 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) Js_obc = Js_obc + 1 ; Je_obc = Je_obc - 1 endif - if (Je_obc>Js_obc) OBC%OBC_segment_number(l_seg)%direction = OBC_DIRECTION_E - if (Je_obcJs_obc) then + OBC%OBC_segment_number(l_seg)%direction = OBC_DIRECTION_E + else if (Je_obcG%HI%IedB) return ! Boundary is not on tile - if (max(Js_obc,Je_obc)G%HI%JedB) return ! Segment is not on tile +! if (max(Js_obc,Je_obc)G%HI%JedB) return ! Segment is not on tile + if (Js_obcG%HI%JedB) return ! Segment is not on tile + + OBC%OBC_segment_number(l_seg)%on_pe = .true. do j=G%HI%jsd, G%HI%jed - if (j>min(Js_obc,Je_obc) .and. j<=max(Js_obc,Je_obc)) then + if (j>Js_obc .and. j<=Je_obc) then + OBC%OBC_segment_u(I_obc,j) = l_seg - if (Je_obc>Js_obc) then ! East is outward + if (OBC%OBC_segment_number(l_seg)%direction == OBC_DIRECTION_E) then ! East is outward if (this_kind == OBC_FLATHER) then ! Set v points outside segment if (OBC%OBC_segment_v(i_obc+1,J) == OBC_NONE) then @@ -324,6 +599,10 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg) ! This returns the global indices for the segment call parse_segment_str(G%ieg, G%jeg, segment_str, J_obc, Is_obc, Ie_obc, action_str ) + + call setup_segment_indices(G, OBC%OBC_segment_number(l_seg),Is_obc,Ie_obc,J_obc,J_obc) + call allocate_OBC_segment_data(OBC, OBC%OBC_segment_number(l_seg)) + J_obc = J_obc - G%jdg_offset ! Convert to local tile indices on this tile Is_obc = Is_obc - G%idg_offset ! Convert to local tile indices on this tile Ie_obc = Ie_obc - G%idg_offset ! Convert to local tile indices on this tile @@ -336,14 +615,25 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg) Is_obc = Is_obc + 1 ; Ie_obc = Ie_obc - 1 endif - if (Ie_obc>Is_obc) OBC%OBC_segment_number(l_seg)%direction = OBC_DIRECTION_S - if (Ie_obcIs_obc) then + OBC%OBC_segment_number(l_seg)%direction = OBC_DIRECTION_S + else if (Ie_obcIs_obc) OBC%OBC_segment_number(l_seg)%direction = OBC_DIRECTION_S +! if (Ie_obcG%HI%JedB) return ! Boundary is not on tile - if (max(Is_obc,Ie_obc)G%HI%IedB) return ! Segment is not on tile + if (Is_obcG%HI%IedB) return ! Segment is not on tile + + OBC%OBC_segment_number(l_seg)%on_pe = .true. + +! if (J_obcG%HI%JedB) return ! Boundary is not on tile +! if (max(Is_obc,Ie_obc)G%HI%IedB) return ! Segment is not on tile do i=G%HI%isd, G%HI%ied - if (i>min(Is_obc,Ie_obc) .and. i<=max(Is_obc,Ie_obc)) then + if (i>Is_obc .and. i<=Ie_obc) then OBC%OBC_segment_v(i,J_obc) = l_seg - if (Is_obc>Ie_obc) then ! North is outward + if (OBC%OBC_segment_number(l_seg)%direction == OBC_DIRECTION_N) then ! North is outward if (this_kind == OBC_FLATHER) then ! Set u points outside segment if (OBC%OBC_segment_u(I,j_obc+1) == OBC_NONE) then @@ -511,6 +806,89 @@ integer function interpret_int_expr(string, imax) end function interpret_int_expr end subroutine parse_segment_str + subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fields, num_fields, debug ) + character(len=*), intent(in) :: segment_str !< A string in form of "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..." + character(len=*), intent(in), optional :: var !< The name of the variable for which parameters are needed + character(len=*), intent(out), optional :: filenam !< The name of the input file if using "file" method + character(len=*), intent(out), optional :: fieldnam !< The name of the variable in the input file if using "file" method + real, intent(out), optional :: value !< A constant value if using the "value" method + character(len=*), dimension(MAX_OBC_FIELDS), intent(out), optional :: fields !< List of fieldnames for each segment + integer, intent(out), optional :: num_fields + logical, intent(in), optional :: debug + ! Local variables + character(len=128) :: word1, word2, word3, method + integer :: lword, nfields, n, m, orient + logical :: continue,dbg + character(len=32), dimension(MAX_OBC_FIELDS) :: flds + + nfields=0 + continue=.true. + dbg=.false. + if (PRESENT(debug)) dbg=debug + + do while (continue) + word1 = extract_word(segment_str,',',nfields+1) + if (trim(word1) == '') exit + nfields=nfields+1 + word2 = extract_word(word1,'=',1) + flds(nfields) = trim(word2) + enddo + + if (PRESENT(fields)) then + do n=1,nfields + fields(n) = flds(n) + enddo + endif + + if (PRESENT(num_fields)) then + num_fields=nfields + return + endif + + m=0 + if (PRESENT(var)) then + do n=1,nfields + if (trim(var)==trim(flds(n))) then + m=n + exit + endif + enddo + if (m==0) then + call abort() + endif + + ! Process first word which will start with the fieldname + word3 = extract_word(segment_str,',',m) + word1 = extract_word(word3,':',1) +! if (trim(word1) == '') exit + word2 = extract_word(word1,'=',1) + if (trim(word2) == trim(var)) then + method=trim(extract_word(word1,'=',2)) + lword=len_trim(method) + if (method(lword-3:lword) == 'file') then + ! raise an error id filename/fieldname not in argument list + word1 = extract_word(word3,':',2) + filenam = extract_word(word1,'(',1) + fieldnam = extract_word(word1,'(',2) + lword=len_trim(fieldnam) + fieldnam = fieldnam(1:lword-1) ! remove trailing parenth + value=-999. + elseif (method(lword-4:lword) == 'value') then + filenam = 'none' + fieldnam = 'none' + word1 = extract_word(word3,':',2) + lword=len_trim(word1) + read(word1(1:lword),*,end=986,err=987) value + endif + endif + endif + + return + 986 call MOM_error(FATAL,'End of record while parsing segment data specification! '//trim(segment_str)) + 987 call MOM_error(FATAL,'Error while parsing segment data specification! '//trim(segment_str)) + + end subroutine parse_segment_data_str + !> Initialize open boundary control structure subroutine open_boundary_init(G, param_file, OBC) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure @@ -546,12 +924,12 @@ subroutine open_boundary_init(G, param_file, OBC) end subroutine open_boundary_init -!> Query the state of open boundary module configuration -logical function open_boundary_query(OBC, apply_specified_OBC, apply_Flather_OBC, apply_nudged_OBC) +logical function open_boundary_query(OBC, apply_specified_OBC, apply_Flather_OBC, apply_nudged_OBC, needs_ext_seg_data) type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure logical, optional, intent(in) :: apply_specified_OBC !< If present, returns True if specified_*_BCs_exist_globally is true logical, optional, intent(in) :: apply_Flather_OBC !< If present, returns True if Flather_*_BCs_exist_globally is true logical, optional, intent(in) :: apply_nudged_OBC !< If present, returns True if nudged_*_BCs_exist_globally is true + logical, optional, intent(in) :: needs_ext_seg_data !< If present, returns True if external segment data needed open_boundary_query = .false. if (.not. associated(OBC)) return if (present(apply_specified_OBC)) open_boundary_query = OBC%specified_u_BCs_exist_globally .or. & @@ -560,6 +938,9 @@ logical function open_boundary_query(OBC, apply_specified_OBC, apply_Flather_OBC OBC%Flather_v_BCs_exist_globally if (present(apply_nudged_OBC)) open_boundary_query = OBC%nudged_u_BCs_exist_globally .or. & OBC%nudged_v_BCs_exist_globally + if (present(needs_ext_seg_data) .and. OBC%number_of_segments > 0 .and. & + OBC%update_OBC) open_boundary_query = OBC%update_OBC + end function open_boundary_query !> Deallocate open boundary data @@ -722,15 +1103,16 @@ end subroutine open_boundary_impose_land_mask !> Apply radiation conditions to 3D u,v (,h) at open boundaries subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, & - h_new, h_old, G) + h_new, h_old, G, dt) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< New u values on open boundaries + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< New u values on open boundaries real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old !< Original unadjusted u real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new !< New v values on open boundaries real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v_old !< Original unadjusted v real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_new !< New h values on open boundaries real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old !< Original h values + real, intent(in) :: dt !< Appropriate timestep ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: grad real :: dhdt, dhdx, dhdy, gamma_u, gamma_h, gamma_v @@ -795,7 +1177,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, & else tau = OBC%OBC_segment_number(OBC%OBC_segment_u(I,j))%Tnudge_out endif - u_new(I,j,k) = u_new(I,j,k) + tau*(OBC%u(I,j,k) - u_old(I,j,k)) + u_new(I,j,k) = u_new(I,j,k) + dt*tau*(OBC%u(I,j,k) - u_old(I,j,k)) endif endif @@ -844,7 +1226,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, & else tau = OBC%OBC_segment_number(OBC%OBC_segment_u(I,j))%Tnudge_out endif - u_new(I,j,k) = u_new(I,j,k) + tau*(OBC%u(I,j,k) - u_old(I,j,k)) + u_new(I,j,k) = u_new(I,j,k) + dt*tau*(OBC%u(I,j,k) - u_old(I,j,k)) endif endif @@ -949,7 +1331,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, & else tau = OBC%OBC_segment_number(OBC%OBC_segment_v(i,J))%Tnudge_out endif - v_new(i,J,k) = v_new(i,J,k) + tau*(OBC%v(i,J,k) - v_old(i,J,k)) + v_new(i,J,k) = v_new(i,J,k) + dt*tau*(OBC%v(i,J,k) - v_old(i,J,k)) endif endif @@ -998,7 +1380,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, & else tau = OBC%OBC_segment_number(OBC%OBC_segment_v(i,J))%Tnudge_out endif - v_new(i,J,k) = v_new(i,J,k) + tau*(OBC%v(i,J,k) - v_old(i,J,k)) + v_new(i,J,k) = v_new(i,J,k) + dt*tau*(OBC%v(i,J,k) - v_old(i,J,k)) endif endif @@ -1163,7 +1545,6 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) ! allocate(OBC%ry_old_h(isd:ied,Jsd:Jed,nz)) ; OBC%ry_old_h(:,:,:) = 0.0 endif - if (associated(tv%T)) then allocate(OBC_T_u(IsdB:IedB,jsd:jed,nz)) ; OBC_T_u(:,:,:) = 0.0 allocate(OBC_S_u(IsdB:IedB,jsd:jed,nz)) ; OBC_S_u(:,:,:) = 0.0 @@ -1298,6 +1679,410 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) end subroutine set_Flather_data +function lookup_seg_field(OBC_seg,field) + type(OBC_segment_type), pointer :: OBC_seg + character(len=32), intent(in) :: field ! The field name + integer :: lookup_seg_field + + integer :: n,m + + lookup_seg_field=-1 + do n=1,OBC_seg%num_fields + if (trim(field) == OBC_seg%field_names(n)) then + lookup_seg_field=n + return + endif + enddo + + return + +end function lookup_seg_field + + +subroutine fill_OBC_halos(G, GV, OBC, tv, h, tracer_Reg) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Thickness + type(tracer_registry_type), pointer :: tracer_Reg !< Tracer registry + ! Local variables + + integer :: i, j, k, isd, ied, jsd, jed, is, ie, js, je + integer :: n, nz + character(len=40) :: mod = "fill_OBC_halos" ! This subroutine's name. + type(OBC_segment_type), pointer :: segment + real, pointer, dimension(:,:,:) :: & + OBC_T_u => NULL(), & ! These arrays should be allocated and set to + OBC_T_v => NULL(), & ! specify the values of T and S that should come + OBC_S_u => NULL(), & + OBC_S_v => NULL() + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + nz=G%ke + + if (.not. associated(OBC) .or. .not. associated(tv%T)) return + + call pass_var(tv%T,G%Domain) + call pass_var(tv%S,G%Domain) + + do n = 1, OBC%number_of_segments + segment => OBC%OBC_segment_number(n) + + if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain + + do j=G%jsdB+1,G%jedB ; do I=G%isdB,G%iedB-1 + if (segment%direction == OBC_DIRECTION_E .and. OBC%OBC_segment_u(I,j) /= OBC_NONE ) then + do k=1,nz + tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k); h(i+1,j,k) = h(i,j,k) + enddo + endif + enddo ; enddo + + do j=G%jsdB+1,G%jedB ; do I=G%isdB+1,G%iedB + if (segment%direction == OBC_DIRECTION_W .and. OBC%OBC_segment_u(I,j) /= OBC_NONE ) then + do k=1,nz + tv%T(i-1,j,k) = tv%T(i,j,k) ; tv%S(i-1,j,k) = tv%S(i,j,k); h(i-1,j,k) = h(i,j,k) + enddo + endif + enddo ; enddo + + do j=G%jsdB,G%jedB-1 ; do I=G%isdB+1,G%iedB + if (segment%direction == OBC_DIRECTION_N .and. OBC%OBC_segment_v(I,j) /= OBC_NONE ) then + do k=1,nz + tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k); h(i,j+1,k) = h(i,j,k) + enddo + endif + enddo ; enddo + + do j=G%jsdB+1,G%jedB ; do I=G%isdB+1,G%iedB + if (segment%direction == OBC_DIRECTION_S .and. OBC%OBC_segment_v(I,j) /= OBC_NONE ) then + do k=1,nz + tv%T(i,j-1,k) = tv%T(i,j,k) ; tv%S(i,j-1,k) = tv%S(i,j,k); h(i,j-1,k) = h(i,j,k) + enddo + endif + enddo ; enddo + + enddo + + if (associated(tv%T)) then + if (.not. associated(OBC_T_u)) then + allocate(OBC_T_u(G%IsdB:G%IedB,G%jsd:G%jed,nz)) ; OBC_T_u(:,:,:) = 0.0 + allocate(OBC_S_u(G%IsdB:G%IedB,G%jsd:G%jed,nz)) ; OBC_S_u(:,:,:) = 0.0 + allocate(OBC_T_v(G%Isd:G%Ied,G%jsdB:G%jedB,nz)) ; OBC_T_v(:,:,:) = 0.0 + allocate(OBC_S_v(G%Isd:G%Ied,G%jsdB:G%jedB,nz)) ; OBC_S_v(:,:,:) = 0.0 + endif + do k=1,nz ; do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 + if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then + if (OBC%OBC_segment_number(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_E) then + OBC_T_u(I,j,k) = tv%T(i,j,k) + OBC_S_u(I,j,k) = tv%S(i,j,k) + elseif (OBC%OBC_segment_number(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_W) then + OBC_T_u(I,j,k) = tv%T(i+1,j,k) + OBC_S_u(I,j,k) = tv%S(i+1,j,k) + elseif (G%mask2dT(i,j) + G%mask2dT(i+1,j) > 0) then + OBC_T_u(I,j,k) = (G%mask2dT(i,j)*tv%T(i,j,k) + G%mask2dT(i+1,j)*tv%T(i+1,j,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + OBC_S_u(I,j,k) = (G%mask2dT(i,j)*tv%S(i,j,k) + G%mask2dT(i+1,j)*tv%S(i+1,j,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + else ! This probably shouldn't happen or maybe it doesn't matter? + OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k)) + OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k)) + endif + else + OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k)) + OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k)) + endif + enddo; enddo ; enddo + + do k=1,nz ; do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied + if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then + if (OBC%OBC_segment_number(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_N) then + OBC_T_v(i,J,k) = tv%T(i,j,k) + OBC_S_v(i,J,k) = tv%S(i,j,k) + elseif (OBC%OBC_segment_number(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_S) then + OBC_T_v(i,J,k) = tv%T(i,j+1,k) + OBC_S_v(i,J,k) = tv%S(i,j+1,k) + elseif (G%mask2dT(i,j) + G%mask2dT(i,j+1) > 0) then + OBC_T_v(i,J,k) = (G%mask2dT(i,j)*tv%T(i,j,k) + G%mask2dT(i,j+1)*tv%T(i,j+1,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + OBC_S_v(i,J,k) = (G%mask2dT(i,j)*tv%S(i,j,k) + G%mask2dT(i,j+1)*tv%S(i,j+1,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + else ! This probably shouldn't happen or maybe it doesn't matter? + OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k)) + OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k)) + endif + else + OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k)) + OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k)) + endif + enddo; enddo ; enddo + + call pass_vector(OBC_T_u, OBC_T_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + call pass_vector(OBC_S_u, OBC_S_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + + call add_tracer_OBC_values("T",tracer_Reg, OBC_in_u=OBC_T_u,OBC_in_v=OBC_T_v) + call add_tracer_OBC_values("S",tracer_Reg, OBC_in_u=OBC_S_u,OBC_in_v=OBC_S_v) + + + endif +end subroutine fill_OBC_halos + +!> Allocate segment data fields +subroutine allocate_OBC_segment_data(OBC, segment) + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(OBC_segment_type), intent(inout) :: segment !< Open boundary segment + ! Local variables + integer :: isd, ied, jsd, jed + integer :: IsdB, IedB, JsdB, JedB + character(len=40) :: mod = "allocate_OBC_segment_data" ! This subroutine's name. + + 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 + + if (.not. segment%on_pe) return + + if (.not. ASSOCIATED(segment%Cg)) then ! finishing allocating storage for segments + allocate(segment%Cg(isd:ied,jsd:jed)); segment%Cg(:,:)=0. + allocate(segment%Htot(isd:ied,jsd:jed)); segment%Htot(:,:)=0.0 + allocate(segment%h(isd:ied,jsd:jed,OBC%ke)); segment%h(:,:,:)=0.0 + allocate(segment%normal_vel(isd:ied,jsd:jed,OBC%ke)); segment%normal_vel(:,:,:)=0.0 + allocate(segment%normal_trans(isd:ied,jsd:jed,OBC%ke)); segment%normal_trans(:,:,:)=0.0 + allocate(segment%normal_vel_bt(isd:ied,jsd:jed)); segment%normal_vel_bt(:,:)=0.0 + allocate(segment%normal_trans_bt(isd:ied,jsd:jed)); segment%normal_trans_bt(:,:)=0.0 + allocate(segment%tangent_vel(IsdB:IedB,JsdB:JedB,OBC%ke));segment%tangent_vel(:,:,:)=0.0 + allocate(segment%tangent_vel_bt(IsdB:IedB,JsdB:JedB)); segment%tangent_vel_bt(:,:)=0.0 + allocate(segment%eta(isd:ied,jsd:jed)); segment%eta(:,:)=0.0 + endif +end subroutine allocate_OBC_segment_data + +subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Thickness +! real, dimension(SZI_(G),SZJ_(G)) , intent(inout) :: eta !< Thickness + type(time_type), intent(in) :: Time + ! Local variables + + integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed + integer :: IsdB, IedB, JsdB, JedB, n, m, nz + character(len=40) :: mod = "set_OBC_segment_data" ! This subroutine's name. + character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path + type(OBC_segment_type), pointer :: segment + integer, dimension(4) :: siz,siz2 + real :: sumh ! column sum of thicknesses (m) + integer :: ni_seg, nj_seg ! number of src gridpoints along the segments + integer :: i2, j2 ! indices for referencing local domain array + integer :: is_obc, ie_obc, js_obc, je_obc ! segment indices within local domain + integer :: ishift, jshift ! offsets for staggered locations + real, dimension(:,:), pointer :: seg_vel => NULL() ! pointer to segment velocity array + real, dimension(:,:), pointer :: seg_trans => NULL() ! pointer to segment transport array + real, dimension(:,:,:), allocatable :: tmp_buffer + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + nz=G%ke + + if (.not. associated(OBC)) return + + do n = 1, OBC%number_of_segments + segment => OBC%OBC_segment_number(n) + + if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain + + ni_seg = segment%ie_obc-segment%is_obc+1 + nj_seg = segment%je_obc-segment%js_obc+1 + is_obc = max(segment%is_obc,is-1) + ie_obc = min(segment%ie_obc,ie) + js_obc = max(segment%js_obc,js-1) + je_obc = min(segment%je_obc,je) + + + if (segment%direction == OBC_DIRECTION_W .or. segment%direction == OBC_DIRECTION_E) then + nj_seg=nj_seg-1 + js_obc=js_obc+1 + else + ni_seg=ni_seg-1 + is_obc=is_obc+1 + endif + +! do j=jsd,jed ; do I=isd,ied-1 +! if (segment%direction == OBC_DIRECTION_E .and. OBC%OBC_segment_u(I,j) /= OBC_NONE ) then +! do k=1,nz +! tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k); h(i+1,j,k) = h(i,j,k) +! enddo +! else if (segment%direction == OBC_DIRECTION_W .and. OBC%OBC_segment_u(I,j) /= OBC_NONE ) then +! tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k); h(i,j,k) = h(i+1,j,k) +! endif +! enddo ; enddo + +! do j=jsd,jed-1 ; do I=isd,ied-1 +! if (segment%direction == OBC_DIRECTION_N .and. OBC%OBC_segment_v(I,j) /= OBC_NONE ) then +! do k=1,nz +! tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k); h(i,j+1,k) = h(i,j,k) +! enddo +! else if (segment%direction == OBC_DIRECTION_S .and. OBC%OBC_segment_v(I,j) /= OBC_NONE ) then +! tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k); h(i,j,k) = h(i,j+1,k) +! endif +! enddo ; enddo + +! Calculate auxiliary fields at staggered locations. +! Segment indices are on q points: +! +! |-----------|------------|-----------|-----------| J_obc +! Is_obc Ie_obc +! +! i2 has to start at Is_obc+1 and end at Ie_obc. +! j2 is J_obc and jshift has to be +1 at both the north and south. + + ! calculate auxiliary fields at staggered locations + ishift=0;jshift=0 + if (segment%direction == OBC_DIRECTION_W .or. segment%direction == OBC_DIRECTION_E) then + if (segment%direction == OBC_DIRECTION_W) ishift=1 + do j=js_obc,je_obc + do I=Is_obc,Ie_obc +! i2 = segment%Is_obc + i - 1 +! j2 = segment%Js_obc + j - 1 +! if ((i2 .gt. ied .or. i2 .lt. isd) .or. (j2 .gt. jed .or. j2 .lt. jsd)) cycle +! if (OBC%OBC_segment_u(i2,j2) /= n) cycle + segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j)) + if (GV%Boussinesq) then + segment%Htot(I,j) = G%bathyT(i+ishift,j)*GV%m_to_H! + eta(i+ishift,j) +! else +! segment%Htot(I,j) = eta(i+ishift,j) + endif + do k=1,G%ke + segment%h(I,j,k) = h(i+ishift,j,k) + enddo + enddo + enddo + + else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S) + if (segment%direction == OBC_DIRECTION_S) jshift=1 + + do J=Js_obc,Je_obc + do i=is_obc,ie_obc + segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift)) + if (GV%Boussinesq) then + segment%Htot(i,J) = G%bathyT(i,j+jshift)*GV%m_to_H! + eta(i,j+jshift) +! else +! segment%Htot(i,J) = eta(i,j+jshift) + endif + do k=1,G%ke + segment%h(i,J,k) = h(i,j+jshift,k) + enddo + enddo + enddo + endif + + do m = 1,segment%num_fields + 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) + siz(3)=size(segment%field(m)%buffer_src,3) + if (.not.associated(segment%field(m)%buffer_dst)) then + if (siz(3) /= segment%field(m)%nk_src) call MOM_error(FATAL,'nk_src inconsistency') + if (segment%field(m)%nk_src > 1) then + allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) + else + allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1)) + endif + segment%field(m)%buffer_dst(:,:,:)=0.0 + if (trim(segment%field(m)%name) == 'UO' .or. trim(segment%field(m)%name) == 'VO') then + allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc)) + segment%field(m)%bt_vel(:,:)=0.0 + endif + endif + ! read source data interpolated to the current model time + if (siz(1)==1) then + allocate(tmp_buffer(1,nj_seg*2+1,segment%field(m)%nk_src)) ! segment data is currrently on supergrid + else + allocate(tmp_buffer(ni_seg*2+1,1,segment%field(m)%nk_src)) ! segment data is currrently on supergrid + endif + call time_interp_external(segment%field(m)%fid,Time, tmp_buffer) + if (siz(1)==1) then + segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,2*(js_obc+G%jdg_offset)-1:2*(je_obc+G%jdg_offset)-1:2,:) + else + segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(2*(is_obc+G%idg_offset)-1:2*(ie_obc+G%idg_offset)-1:2,1,:) + endif + if (segment%field(m)%nk_src > 1) then + call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer) + if (siz(1)==1) then + segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,2*(js_obc+G%jdg_offset)-1:2*(je_obc+G%jdg_offset)-1:2,:) + else + segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(2*(is_obc+G%idg_offset)-1:2*(ie_obc+G%idg_offset)-1:2,1,:) + endif + do j=js_obc,je_obc + do i=is_obc,ie_obc + + ! Using the h remapping approach + ! Pretty sure we need to check for source/target grid consistency here + segment%field(m)%buffer_dst(i,j,:)=0.0 ! initialize remap destination buffer + if (G%mask2dT(i,j)>0.) then + call remapping_core_h(segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:),& + segment%field(m)%buffer_src(i,j,:),G%ke, h(i,j,:),& + segment%field(m)%buffer_dst(i,j,:),OBC%remap_CS) + endif + enddo + enddo + else ! 2d data + segment%field(m)%buffer_dst(:,:,1)=segment%field(m)%buffer_src(:,:,1) ! initialize remap destination buffer + endif + deallocate(tmp_buffer) + else ! fid <= 0 + if (.not. ASSOCIATED(segment%field(m)%buffer_dst)) then + allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) + segment%field(m)%buffer_dst(:,:,:)=segment%field(m)%value + if (trim(segment%field(m)%name) == 'UO' .or. trim(segment%field(m)%name) == 'VO') then + allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc)) + segment%field(m)%bt_vel(:,:)=segment%field(m)%value + endif + endif + endif + + if (trim(segment%field(m)%name) == 'UO' .or. trim(segment%field(m)%name) == 'VO') then + if (segment%field(m)%fid>0) then ! calculate external BT velocity and transport if needed + if((trim(segment%field(m)%name) == 'UO' .and. (segment%direction == OBC_DIRECTION_W .or. & + segment%direction == OBC_DIRECTION_E)) .or. & + (trim(segment%field(m)%name) == 'VO' .and. (segment%direction == OBC_DIRECTION_N .or. & + segment%direction == OBC_DIRECTION_S))) then + do j=js_obc,je_obc + do i=is_obc,ie_obc + segment%normal_trans_bt(i,j) = 0.0 + do k=1,G%ke + segment%normal_vel(i,j,k) = segment%field(m)%buffer_dst(i,j,k) + segment%normal_trans(i,j,k) = segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k) + segment%normal_trans_bt(i,j)= segment%normal_trans_bt(i,j)+segment%normal_trans(i,j,k) + enddo + segment%normal_vel_bt(i,j) = segment%normal_trans_bt(i,j)/max(segment%Htot(i,j),1.e-12) + enddo + enddo + endif + endif + endif + + if (trim(segment%field(m)%name) == 'ZOS') then + do j=js_obc,je_obc + do i=is_obc,ie_obc + segment%eta(i,j) = segment%field(m)%buffer_dst(i,j,1) + enddo + enddo + endif + enddo + + + enddo ! end segment loop + +end subroutine update_OBC_segment_data + !> \namespace mom_open_boundary !! This module implements some aspects of internal open boundary !! conditions in MOM. diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 7d91be9f72..e713c5265b 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -33,7 +33,6 @@ module MOM_fixed_initialization use sloshing_initialization, only : sloshing_initialize_topography use seamount_initialization, only : seamount_initialize_topography use Phillips_initialization, only : Phillips_initialize_topography -use supercritical_initialization, only : supercritical_initialize_topography use netcdf @@ -91,7 +90,6 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) case ("none") case ("DOME") ! Avoid FATAL when using segments case ("tidal_bay") ; !Using segments now - case ("supercritical") ; !Using segments now case ("USER") ! Avoid FATAL when using segments case default ; call MOM_error(FATAL, "MOM_initialize_fixed: "// & "The open boundary positions specified by OBC_CONFIG="//& @@ -207,7 +205,6 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF) " \t DOME2D - use a shelf and slope configuration for the \n"//& " \t\t DOME2D gravity current/overflow test case. \n"//& " \t seamount - Gaussian bump for spontaneous motion test case.\n"//& - " \t supercritical - flat but with 8.95 degree land mask.\n"//& " \t Phillips - ACC-like idealized topography used in the Phillips config.\n"//& " \t USER - call a user modified routine.", & fail_if_missing=.true.) @@ -225,7 +222,6 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF) case ("sloshing"); call sloshing_initialize_topography(D, G, PF, max_depth) case ("seamount"); call seamount_initialize_topography(D, G, PF, max_depth) case ("Phillips"); call Phillips_initialize_topography(D, G, PF, max_depth) - case ("supercritical"); call supercritical_initialize_topography(D, G, PF, max_depth) case ("USER"); call user_initialize_topography(D, G, PF, max_depth) case default ; call MOM_error(FATAL,"MOM_initialize_topography: "// & "Unrecognized topography setup '"//trim(config)//"'") diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index c6a7728110..aa197618a7 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -449,6 +449,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & if (open_boundary_query(OBC, apply_specified_OBC=.true.)) then if (trim(config) == "DOME") then call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) + OBC%update_OBC = .false. elseif (trim(config) == "USER") then call user_set_OBC_data(OBC, tv, G, PF, tracer_Reg) elseif (.not. trim(config) == "none") then diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index d7b96bd20c..98cbd7669b 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -26,6 +26,7 @@ module DOME_initialization use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE +use MOM_open_boundary, only : OBC_segment_type use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -269,6 +270,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) character(len=40) :: mod = "DOME_set_OBC_data" ! This subroutine's name. integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB + type(OBC_segment_type), pointer :: segment is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -288,6 +290,13 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) Def_Rad = sqrt(D_edge*g_prime_tot) / (1.0e-4*1000.0) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5e3*Def_Rad) * GV%m_to_H + if (OBC%number_of_segments .ne. 1) then + print *, 'Error in DOME OBC segment setup' + return !!! Need a better error message here + endif + segment => OBC%OBC_segment_number(1) + if (.not. segment%on_pe) return + do k=1,nz rst = -1.0 if (k>1) rst = -1.0 + (real(k-1)-0.5)/real(nz-1) @@ -305,6 +314,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) (2.0 - Ri_trans)) if (k == nz) tr_k = tr_k + tr_0 * (2.0/(Ri_trans*(2.0+Ri_trans))) * & log((2.0+Ri_trans)/(2.0-Ri_trans)) + ! Old way do J=JsdB,JedB ; do i=isd,ied if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then ! This needs to be unneccesarily complicated without symmetric memory. @@ -317,6 +327,16 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) OBC%vh(i,J,k) = 0.0 ; OBC%v(i,J,k) = 0.0 endif enddo ; enddo + ! New way + isd = segment%HI%isd ; ied = segment%HI%ied + JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB + print *, 'DOME segment indices', isd, ied, JsdB, JedB +! do J=JsdB,JedB ; do i=isd,ied +! lon_im1 = 2.0*G%geoLonCv(i,J) - G%geoLonBu(I,J) +! segment%normal_trans(i,J,k) = tr_k * (exp(-2.0*(lon_im1 - 1000.0)/Def_Rad) -& +! exp(-2.0*(G%geoLonBu(I,J) - 1000.0)/Def_Rad)) +! segment%normal_vel(i,J,k) = v_k * exp(-2.0*(G%geoLonCv(i,J) - 1000.0)/Def_Rad) +! enddo ; enddo enddo ! The inflow values of temperature and salinity also need to be set here if diff --git a/src/user/supercritical_initialization.F90 b/src/user/supercritical_initialization.F90 index 99c272dce1..d4fcba29f9 100644 --- a/src/user/supercritical_initialization.F90 +++ b/src/user/supercritical_initialization.F90 @@ -110,6 +110,9 @@ subroutine supercritical_set_OBC_data(OBC, G, param_file) allocate(OBC%u(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) ; OBC%u(:,:,:) = 0.0 allocate(OBC%uh(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) ; OBC%uh(:,:,:) = 0.0 + if (.not.associated(OBC%ubt_outer)) then + allocate(OBC%ubt_outer(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%ubt_outer(:,:) = 0.0 + endif do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB if (OBC%OBC_segment_u(I,j)>0) then diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 index 95af6eeddd..c7f23b9c94 100644 --- a/src/user/tidal_bay_initialization.F90 +++ b/src/user/tidal_bay_initialization.F90 @@ -24,6 +24,7 @@ module tidal_bay_initialization use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE +use MOM_open_boundary, only : OBC_segment_type use MOM_verticalGrid, only : verticalGrid_type use MOM_time_manager, only : time_type, set_time, time_type_to_real @@ -49,8 +50,9 @@ subroutine tidal_bay_set_OBC_data(OBC, G, h, Time) real :: my_area, my_flux real :: PI character(len=40) :: mod = "tidal_bay_set_OBC_data" ! This subroutine's name. - integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n integer :: IsdB, IedB, JsdB, JedB + type(OBC_segment_type), pointer :: segment is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -75,6 +77,7 @@ subroutine tidal_bay_set_OBC_data(OBC, G, h, Time) enddo ; enddo my_flux = -tide_flow*SIN(2.0*PI*time_sec/(12.0*3600.0)) + ! Old way do j=jsd,jed ; do I=IsdB,IedB if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then OBC%eta_outer_u(I,j) = cff @@ -88,6 +91,17 @@ subroutine tidal_bay_set_OBC_data(OBC, G, h, Time) endif enddo ; enddo + ! New way +! do n = 1, OBC%number_of_segments +! segment => OBC%OBC_segment_number(n) + +! if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain + +! segment%normal_vel_bt(:,:) = my_flux/my_area +! segment%eta(:,:) = cff + +! enddo ! end segment loop + end subroutine tidal_bay_set_OBC_data !> \class tidal_bay_Initialization