From ec8ba6f13151b40656091a25d952db3f577e554e Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 16:20:21 -0400 Subject: [PATCH 01/13] Moved open_boundary_CS parameters into ocean_OBC_type - Although open_boundary_CS was a true "control structure" and ocean_OBC_type was a data container with public members I am moving the three parameters in open_boundary_CS into ocean_OBC_type and removed open_boundary_CS as part of a larger refactoring. - No answer changes. --- src/core/MOM_dynamics_legacy_split.F90 | 8 ++-- src/core/MOM_dynamics_split_RK2.F90 | 8 ++-- src/core/MOM_dynamics_unsplit.F90 | 4 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 4 +- src/core/MOM_open_boundary.F90 | 61 +++++++++++--------------- 5 files changed, 33 insertions(+), 52 deletions(-) diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 3ad5c3363c..944fe27e9a 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -111,7 +111,6 @@ module MOM_dynamics_legacy_split use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant @@ -241,7 +240,6 @@ module MOM_dynamics_legacy_split type(legacy_barotropic_CS), pointer :: barotropic_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to an open boundary ! condition type that specifies whether, where, and what open boundary ! conditions are used. If no open BCs are used, this pointer stays @@ -739,7 +737,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, CS%open_boundary_CSp) + v_old_rad_OBC, hp, h_old_rad_OBC, G) endif ! h_av = (h + hp)/2 @@ -998,7 +996,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, CS%open_boundary_CSp) + v_old_rad_OBC, h, h_old_rad_OBC, G) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. @@ -1382,7 +1380,7 @@ subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, GV, param_ if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp if (associated(OBC)) then CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) + call open_boundary_init(Time, G, param_file, diag, CS%OBC) endif if (.not. query_initialized(CS%eta,"sfc",restart_CS)) then diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index b65ab8281c..e107e98f2c 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -47,7 +47,6 @@ module MOM_dynamics_split_RK2 use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -161,7 +160,6 @@ module MOM_dynamics_split_RK2 type(barotropic_CS), pointer :: barotropic_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(tidal_forcing_CS), pointer :: tides_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() !< A pointer to an open boundary @@ -640,7 +638,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, CS%open_boundary_CSp) + v_old_rad_OBC, hp, h_old_rad_OBC, G) endif ! h_av = (h + hp)/2 @@ -852,7 +850,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, CS%open_boundary_CSp) + v_old_rad_OBC, h, h_old_rad_OBC, G) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. @@ -1143,7 +1141,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, param_fil if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp if (associated(OBC)) then CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) + call open_boundary_init(Time, G, param_file, diag, CS%OBC) endif if (.not. query_initialized(CS%eta,"sfc",restart_CS)) then diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 6c50d2b80a..d13d506072 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -104,7 +104,6 @@ module MOM_dynamics_unsplit use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -156,7 +155,6 @@ module MOM_dynamics_unsplit type(PressureForce_CS), pointer :: PressureForce_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to an open boundary ! condition type that specifies whether, where, and what open boundary ! conditions are used. If no open BCs are used, this pointer stays @@ -682,7 +680,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, param_file, diag, CS, & if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp if (associated(OBC)) then CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) + call open_boundary_init(Time, G, param_file, diag, CS%OBC) endif flux_units = get_flux_units(GV) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 7055547e0a..f6582787c3 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -102,7 +102,6 @@ module MOM_dynamics_unsplit_RK2 use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -162,7 +161,6 @@ module MOM_dynamics_unsplit_RK2 type(PressureForce_CS), pointer :: PressureForce_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to an open boundary ! condition type that specifies whether, where, and what open boundary ! conditions are used. If no open BCs are used, this pointer stays @@ -646,7 +644,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, param_file, diag, CS if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp if (associated(OBC)) then CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) + call open_boundary_init(Time, G, param_file, diag, CS%OBC) endif flux_units = get_flux_units(GV) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index ea419ae4c6..f0bffa6d49 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -14,21 +14,6 @@ module MOM_open_boundary public Radiation_Open_Bdry_Conds, open_boundary_init, open_boundary_end -!> The control structure for open-boundaries -type, public :: open_boundary_CS ; private - real :: gamma_uv !< The relative weighting for the baroclinic radiation - !! velocities (or speed of characteristics) at the - !! new time level (1) or the running mean (0) for velocities. - !! Valid values range from 0 to 1, with a default of 0.3. - real :: gamma_h !< The relative weighting for the baroclinic radiation - !! velocities (or speed of characteristics) at the - !! new time level (1) or the running mean (0) for thicknesses. - !! Valid values range from 0 to 1, with a default of 0.2. - real :: rx_max !< The maximum magnitude of the baroclinic radiation - !! velocity (or speed of characteristics), in m s-1. The - !! default value is 10 m s-1. -end type open_boundary_CS - integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 integer, parameter, public :: OBC_FLATHER_E = 4, OBC_FLATHER_W = 5 integer, parameter, public :: OBC_FLATHER_N = 6, OBC_FLATHER_S = 7 @@ -78,6 +63,19 @@ module MOM_open_boundary v => NULL(), & !< The prescribed values of the meridional velocity (v) at OBC points. uh => NULL(), & !< The prescribed values of the zonal volume transport (uh) at OBC points. vh => NULL() !< The prescribed values of the meridional volume transport (vh) at OBC points. + + ! The following parameters are used in the baroclinic radiation code: + real :: gamma_uv !< The relative weighting for the baroclinic radiation + !! velocities (or speed of characteristics) at the + !! new time level (1) or the running mean (0) for velocities. + !! Valid values range from 0 to 1, with a default of 0.3. + real :: gamma_h !< The relative weighting for the baroclinic radiation + !! velocities (or speed of characteristics) at the + !! new time level (1) or the running mean (0) for thicknesses. + !! Valid values range from 0 to 1, with a default of 0.2. + real :: rx_max !< The maximum magnitude of the baroclinic radiation + !! velocity (or speed of characteristics), in m s-1. The + !! default value is 10 m s-1. end type ocean_OBC_type integer :: id_clock_pass @@ -90,16 +88,15 @@ module MOM_open_boundary !> Diagnose radiation conditions at open boundaries subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & - h_new, h_old, G, CS) + h_new, h_old, G) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure - type(ocean_OBC_type), pointer :: OBC !< Open boundary data + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v_old real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_new real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old - type(open_boundary_CS), pointer :: CS !< Open boundary control structure ! Local variables real :: dhdt, dhdx, gamma_u, gamma_h, gamma_v real :: rx_max, ry_max ! coefficients for radiation @@ -112,11 +109,9 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & if (.not.(OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west .or. & OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south)) & return - if (.not.associated(CS)) call MOM_error(FATAL, & - "MOM_open_boundary: Module must be initialized before it is used.") - gamma_u = CS%gamma_uv ; gamma_v = CS%gamma_uv ; gamma_h = CS%gamma_h - rx_max = CS%rx_max ; ry_max = CS%rx_max + gamma_u = OBC%gamma_uv ; gamma_v = OBC%gamma_uv ; gamma_h = OBC%gamma_h + rx_max = OBC%rx_max ; ry_max = OBC%rx_max if (OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west) then do k=1,nz ; do j=js,je ; do I=is-1,ie ; if (OBC%OBC_mask_u(I,j)) then @@ -206,20 +201,15 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & end subroutine Radiation_Open_Bdry_Conds !> Initialize open boundary control structure -subroutine open_boundary_init(Time, G, param_file, diag, CS) +subroutine open_boundary_init(Time, G, param_file, diag, OBC) type(time_type), target, intent(in) :: Time !< Current model time type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure - type(open_boundary_CS), pointer :: CS !< Open boundary control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables logical :: flather_east, flather_west, flather_north, flather_south - if (associated(CS)) then - call MOM_error(WARNING, "MOM_open_boundary: open_boundary_init called with associated control structure.") - return - endif - call log_version(param_file, mod, version) call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", flather_east, & "If true, some zonal velocity points use Flather open \n"//& @@ -240,20 +230,19 @@ subroutine open_boundary_init(Time, G, param_file, diag, CS) if (.not.(flather_east .or. flather_west .or. flather_north .or. & flather_south)) return - allocate(CS) - call get_param(param_file, mod, "OBC_RADIATION_MAX", CS%rx_max, & + call get_param(param_file, mod, "OBC_RADIATION_MAX", OBC%rx_max, & "The maximum magnitude of the baroclinic radiation \n"//& "velocity (or speed of characteristics). This is only \n"//& "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & units="m s-1", default=10.0) - call get_param(param_file, mod, "OBC_RAD_VEL_WT", CS%gamma_uv, & + call get_param(param_file, mod, "OBC_RAD_VEL_WT", OBC%gamma_uv, & "The relative weighting for the baroclinic radiation \n"//& "velocities (or speed of characteristics) at the new \n"//& "time level (1) or the running mean (0) for velocities. \n"//& "Valid values range from 0 to 1. This is only used if \n"//& "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & units="nondim", default=0.3) - call get_param(param_file, mod, "OBC_RAD_THICK_WT", CS%gamma_h, & + call get_param(param_file, mod, "OBC_RAD_THICK_WT", OBC%gamma_h, & "The relative weighting for the baroclinic radiation \n"//& "velocities (or speed of characteristics) at the new \n"//& "time level (1) or the running mean (0) for thicknesses. \n"//& @@ -266,9 +255,9 @@ subroutine open_boundary_init(Time, G, param_file, diag, CS) end subroutine open_boundary_init !> Deallocate open boundary data -subroutine open_boundary_end(CS) - type(open_boundary_CS), pointer :: CS !< Open boundary control structure - deallocate(CS) +subroutine open_boundary_end(OBC) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + deallocate(OBC) end subroutine open_boundary_end !> \namespace mom_open_boundary From 98c4918e25ca33fe0c48b0e74f7f78a362eb5f24 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 18:07:48 -0400 Subject: [PATCH 02/13] Moved calls to open_boundary_init() from MOM_dynamics_* to MOM.F90 - Each algorithm had it's own call to open_boundary_init() which simply set three parameters. I've moved this up to early in the main initialization sequence ready to re-purpose the s/r. - No answer changes. --- src/core/MOM.F90 | 2 +- src/core/MOM_dynamics_legacy_split.F90 | 7 ++----- src/core/MOM_dynamics_split_RK2.F90 | 7 ++----- src/core/MOM_dynamics_unsplit.F90 | 7 ++----- src/core/MOM_dynamics_unsplit_RK2.F90 | 7 ++----- 5 files changed, 9 insertions(+), 21 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index cc2f199989..3221d0371d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1799,6 +1799,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, param_file, & dirs, CS%restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) + call open_boundary_init(Time, G, param_file, diag, CS%OBC) call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") @@ -1916,7 +1917,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call wave_speed_init(Time, G, param_file, diag, CS%wave_speed_CSp) call VarMix_init(Time, G, param_file, diag, CS%VarMix, CS%wave_speed_CSp) call set_visc_init(Time, G, GV, param_file, diag, CS%visc, CS%set_visc_CSp) - if (CS%split) then allocate(eta(SZI_(G),SZJ_(G))) ; eta(:,:) = 0.0 if (CS%legacy_split) then diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 944fe27e9a..77b9f3c6f5 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -110,7 +110,7 @@ module MOM_dynamics_legacy_split use MOM_interface_heights, only : find_eta use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant @@ -1378,10 +1378,7 @@ subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, GV, param_ CS%set_visc_CSp => setVisc_CSp if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%OBC) - endif + if (associated(OBC)) CS%OBC => OBC if (.not. query_initialized(CS%eta,"sfc",restart_CS)) then ! Estimate eta based on the layer thicknesses - h. With the Boussinesq diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index e107e98f2c..2e2173a502 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -46,7 +46,7 @@ module MOM_dynamics_split_RK2 use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -1139,10 +1139,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, param_fil (LEN_TRIM(dirs%input_filename) == 1)) ) if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%OBC) - endif + if (associated(OBC)) CS%OBC => OBC if (.not. query_initialized(CS%eta,"sfc",restart_CS)) then ! Estimate eta based on the layer thicknesses - h. With the Boussinesq diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index d13d506072..fc392cbb88 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -103,7 +103,7 @@ module MOM_dynamics_unsplit use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -678,10 +678,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, param_file, diag, CS, & CS%set_visc_CSp => setVisc_CSp if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%OBC) - endif + if (associated(OBC)) CS%OBC => OBC flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index f6582787c3..8e228a9189 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -101,7 +101,7 @@ module MOM_dynamics_unsplit_RK2 use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -642,10 +642,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, param_file, diag, CS CS%set_visc_CSp => setVisc_CSp if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%OBC) - endif + if (associated(OBC)) CS%OBC => OBC flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & From 89b858f1d8d4b50927008c9ef85949c9be5d535a Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 18:14:31 -0400 Subject: [PATCH 03/13] Re-ordered s/r in MOM_open_boundary.F90 - Trivial re-factor: changed order of subroutines to the same order as listed in the public list, ready for adding new subroutines in the logical place in the file. - No answer changes. --- src/core/MOM_open_boundary.F90 | 124 +++++++++++++++++---------------- 1 file changed, 63 insertions(+), 61 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f0bffa6d49..1e769046e5 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -12,7 +12,9 @@ module MOM_open_boundary #include -public Radiation_Open_Bdry_Conds, open_boundary_init, open_boundary_end +public open_boundary_init +public open_boundary_end +public Radiation_Open_Bdry_Conds integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 integer, parameter, public :: OBC_FLATHER_E = 4, OBC_FLATHER_W = 5 @@ -86,6 +88,66 @@ module MOM_open_boundary contains +!> Initialize open boundary control structure +subroutine open_boundary_init(Time, G, param_file, diag, OBC) + type(time_type), target, intent(in) :: Time !< Current model time + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + ! Local variables + logical :: flather_east, flather_west, flather_north, flather_south + + call log_version(param_file, mod, version) + call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", flather_east, & + "If true, some zonal velocity points use Flather open \n"//& + "boundary conditions on the east side of the ocean.", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", flather_west, & + "If true, some zonal velocity points use Flather open \n"//& + "boundary conditions on the west side of the ocean.", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", flather_north, & + "If true, some meridional velocity points use Flather \n"//& + "open boundary conditions on the north side of the ocean.", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", flather_south, & + "If true, some meridional velocity points use Flather \n"//& + "open boundary conditions on the north side of the ocean.", & + default=.false.) + if (.not.(flather_east .or. flather_west .or. flather_north .or. & + flather_south)) return + + call get_param(param_file, mod, "OBC_RADIATION_MAX", OBC%rx_max, & + "The maximum magnitude of the baroclinic radiation \n"//& + "velocity (or speed of characteristics). This is only \n"//& + "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="m s-1", default=10.0) + call get_param(param_file, mod, "OBC_RAD_VEL_WT", OBC%gamma_uv, & + "The relative weighting for the baroclinic radiation \n"//& + "velocities (or speed of characteristics) at the new \n"//& + "time level (1) or the running mean (0) for velocities. \n"//& + "Valid values range from 0 to 1. This is only used if \n"//& + "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="nondim", default=0.3) + call get_param(param_file, mod, "OBC_RAD_THICK_WT", OBC%gamma_h, & + "The relative weighting for the baroclinic radiation \n"//& + "velocities (or speed of characteristics) at the new \n"//& + "time level (1) or the running mean (0) for thicknesses. \n"//& + "Valid values range from 0 to 1. This is only used if \n"//& + "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="nondim", default=0.2) + + id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) + +end subroutine open_boundary_init + +!> Deallocate open boundary data +subroutine open_boundary_end(OBC) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + deallocate(OBC) +end subroutine open_boundary_end + !> Diagnose radiation conditions at open boundaries subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & h_new, h_old, G) @@ -200,66 +262,6 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & end subroutine Radiation_Open_Bdry_Conds -!> Initialize open boundary control structure -subroutine open_boundary_init(Time, G, param_file, diag, OBC) - type(time_type), target, intent(in) :: Time !< Current model time - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(param_file_type), intent(in) :: param_file !< Parameter file handle - type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure - type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - ! Local variables - logical :: flather_east, flather_west, flather_north, flather_south - - call log_version(param_file, mod, version) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", flather_east, & - "If true, some zonal velocity points use Flather open \n"//& - "boundary conditions on the east side of the ocean.", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", flather_west, & - "If true, some zonal velocity points use Flather open \n"//& - "boundary conditions on the west side of the ocean.", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", flather_north, & - "If true, some meridional velocity points use Flather \n"//& - "open boundary conditions on the north side of the ocean.", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", flather_south, & - "If true, some meridional velocity points use Flather \n"//& - "open boundary conditions on the north side of the ocean.", & - default=.false.) - if (.not.(flather_east .or. flather_west .or. flather_north .or. & - flather_south)) return - - call get_param(param_file, mod, "OBC_RADIATION_MAX", OBC%rx_max, & - "The maximum magnitude of the baroclinic radiation \n"//& - "velocity (or speed of characteristics). This is only \n"//& - "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="m s-1", default=10.0) - call get_param(param_file, mod, "OBC_RAD_VEL_WT", OBC%gamma_uv, & - "The relative weighting for the baroclinic radiation \n"//& - "velocities (or speed of characteristics) at the new \n"//& - "time level (1) or the running mean (0) for velocities. \n"//& - "Valid values range from 0 to 1. This is only used if \n"//& - "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="nondim", default=0.3) - call get_param(param_file, mod, "OBC_RAD_THICK_WT", OBC%gamma_h, & - "The relative weighting for the baroclinic radiation \n"//& - "velocities (or speed of characteristics) at the new \n"//& - "time level (1) or the running mean (0) for thicknesses. \n"//& - "Valid values range from 0 to 1. This is only used if \n"//& - "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="nondim", default=0.2) - - id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) - -end subroutine open_boundary_init - -!> Deallocate open boundary data -subroutine open_boundary_end(OBC) - type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - deallocate(OBC) -end subroutine open_boundary_end - !> \namespace mom_open_boundary !! This module implements some aspects of internal open boundary !! conditions in MOM. From 83e50eb0dfe8454a0431a079bf716c1bac9a9b37 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 20:26:35 -0400 Subject: [PATCH 04/13] Split out DOME OBC positions - Previously, s/r DOME_set_Open_Bdry_Conds() set both the OBC masks and the OBC data (inflow conditions). - DOME_set_OBC_positions() now sets the OBC masks. - DOME_set_OBC_data() replaces DOME_set_Open_Bdry_Conds() and only sets inflow data. - This is one step in a larger re-factor of OBC code. --- .../MOM_state_initialization.F90 | 5 +- src/user/DOME_initialization.F90 | 138 ++++++++---------- 2 files changed, 67 insertions(+), 76 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0f9568cd2c..fe310f51a3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -42,7 +42,7 @@ module MOM_state_initialization use user_initialization, only : user_init_temperature_salinity use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges use DOME_initialization, only : DOME_initialize_thickness -use DOME_initialization, only : DOME_set_Open_Bdry_Conds +use DOME_initialization, only : DOME_set_OBC_positions, DOME_set_OBC_data use DOME_initialization, only : DOME_initialize_sponges use ISOMIP_initialization, only : ISOMIP_initialize_thickness use ISOMIP_initialization, only : ISOMIP_initialize_sponges @@ -440,7 +440,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & " \t USER - call a user modified routine.", default="file", & fail_if_missing=.true.) if (trim(config) == "DOME") then - call DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg) + call DOME_set_OBC_positions(G, PF, OBC) + call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) elseif (trim(config) == "USER") then call user_set_Open_Bdry_Conds(OBC, tv, G, PF, tracer_Reg) else diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index b94fb69a18..455cb559c5 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -37,7 +37,8 @@ module DOME_initialization public DOME_initialize_topography public DOME_initialize_thickness public DOME_initialize_sponges -public DOME_set_Open_Bdry_Conds +public DOME_set_OBC_positions +public DOME_set_OBC_data contains @@ -230,31 +231,69 @@ subroutine DOME_initialize_sponges(G, GV, tv, PF, CSp) endif end subroutine DOME_initialize_sponges -! ----------------------------------------------------------------------------- -! ----------------------------------------------------------------------------- +!> Set the positions of the open boundary needed for the DOME experiment. +subroutine DOME_set_OBC_positions(G, param_file, OBC) + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(param_file_type), intent(in) :: param_file !< Parameter file handle. + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure. + ! Local variables + character(len=40) :: mod = "DOME_set_OBC_positions" ! This subroutine's name. + integer :: i, j + logical :: any_OBC ! Set to true if any points in this subdomain use OBCs + + if (.not.associated(OBC)) allocate(OBC) + + call get_param(param_file, mod, "APPLY_OBC_U", OBC%apply_OBC_u, & + "If true, open boundary conditions may be set at some \n"//& + "u-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V", OBC%apply_OBC_v, & + "If true, open boundary conditions may be set at some \n"//& + "v-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + + any_OBC = .false. + if (OBC%apply_OBC_u) then + ! Set where u points are determined by OBCs. + !allocate(OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC_mask_u(:,:) = .false. + call MOM_error(FATAL,"DOME_initialization, DOME_set_OBC_positions: "//& + "APPLY_OBC_U=True is not coded for the DOME experiment") + endif + if (OBC%apply_OBC_v) then + ! Set where v points are determined by OBCs. + allocate(OBC%OBC_mask_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_mask_v(:,:) = .false. + do J=G%JsdB,G%JedB ; do i=G%isd,G%ied + if ((G%geoLonCv(i,J) > 1000.0) .and. (G%geoLonCv(i,J) < 1100.0) .and. & + (abs(G%geoLatCv(i,J) - G%gridLatB(G%JegB)) < 0.1)) then + OBC%OBC_mask_v(i,J) = .true. + any_OBC = .true. + endif + enddo ; enddo + endif + if (.not.any_OBC) then + ! If this PE does not have any OBC points then we do not need the mask + OBC%apply_OBC_v = .false. + deallocate(OBC%OBC_mask_v) + endif +end subroutine DOME_set_OBC_positions + !> This subroutine sets the properties of flow at open boundary conditions. !! This particular example is for the DOME inflow describe in Legg et al. 2006. -subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) +subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies !! whether, where, and what open boundary !! conditions are used. type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any !! available thermodynamic fields, including potential - !! temperature and salinity or mixed layer density. Absent - !! fields have NULL ptrs. + !! temperature and salinity or mixed layer density. Absent + !! fields have NULL ptrs. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry. - logical :: any_OBC ! Set to true if any points in this subdomain use - ! open boundary conditions. - logical, pointer, dimension(:,:) :: & - OBC_mask_u => NULL(), & ! These arrays are true at zonal or meridional - OBC_mask_v => NULL() ! velocity points that have prescribed open boundary - ! conditions. 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 @@ -276,7 +315,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) ! thickness D_edge, in the same units as lat. real :: Ri_trans ! The shear Richardson number in the transition ! region of the specified shear profile. - character(len=40) :: mod = "DOME_set_Open_Bdry_Conds" ! This subroutine's name. + 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 @@ -289,56 +328,10 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) Ri_trans = 1.0/3.0 ! The shear Richardson number in the transition region ! region of the specified shear profile. - call get_param(param_file, mod, "APPLY_OBC_U", apply_OBC_u, & - "If true, open boundary conditions may be set at some \n"//& - "u-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V", apply_OBC_v, & - "If true, open boundary conditions may be set at some \n"//& - "v-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - - if (apply_OBC_u) then - ! Determine where u points are applied. - allocate(OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC_mask_u(:,:) = .false. - any_OBC = .false. - do j=jsd,jed ; do I=IsdB,IedB - ! if (SOME_TEST_FOR_U_OPEN_BCS) then - ! OBC_mask_u(I,j) = .true. ; any_OBC = .true. - ! endif - enddo ; enddo - if (.not.any_OBC) then - ! This processor has no u points at which open boundary conditions are - ! to be applied. - apply_OBC_u = .false. - deallocate(OBC_mask_u) - endif - endif - if (apply_OBC_v) then - ! Determine where v points are applied. - allocate(OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC_mask_v(:,:) = .false. - any_OBC = .false. - do J=JsdB,JedB ; do i=isd,ied - if ((G%geoLonCv(i,J) > 1000.0) .and. (G%geoLonCv(i,J) < 1100.0) .and. & - (abs(G%geoLatCv(i,J) - G%gridLatB(G%JegB)) < 0.1)) then - OBC_mask_v(i,J) = .true. ; any_OBC = .true. - endif - enddo ; enddo - if (.not.any_OBC) then - ! This processor has no v points at which open boundary conditions are - ! to be applied. - apply_OBC_v = .false. - deallocate(OBC_mask_v) - endif - endif - - if (.not.(apply_OBC_u .or. apply_OBC_v)) return + if (.not.associated(OBC)) return + if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v)) return - if (.not.associated(OBC)) allocate(OBC) - - if (apply_OBC_u) then - OBC%apply_OBC_u = .true. - OBC%OBC_mask_u => OBC_mask_u + if (OBC%apply_OBC_u) then allocate(OBC%u(IsdB:IedB,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0 allocate(OBC%uh(IsdB:IedB,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0 allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE @@ -346,9 +339,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) if (OBC%OBC_mask_u(I,j)) OBC%OBC_kind_u(I,j) = OBC_SIMPLE enddo ; enddo endif - if (apply_OBC_v) then - OBC%apply_OBC_v = .true. - OBC%OBC_mask_v => OBC_mask_v + if (OBC%apply_OBC_v) then allocate(OBC%v(isd:ied,JsdB:JedB,nz)) ; OBC%v(:,:,:) = 0.0 allocate(OBC%vh(isd:ied,JsdB:JedB,nz)) ; OBC%vh(:,:,:) = 0.0 allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE @@ -357,7 +348,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) enddo ; enddo endif - if (apply_OBC_v) then + if (OBC%apply_OBC_v) then g_prime_tot = (G%g_Earth/GV%Rho0)*2.0 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 @@ -380,7 +371,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) 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)) do J=JsdB,JedB ; do i=isd,ied - if (OBC_mask_v(i,J)) then + if (OBC%OBC_mask_v(i,J)) then ! This needs to be unneccesarily complicated without symmetric memory. lon_im1 = 2.0*G%geoLonCv(i,J) - G%geoLonBu(I,J) ! if (isd > IsdB) lon_im1 = G%geoLonBu(I-1,J) @@ -394,9 +385,9 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) enddo endif - if (apply_OBC_u) then + if (OBC%apply_OBC_u) then do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB - if (OBC_mask_u(I,j)) then + if (OBC%OBC_mask_u(I,j)) then ! An appropriate expression for the zonal inflow velocities and ! transports should go here. OBC%uh(I,j,k) = 0.0 * GV%m_to_H ; OBC%u(I,j,k) = 0.0 @@ -408,7 +399,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) ! The inflow values of temperature and salinity also need to be set here if ! these variables are used. The following code is just a naive example. - if (apply_OBC_u .or. apply_OBC_v) then + if (OBC%apply_OBC_u .or. OBC%apply_OBC_v) then if (associated(tv%S)) then ! In this example, all S inflows have values of 35 psu. call add_tracer_OBC_values("S", tr_Reg, OBC_inflow=35.0) @@ -428,13 +419,13 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) do k=1,nz ; T0(k) = T0(k) + (GV%Rlay(k)-rho_guess(k)) / drho_dT(k) ; enddo enddo - if (apply_OBC_u) then + if (OBC%apply_OBC_u) then allocate(OBC_T_u(IsdB:IedB,jsd:jed,nz)) do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB OBC_T_u(I,j,k) = T0(k) enddo ; enddo ; enddo endif - if (apply_OBC_v) then + if (OBC%apply_OBC_v) then allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied OBC_T_v(i,J,k) = T0(k) @@ -445,8 +436,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) endif endif -end subroutine DOME_set_Open_Bdry_Conds -! ----------------------------------------------------------------------------- +end subroutine DOME_set_OBC_data !> \class DOME_initialization !! From 2d3856f811041ab88b5368872c3db9630064dead Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 20:42:18 -0400 Subject: [PATCH 05/13] Split out OBC positions in USER_initialization - As for commit 83e50eb0dfe8454a0431a079bf716c1bac9a9b37 the setting OBC masks and inflow data are now separated. - This commit is mostly a no-op since the USER_initialization module is blank and only provides example APIs. - This is one step in a larger re-factor of OBC code. --- .../MOM_state_initialization.F90 | 6 +++-- src/user/user_initialization.F90 | 25 ++++++++++++++++--- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index fe310f51a3..7bda870db3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -40,7 +40,8 @@ module MOM_state_initialization use MOM_EOS, only : int_specific_vol_dp use user_initialization, only : user_initialize_thickness, user_initialize_velocity use user_initialization, only : user_init_temperature_salinity -use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges +use user_initialization, only : user_set_OBC_positions, user_set_OBC_data +use user_initialization, only : user_initialize_sponges use DOME_initialization, only : DOME_initialize_thickness use DOME_initialization, only : DOME_set_OBC_positions, DOME_set_OBC_data use DOME_initialization, only : DOME_initialize_sponges @@ -443,7 +444,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & call DOME_set_OBC_positions(G, PF, OBC) call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) elseif (trim(config) == "USER") then - call user_set_Open_Bdry_Conds(OBC, tv, G, PF, tracer_Reg) + call user_set_OBC_positions(G, PF, OBC) + call user_set_OBC_data(OBC, tv, G, PF, tracer_Reg) else call MOM_error(FATAL, "The open boundary conditions specified by "//& "OBC_CONFIG = "//trim(config)//" have not been fully implemented.") diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 0364ef884c..33e1090617 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -40,7 +40,7 @@ module user_initialization public USER_set_coord, USER_initialize_topography, USER_initialize_thickness public USER_initialize_velocity, USER_init_temperature_salinity public USER_init_mixed_layer_density, USER_initialize_sponges -public USER_set_Open_Bdry_Conds, USER_set_rotation +public USER_set_OBC_positions, USER_set_OBC_data, USER_set_rotation logical :: first_call = .true. @@ -197,8 +197,25 @@ subroutine USER_initialize_sponges(G, use_temperature, tv, param_file, CSp, h) end subroutine USER_initialize_sponges +!> This subroutine sets the location of open boundaries. +subroutine USER_set_OBC_positions(G, param_file, OBC) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(param_file_type), intent(in) :: param_file !< A structure indicating the + !! open file to parse for model + !! parameter values. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies + !! whether, where, and what open boundary + !! conditions are used. +! call MOM_error(FATAL, & +! "USER_initialization.F90, USER_set_OBC_positions: " // & +! "Unmodified user routine called - you must edit the routine to use it") + + if (first_call) call write_user_log(param_file) + +end subroutine USER_set_OBC_positions + !> This subroutine sets the properties of flow at open boundary conditions. -subroutine USER_set_Open_Bdry_Conds(OBC, tv, G, param_file, tr_Reg) +subroutine USER_set_OBC_data(OBC, tv, G, param_file, tr_Reg) type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies !! whether, where, and what open boundary !! conditions are used. @@ -212,12 +229,12 @@ subroutine USER_set_Open_Bdry_Conds(OBC, tv, G, param_file, tr_Reg) !! parameter values. type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry. ! call MOM_error(FATAL, & -! "USER_initialization.F90, USER_set_Open_Bdry_Conds: " // & +! "USER_initialization.F90, USER_set_OBC_data: " // & ! "Unmodified user routine called - you must edit the routine to use it") if (first_call) call write_user_log(param_file) -end subroutine USER_set_Open_Bdry_Conds +end subroutine USER_set_OBC_data subroutine USER_set_rotation(G, param_file) type(ocean_grid_type), intent(inout) :: G From d5efae284df5f1c7d92bde2e1736b1d49f0a5b67 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 22:11:09 -0400 Subject: [PATCH 06/13] Moved calls that set OBC masks into MOM_fixed_initialization - Prior to this commit, OBC masks were set at the same time as OBC inflow data was set, and all called from MOM_state_initialization, which is too late to correct bathymetry etc. - This is part of a larger OBC re-factor. - No answer changes. --- src/core/MOM.F90 | 2 +- .../MOM_fixed_initialization.F90 | 37 +++++++++++++++++-- .../MOM_state_initialization.F90 | 10 +---- 3 files changed, 35 insertions(+), 14 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3221d0371d..c2156d9cbf 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1785,7 +1785,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call callTree_waypoint("restart registration complete (initialize_MOM)") call cpu_clock_begin(id_clock_MOM_init) - call MOM_initialize_fixed(G, param_file, write_geom_files, dirs%output_directory) + call MOM_initialize_fixed(G, CS%OBC, param_file, write_geom_files, dirs%output_directory) call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)") call MOM_initialize_coord(GV, param_file, write_geom_files, & dirs%output_directory, CS%tv, G%max_depth) diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 5ebe25a844..1d762e7376 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -18,9 +18,10 @@ module MOM_fixed_initialization use MOM_io, only : slasher, vardesc, write_field, var_desc use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_grid_initialize, only : initialize_masks, set_grid_metrics +use MOM_open_boundary, only : ocean_OBC_type use MOM_string_functions, only : uppercase -use user_initialization, only : user_initialize_topography -use DOME_initialization, only : DOME_initialize_topography +use user_initialization, only : user_initialize_topography, USER_set_OBC_positions +use DOME_initialization, only : DOME_initialize_topography, DOME_set_OBC_positions use ISOMIP_initialization, only : ISOMIP_initialize_topography use benchmark_initialization, only : benchmark_initialize_topography use DOME2d_initialization, only : DOME2d_initialize_topography @@ -43,8 +44,9 @@ module MOM_fixed_initialization ! ----------------------------------------------------------------------------- !> MOM_initialize_fixed sets up time-invariant quantities related to MOM6's !! horizontal grid, bathymetry, and the Coriolis parameter. -subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir) +subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure. type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: write_geom !< If true, write grid geometry files. @@ -53,7 +55,7 @@ subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir) ! Local character(len=200) :: inputdir ! The directory where NetCDF input files are. character(len=200) :: config - logical :: debug + logical :: debug, apply_OBC_u, apply_OBC_v ! This include declares and sets the variable "version". #include "version_variable.h" @@ -78,6 +80,33 @@ subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir) ! masks, and Coriolis parameter. ! ==================================================================== +! Determine the position of any open boundaries + call get_param(PF, mod, "APPLY_OBC_U", apply_OBC_u, & + "If true, open boundary conditions may be set at some \n"//& + "u-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + call get_param(PF, mod, "APPLY_OBC_V", apply_OBC_v, & + "If true, open boundary conditions may be set at some \n"//& + "v-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + if (apply_OBC_u .or. apply_OBC_v) then + call get_param(PF, mod, "OBC_CONFIG", config, & + "A string that sets how the open boundary conditions are \n"//& + " configured: \n"//& + " \t DOME - use a slope and channel configuration for the \n"//& + " \t\t DOME sill-overflow test case. \n"//& + " \t USER - call a user modified routine.", default="file", & + fail_if_missing=.true.) + select case ( trim(config) ) + case ("none") + case ("DOME") ; call DOME_set_OBC_positions(G, PF, OBC) + case ("USER") ; call user_set_OBC_positions(G, PF, OBC) + case default ; call MOM_error(FATAL, "MOM_initialize_fixed: "// & + "The open boundary positions specified by OBC_CONFIG="//& + trim(config)//" have not been fully implemented.") + end select + endif + ! This call sets seamasks that prohibit flow over any point with ! ! a bottom that is shallower than min_depth from PF. ! call initialize_masks(G, PF) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7bda870db3..a846988fde 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -433,18 +433,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & "v-points, with the configuration controlled by OBC_CONFIG", & default=.false.) if (apply_OBC_u .or. apply_OBC_v) then - call get_param(PF, mod, "OBC_CONFIG", config, & - "A string that sets how the open boundary conditions are \n"//& - " configured: \n"//& - " \t DOME - use a slope and channel configuration for the \n"//& - " \t\t DOME sill-overflow test case. \n"//& - " \t USER - call a user modified routine.", default="file", & - fail_if_missing=.true.) + call get_param(PF, mod, "OBC_CONFIG", config, fail_if_missing=.true., do_not_log=.true.) if (trim(config) == "DOME") then - call DOME_set_OBC_positions(G, PF, OBC) call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) elseif (trim(config) == "USER") then - call user_set_OBC_positions(G, PF, OBC) call user_set_OBC_data(OBC, tv, G, PF, tracer_Reg) else call MOM_error(FATAL, "The open boundary conditions specified by "//& From c90bec4c9438237e05fc2e80e5a49ebe9a75322d Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 23:35:02 -0400 Subject: [PATCH 07/13] Moved allocation of OBC type to MOM_open_boundary - Previously, allocation of the OBC type was handled in user code rather than centrally. This meant duplication and inconsistent get_param for the controlling run-time parameters. - I've created an open_boundary_config() that allocates and reads the main run-time parameters that can be called from MOM_fixed_init. - To keep things succinct there is also now an open_boundary_query() that return true if collections of flags were set. - No answer changes. --- src/core/MOM_open_boundary.F90 | 120 ++++++++++++------ .../MOM_fixed_initialization.F90 | 14 +- src/user/DOME_initialization.F90 | 12 +- 3 files changed, 88 insertions(+), 58 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 1e769046e5..c7196fe710 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -12,7 +12,9 @@ module MOM_open_boundary #include +public open_boundary_config public open_boundary_init +public open_boundary_query public open_boundary_end public Radiation_Open_Bdry_Conds @@ -88,60 +90,102 @@ module MOM_open_boundary contains -!> Initialize open boundary control structure -subroutine open_boundary_init(Time, G, param_file, diag, OBC) - type(time_type), target, intent(in) :: Time !< Current model time +!> Enables OBC module and reads configuration parameters +subroutine open_boundary_config(G, param_file, OBC) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle - type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables logical :: flather_east, flather_west, flather_north, flather_south + allocate(OBC) + call log_version(param_file, mod, version) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", flather_east, & - "If true, some zonal velocity points use Flather open \n"//& - "boundary conditions on the east side of the ocean.", & + call get_param(param_file, mod, "APPLY_OBC_U", OBC%apply_OBC_u, & + "If true, open boundary conditions may be set at some \n"//& + "u-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V", OBC%apply_OBC_v, & + "If true, open boundary conditions may be set at some \n"//& + "v-points, with the configuration controlled by OBC_CONFIG", & default=.false.) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", flather_west, & - "If true, some zonal velocity points use Flather open \n"//& - "boundary conditions on the west side of the ocean.", & + call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", OBC%apply_OBC_u_flather_east, & + "Apply a Flather open boundary condition on the eastern\n"//& + "side of the global domain", & default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", flather_north, & - "If true, some meridional velocity points use Flather \n"//& - "open boundary conditions on the north side of the ocean.", & + call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", OBC%apply_OBC_u_flather_west, & + "Apply a Flather open boundary condition on the western\n"//& + "side of the global domain", & default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", flather_south, & - "If true, some meridional velocity points use Flather \n"//& - "open boundary conditions on the north side of the ocean.", & + call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", OBC%apply_OBC_v_flather_north, & + "Apply a Flather open boundary condition on the northern\n"//& + "side of the global domain", & default=.false.) - if (.not.(flather_east .or. flather_west .or. flather_north .or. & - flather_south)) return - - call get_param(param_file, mod, "OBC_RADIATION_MAX", OBC%rx_max, & - "The maximum magnitude of the baroclinic radiation \n"//& - "velocity (or speed of characteristics). This is only \n"//& - "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="m s-1", default=10.0) - call get_param(param_file, mod, "OBC_RAD_VEL_WT", OBC%gamma_uv, & - "The relative weighting for the baroclinic radiation \n"//& - "velocities (or speed of characteristics) at the new \n"//& - "time level (1) or the running mean (0) for velocities. \n"//& - "Valid values range from 0 to 1. This is only used if \n"//& - "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="nondim", default=0.3) - call get_param(param_file, mod, "OBC_RAD_THICK_WT", OBC%gamma_h, & - "The relative weighting for the baroclinic radiation \n"//& - "velocities (or speed of characteristics) at the new \n"//& - "time level (1) or the running mean (0) for thicknesses. \n"//& - "Valid values range from 0 to 1. This is only used if \n"//& - "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="nondim", default=0.2) + call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", OBC%apply_OBC_v_flather_south, & + "Apply a Flather open boundary condition on the southern\n"//& + "side of the global domain", & + default=.false.) + if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v .or. & + OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south .or. & + OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west)) then + ! No open boundaries have been requested + deallocate(OBC) + endif + +end subroutine open_boundary_config + +!> Initialize open boundary control structure +subroutine open_boundary_init(Time, G, param_file, diag, OBC) + type(time_type), target, intent(in) :: Time !< Current model time + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + ! Local variables + + if (.not.associated(OBC)) return + + if ( OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south .or. & + OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west ) then + call get_param(param_file, mod, "OBC_RADIATION_MAX", OBC%rx_max, & + "The maximum magnitude of the baroclinic radiation \n"//& + "velocity (or speed of characteristics). This is only \n"//& + "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="m s-1", default=10.0) + call get_param(param_file, mod, "OBC_RAD_VEL_WT", OBC%gamma_uv, & + "The relative weighting for the baroclinic radiation \n"//& + "velocities (or speed of characteristics) at the new \n"//& + "time level (1) or the running mean (0) for velocities. \n"//& + "Valid values range from 0 to 1. This is only used if \n"//& + "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="nondim", default=0.3) + call get_param(param_file, mod, "OBC_RAD_THICK_WT", OBC%gamma_h, & + "The relative weighting for the baroclinic radiation \n"//& + "velocities (or speed of characteristics) at the new \n"//& + "time level (1) or the running mean (0) for thicknesses. \n"//& + "Valid values range from 0 to 1. This is only used if \n"//& + "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="nondim", default=0.2) + endif id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) end subroutine open_boundary_init +!> Query the state of open boundary module configuration +logical function open_boundary_query(OBC, apply_orig_OBCs, apply_orig_Flather) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + logical, optional, intent(in) :: apply_orig_OBCs !< If present, returns True if APPLY_OBC_U/V was set + logical, optional, intent(in) :: apply_orig_Flather !< If present, returns True if APPLY_OBC_*_FLATHER_* was set + open_boundary_query = .false. + if (.not. associated(OBC)) return + if (present(apply_orig_OBCs)) open_boundary_query = OBC%apply_OBC_u .or. OBC%apply_OBC_v + if (present(apply_orig_Flather)) open_boundary_query = OBC%apply_OBC_v_flather_north .or. & + OBC%apply_OBC_v_flather_south .or. & + OBC%apply_OBC_u_flather_east .or. & + OBC%apply_OBC_u_flather_west +end function open_boundary_query + !> Deallocate open boundary data subroutine open_boundary_end(OBC) type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 1d762e7376..175bc2d076 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -19,6 +19,7 @@ module MOM_fixed_initialization use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : open_boundary_config, open_boundary_query use MOM_string_functions, only : uppercase use user_initialization, only : user_initialize_topography, USER_set_OBC_positions use DOME_initialization, only : DOME_initialize_topography, DOME_set_OBC_positions @@ -55,7 +56,7 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) ! Local character(len=200) :: inputdir ! The directory where NetCDF input files are. character(len=200) :: config - logical :: debug, apply_OBC_u, apply_OBC_v + logical :: debug ! This include declares and sets the variable "version". #include "version_variable.h" @@ -81,15 +82,8 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) ! ==================================================================== ! Determine the position of any open boundaries - call get_param(PF, mod, "APPLY_OBC_U", apply_OBC_u, & - "If true, open boundary conditions may be set at some \n"//& - "u-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(PF, mod, "APPLY_OBC_V", apply_OBC_v, & - "If true, open boundary conditions may be set at some \n"//& - "v-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - if (apply_OBC_u .or. apply_OBC_v) then + call open_boundary_config(G, PF, OBC) + if (open_boundary_query(OBC, apply_orig_OBCs=.true.)) then call get_param(PF, mod, "OBC_CONFIG", config, & "A string that sets how the open boundary conditions are \n"//& " configured: \n"//& diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 455cb559c5..ed0ad2eb98 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -242,16 +242,8 @@ subroutine DOME_set_OBC_positions(G, param_file, OBC) integer :: i, j logical :: any_OBC ! Set to true if any points in this subdomain use OBCs - if (.not.associated(OBC)) allocate(OBC) - - call get_param(param_file, mod, "APPLY_OBC_U", OBC%apply_OBC_u, & - "If true, open boundary conditions may be set at some \n"//& - "u-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V", OBC%apply_OBC_v, & - "If true, open boundary conditions may be set at some \n"//& - "v-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) + if (.not.associated(OBC)) call MOM_error(FATAL, & + "DOME_initialization, DOME_set_OBC_positions: OBC type was not allocated!") any_OBC = .false. if (OBC%apply_OBC_u) then From 59c567640c355279bc33cdde5f5f02c881a7b51c Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 00:17:09 -0400 Subject: [PATCH 08/13] Moved set_Flather_Bdry_Conds() into MOM_open_boundary - The initialization routine set_Flather_Bdry_Conds() now lives in MOM_open_boundary.F90, so that we can later split it into fixed and state related parts. - This is part of a larger re-factor of OBC code. - No answer changes. --- src/core/MOM_open_boundary.F90 | 323 ++++++++++++++- .../MOM_state_initialization.F90 | 386 +----------------- 2 files changed, 326 insertions(+), 383 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c7196fe710..400c414698 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4,9 +4,14 @@ module MOM_open_boundary use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : pass_var, pass_vector +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 +use MOM_file_parser, only : get_param, log_version, param_file_type, log_param use MOM_grid, only : ocean_grid_type +use MOM_io, only : EAST_FACE, NORTH_FACE +use MOM_io, only : slasher, read_data +use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type +use MOM_variables, only : thermo_var_ptrs implicit none ; private @@ -17,6 +22,7 @@ module MOM_open_boundary public open_boundary_query public open_boundary_end public Radiation_Open_Bdry_Conds +public set_Flather_Bdry_Conds integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 integer, parameter, public :: OBC_FLATHER_E = 4, OBC_FLATHER_W = 5 @@ -306,6 +312,321 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & end subroutine Radiation_Open_Bdry_Conds +!> Sets the initial definitions of the characteristic open boundary conditions. +!! \author Mehmet Ilicak +subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) + type(ocean_grid_type), intent(inout) :: G !< Ocean 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(param_file_type), intent(in) :: PF !< Parameter file handle + type(tracer_registry_type), pointer :: tracer_Reg !< Tracer registry + ! Local variables + logical :: read_OBC_eta = .false. + logical :: read_OBC_uv = .false. + logical :: read_OBC_TS = .false. + integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: isd_off, jsd_off + integer :: IsdB, IedB, JsdB, JedB + integer :: east_boundary, west_boundary, north_boundary, south_boundary + character(len=40) :: mod = "set_Flather_Bdry_Conds" ! This subroutine's name. + character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path + + real :: temp_u(G%domain%niglobal+1,G%domain%njglobal) + real :: temp_v(G%domain%niglobal,G%domain%njglobal+1) + + 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 + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + call get_param(PF, mod, "READ_OBC_UV", read_OBC_uv, & + "If true, read the values for the velocity open boundary \n"//& + "conditions from the file specified by OBC_FILE.", & + default=.false.) + call get_param(PF, mod, "READ_OBC_ETA", read_OBC_eta, & + "If true, read the values for the sea surface height \n"//& + "open boundary conditions from the file specified by \n"//& + "OBC_FILE.", default=.false.) + call get_param(PF, mod, "READ_OBC_TS", read_OBC_TS, & + "If true, read the values for the temperature and \n"//& + "salinity open boundary conditions from the file \n"//& + "specified by OBC_FILE.", default=.false.) + if (read_OBC_uv .or. read_OBC_eta .or. read_OBC_TS) then + call get_param(PF, mod, "OBC_FILE", OBC_file, & + "The file from which the appropriate open boundary \n"//& + "condition values are read.", default="MOM_OBC_FILE.nc") + call get_param(PF, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + filename = trim(inputdir)//trim(OBC_file) + call log_param(PF, mod, "INPUTDIR/OBC_FILE", filename) + endif + + if (G%symmetric) then + east_boundary = G%ieg + west_boundary = G%isg-1 + north_boundary = G%jeg + south_boundary = G%jsg-1 + else + ! I am not entirely sure that this works properly. -RWH + east_boundary = G%ieg-1 + west_boundary = G%isg + north_boundary = G%jeg-1 + south_boundary = G%jsg + endif + + if (.not.associated(OBC%OBC_mask_u)) then + allocate(OBC%OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_mask_u(:,:) = .false. + endif + if (.not.associated(OBC%OBC_kind_u)) then + allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE + endif + if (.not.associated(OBC%OBC_mask_v)) then + allocate(OBC%OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC%OBC_mask_v(:,:) = .false. + endif + if (.not.associated(OBC%OBC_kind_v)) then + allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE + endif + + if (.not.associated(OBC%vbt_outer)) then + allocate(OBC%vbt_outer(isd:ied,JsdB:JedB)) ; OBC%vbt_outer(:,:) = 0.0 + endif + + if (.not.associated(OBC%ubt_outer)) then + allocate(OBC%ubt_outer(IsdB:IedB,jsd:jed)) ; OBC%ubt_outer(:,:) = 0.0 + endif + + if (.not.associated(OBC%eta_outer_u)) then + allocate(OBC%eta_outer_u(IsdB:IedB,jsd:jed)) ; OBC%eta_outer_u(:,:) = 0.0 + endif + + if (.not.associated(OBC%eta_outer_v)) then + allocate(OBC%eta_outer_v(isd:ied,JsdB:JedB)) ; OBC%eta_outer_v(:,:) = 0.0 + endif + + if (read_OBC_uv) then + call read_data(filename, 'ubt', OBC%ubt_outer, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + call read_data(filename, 'vbt', OBC%vbt_outer, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + endif + + if (read_OBC_eta) then + call read_data(filename, 'eta_outer_u', OBC%eta_outer_u, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + call read_data(filename, 'eta_outer_v', OBC%eta_outer_v, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + endif + + call pass_vector(OBC%eta_outer_u,OBC%eta_outer_v,G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + call pass_vector(OBC%ubt_outer,OBC%vbt_outer,G%Domain) + + ! This code should be modified to allow OBCs to be applied anywhere. + + if (OBC%apply_OBC_u_flather_east) then + ! Determine where u points are applied at east side + do j=jsd,jed ; do I=IsdB,IedB + if ((I+G%idg_offset) == east_boundary) then !eastern side + if (G%mask2dCu(I,j) > 0.50) then + OBC%OBC_mask_u(I,j) = .true. + OBC%OBC_kind_u(I,j) = OBC_FLATHER_E + if (G%mask2dCv(i+1,J) > 0.50) then + OBC%OBC_mask_v(i+1,J) = .true. + if (OBC%OBC_kind_v(i+1,J) == OBC_NONE) OBC%OBC_kind_v(i+1,J) = OBC_FLATHER_E + endif + if (G%mask2dCv(i+1,J-1) > 0.50) then + OBC%OBC_mask_v(i+1,J-1) = .true. + if (OBC%OBC_kind_v(i+1,J-1) == OBC_NONE) OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER_E + endif + endif + endif + enddo ; enddo + endif + + if (OBC%apply_OBC_u_flather_west) then + ! Determine where u points are applied at west side + do j=jsd,jed ; do I=IsdB,IedB + if ((I+G%idg_offset) == west_boundary) then !western side + if (G%mask2dCu(I,j) > 0.50) then + OBC%OBC_mask_u(I,j) = .true. + OBC%OBC_kind_u(I,j) = OBC_FLATHER_W + if (G%mask2dCv(i,J) > 0.50) then + OBC%OBC_mask_v(i,J) = .true. + if (OBC%OBC_kind_v(i,J) == OBC_NONE) OBC%OBC_kind_v(i,J) = OBC_FLATHER_W + endif + if (G%mask2dCv(i,J-1) > 0.50) then + OBC%OBC_mask_v(i,J-1) = .true. + if (OBC%OBC_kind_v(i,J-1) == OBC_NONE) OBC%OBC_kind_v(i,J-1) = OBC_FLATHER_W + endif + endif + endif + enddo ; enddo + endif + + + if (OBC%apply_OBC_v_flather_north) then + ! Determine where v points are applied at north side + do J=JsdB,JedB ; do i=isd,ied + if ((J+G%jdg_offset) == north_boundary) then !northern side + if (G%mask2dCv(i,J) > 0.50) then + OBC%OBC_mask_v(i,J) = .true. + OBC%OBC_kind_v(i,J) = OBC_FLATHER_N + if (G%mask2dCu(I,j+1) > 0.50) then + OBC%OBC_mask_u(I,j+1) = .true. + if (OBC%OBC_kind_u(I,j+1) == OBC_NONE) OBC%OBC_kind_u(I,j+1) = OBC_FLATHER_N + endif + if (G%mask2dCu(I-1,j+1) > 0.50) then + OBC%OBC_mask_u(I-1,j+1) = .true. + if (OBC%OBC_kind_u(I-1,j+1) == OBC_NONE) OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER_N + endif + endif + endif + enddo ; enddo + endif + + if (OBC%apply_OBC_v_flather_south) then + ! Determine where v points are applied at south side + do J=JsdB,JedB ; do i=isd,ied + if ((J+G%jdg_offset) == south_boundary) then !southern side + if (G%mask2dCv(i,J) > 0.50) then + OBC%OBC_mask_v(i,J) = .true. + OBC%OBC_kind_v(i,J) = OBC_FLATHER_S + if (G%mask2dCu(I,j) > 0.50) then + OBC%OBC_mask_u(I,j) = .true. + if (OBC%OBC_kind_u(I,j) == OBC_NONE) OBC%OBC_kind_u(I,j) = OBC_FLATHER_S + endif + if (G%mask2dCu(I-1,j) > 0.50) then + OBC%OBC_mask_u(I-1,j) = .true. + if (OBC%OBC_kind_u(I-1,j) == OBC_NONE) OBC%OBC_kind_u(I-1,j) = OBC_FLATHER_S + endif + endif + endif + enddo ; enddo + endif + + ! If there are no OBC points on this PE, there is no reason to keep the OBC + ! type, and it could be deallocated. + + + ! Define radiation coefficients r[xy]_old_[uvh] as needed. For now, there are + ! no radiation conditions applied to the thicknesses, since the thicknesses + ! might not be physically motivated. Instead, sponges should be used to + ! enforce the near-boundary layer structure. + if (OBC%apply_OBC_u_flather_west .or. OBC%apply_OBC_u_flather_east) then + allocate(OBC%rx_old_u(IsdB:IedB,jsd:jed,nz)) ; OBC%rx_old_u(:,:,:) = 0.0 + ! allocate(OBC%rx_old_h(Isd:Ied,jsd:jed,nz)) ; OBC%rx_old_h(:,:,:) = 0.0 + endif + if (OBC%apply_OBC_v_flather_south .or. OBC%apply_OBC_v_flather_north) then + allocate(OBC%ry_old_v(isd:ied,JsdB:JedB,nz)) ; OBC%ry_old_v(:,:,:) = 0.0 + ! 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 + allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) ; OBC_T_v(:,:,:) = 0.0 + allocate(OBC_S_v(isd:ied,JsdB:JedB,nz)) ; OBC_S_v(:,:,:) = 0.0 + + if (read_OBC_TS) then + call read_data(filename, 'OBC_T_u', OBC_T_u, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + call read_data(filename, 'OBC_S_u', OBC_S_u, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + + call read_data(filename, 'OBC_T_v', OBC_T_v, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + call read_data(filename, 'OBC_S_v', OBC_S_v, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + else + call pass_var(tv%T, G%Domain) + call pass_var(tv%S, G%Domain) + do k=1,nz ; do j=js,je ; do I=is-1,ie + if (OBC%OBC_mask_u(I,j)) then + if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_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_kind_u(I,j) == OBC_FLATHER_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=js-1,je ; do i=is,ie + if (OBC%OBC_mask_v(i,J)) then + if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_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_kind_v(i,J) == OBC_FLATHER_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 + endif + + 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) + do k=1,nz ; do j=js,je ; do I=is-1,ie + if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) then + tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k) + elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then + tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k) + endif + enddo ; enddo ; enddo + do k=1,nz ; do J=js-1,je ; do i=is,ie + if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) then + tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k) + elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then + tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k) + endif + enddo ; enddo ; enddo + endif + + do k=1,nz ; do j=js-1,je+1 ; do I=is-1,ie+1 + if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) h(i+1,j,k) = h(i,j,k) + if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) h(i,j,k) = h(i+1,j,k) + enddo ; enddo ; enddo + do k=1,nz ; do J=js-1,je+1 ; do i=is-1,ie+1 + if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) h(i,j+1,k) = h(i,j,k) + if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) h(i,j,k) = h(i,j+1,k) + enddo ; enddo ; enddo + +end subroutine set_Flather_Bdry_Conds + !> \namespace mom_open_boundary !! This module implements some aspects of internal open boundary !! conditions in MOM. diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index a846988fde..cf77120ee3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -22,8 +22,8 @@ module MOM_state_initialization use MOM_io, only : slasher, vardesc, write_field use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_open_boundary, only : ocean_OBC_type -use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE, OBC_FLATHER_E, OBC_FLATHER_W -use MOM_open_boundary, only : OBC_FLATHER_N, OBC_FLATHER_S +use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE +use MOM_open_boundary, only : open_boundary_query, set_Flather_Bdry_Conds use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_restart, only : restore_state, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density @@ -142,9 +142,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced ! by a large surface pressure, such as with an ice sheet. logical :: Analytic_FV_PGF, obsol_test - logical :: apply_OBC_u, apply_OBC_v - logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west - logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south logical :: convert type(EOS_type), pointer :: eos => NULL() logical :: debug ! indicates whether to write debugging output @@ -424,15 +421,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & endif ! This subroutine call sets optional open boundary conditions. - call get_param(PF, mod, "APPLY_OBC_U", apply_OBC_u, & - "If true, open boundary conditions may be set at some \n"//& - "u-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(PF, mod, "APPLY_OBC_V", apply_OBC_v, & - "If true, open boundary conditions may be set at some \n"//& - "v-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - if (apply_OBC_u .or. apply_OBC_v) then + if (open_boundary_query(OBC, apply_orig_OBCs=.true.)) then call get_param(PF, mod, "OBC_CONFIG", config, fail_if_missing=.true., do_not_log=.true.) if (trim(config) == "DOME") then call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) @@ -445,19 +434,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & endif endif - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,& - "Apply a Flather open boundary condition on the eastern \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,& - "Apply a Flather open boundary condition on the western \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,& - "Apply a Flather open boundary condition on the northern \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,& - "Apply a Flather open boundary condition on the southern \n"//& - "side of the global domain", default=.false.) - if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) then + if (open_boundary_query(OBC, apply_orig_Flather=.true.)) then call set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) endif @@ -1660,361 +1637,6 @@ subroutine set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tracer_Reg) end subroutine set_Open_Bdry_Conds ! ----------------------------------------------------------------------------- -! ----------------------------------------------------------------------------- -subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) - type(ocean_grid_type), intent(inout) :: G - type(ocean_OBC_type), pointer :: OBC - type(thermo_var_ptrs), intent(inout) :: tv - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h - type(param_file_type), intent(in) :: PF - type(tracer_registry_type), pointer :: tracer_Reg -! This subroutine sets the initial definitions of the characteristic open boundary -! conditions. Written by Mehmet Ilicak - -! Arguments: OBC - This open boundary condition type specifies whether, where, -! and what open boundary conditions are used. -! (out) tv - A structure containing pointers to any available -! thermodynamic fields, including potential temperature and -! salinity or mixed layer density. Absent fields have NULL ptrs. -! (in) G - The ocean's grid structure. -! (in) PF - A structure indicating the open file to parse for -! model parameter values. - - logical :: any_OBC ! Set to true if any points in this subdomain use - ! open boundary conditions. - - logical :: apply_OBC_u_flather_east = .false., apply_OBC_u_flather_west = .false. - logical :: apply_OBC_v_flather_north = .false., apply_OBC_v_flather_south = .false. - logical :: read_OBC_eta = .false. - logical :: read_OBC_uv = .false. - logical :: read_OBC_TS = .false. - - integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz - integer :: isd_off, jsd_off - integer :: IsdB, IedB, JsdB, JedB - integer :: east_boundary, west_boundary, north_boundary, south_boundary - character(len=40) :: mod = "set_Flather_Bdry_Conds" ! This subroutine's name. - character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path - - real :: temp_u(G%domain%niglobal+1,G%domain%njglobal) - real :: temp_v(G%domain%niglobal,G%domain%njglobal+1) - - 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 - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,& - "Apply a Flather open boundary condition on the eastern \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,& - "Apply a Flather open boundary condition on the western \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,& - "Apply a Flather open boundary condition on the northern \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,& - "Apply a Flather open boundary condition on the southern \n"//& - "side of the global domain", default=.false.) - - if (.not.(apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. & - apply_OBC_v_flather_north .or. apply_OBC_v_flather_south)) return - - if (.not.associated(OBC)) allocate(OBC) - - OBC%apply_OBC_u_flather_east = apply_OBC_u_flather_east - OBC%apply_OBC_u_flather_west = apply_OBC_u_flather_west - OBC%apply_OBC_v_flather_north = apply_OBC_v_flather_north - OBC%apply_OBC_v_flather_south = apply_OBC_v_flather_south - - call get_param(PF, mod, "READ_OBC_UV", read_OBC_uv, & - "If true, read the values for the velocity open boundary \n"//& - "conditions from the file specified by OBC_FILE.", & - default=.false.) - call get_param(PF, mod, "READ_OBC_ETA", read_OBC_eta, & - "If true, read the values for the sea surface height \n"//& - "open boundary conditions from the file specified by \n"//& - "OBC_FILE.", default=.false.) - call get_param(PF, mod, "READ_OBC_TS", read_OBC_TS, & - "If true, read the values for the temperature and \n"//& - "salinity open boundary conditions from the file \n"//& - "specified by OBC_FILE.", default=.false.) - if (read_OBC_uv .or. read_OBC_eta .or. read_OBC_TS) then - call get_param(PF, mod, "OBC_FILE", OBC_file, & - "The file from which the appropriate open boundary \n"//& - "condition values are read.", default="MOM_OBC_FILE.nc") - call get_param(PF, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - filename = trim(inputdir)//trim(OBC_file) - call log_param(PF, mod, "INPUTDIR/OBC_FILE", filename) - endif - - if (G%symmetric) then - east_boundary = G%ieg - west_boundary = G%isg-1 - north_boundary = G%jeg - south_boundary = G%jsg-1 - else - ! I am not entirely sure that this works properly. -RWH - east_boundary = G%ieg-1 - west_boundary = G%isg - north_boundary = G%jeg-1 - south_boundary = G%jsg - endif - - if (.not.associated(OBC%OBC_mask_u)) then - allocate(OBC%OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_mask_u(:,:) = .false. - endif - if (.not.associated(OBC%OBC_kind_u)) then - allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE - endif - if (.not.associated(OBC%OBC_mask_v)) then - allocate(OBC%OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC%OBC_mask_v(:,:) = .false. - endif - if (.not.associated(OBC%OBC_kind_v)) then - allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE - endif - - if (.not.associated(OBC%vbt_outer)) then - allocate(OBC%vbt_outer(isd:ied,JsdB:JedB)) ; OBC%vbt_outer(:,:) = 0.0 - endif - - if (.not.associated(OBC%ubt_outer)) then - allocate(OBC%ubt_outer(IsdB:IedB,jsd:jed)) ; OBC%ubt_outer(:,:) = 0.0 - endif - - if (.not.associated(OBC%eta_outer_u)) then - allocate(OBC%eta_outer_u(IsdB:IedB,jsd:jed)) ; OBC%eta_outer_u(:,:) = 0.0 - endif - - if (.not.associated(OBC%eta_outer_v)) then - allocate(OBC%eta_outer_v(isd:ied,JsdB:JedB)) ; OBC%eta_outer_v(:,:) = 0.0 - endif - - if (read_OBC_uv) then - call read_data(filename, 'ubt', OBC%ubt_outer, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - call read_data(filename, 'vbt', OBC%vbt_outer, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - endif - - if (read_OBC_eta) then - call read_data(filename, 'eta_outer_u', OBC%eta_outer_u, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - call read_data(filename, 'eta_outer_v', OBC%eta_outer_v, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - endif - - call pass_vector(OBC%eta_outer_u,OBC%eta_outer_v,G%Domain, To_All+SCALAR_PAIR, CGRID_NE) - call pass_vector(OBC%ubt_outer,OBC%vbt_outer,G%Domain) - - ! This code should be modified to allow OBCs to be applied anywhere. - - if (apply_OBC_u_flather_east) then - ! Determine where u points are applied at east side - do j=jsd,jed ; do I=IsdB,IedB - if ((I+G%idg_offset) == east_boundary) then !eastern side - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_kind_u(I,j) = OBC_FLATHER_E - if (G%mask2dCv(i+1,J) > 0.50) then - OBC%OBC_mask_v(i+1,J) = .true. - if (OBC%OBC_kind_v(i+1,J) == OBC_NONE) OBC%OBC_kind_v(i+1,J) = OBC_FLATHER_E - endif - if (G%mask2dCv(i+1,J-1) > 0.50) then - OBC%OBC_mask_v(i+1,J-1) = .true. - if (OBC%OBC_kind_v(i+1,J-1) == OBC_NONE) OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER_E - endif - endif - endif - enddo ; enddo - endif - - if (apply_OBC_u_flather_west) then - ! Determine where u points are applied at west side - do j=jsd,jed ; do I=IsdB,IedB - if ((I+G%idg_offset) == west_boundary) then !western side - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_kind_u(I,j) = OBC_FLATHER_W - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - if (OBC%OBC_kind_v(i,J) == OBC_NONE) OBC%OBC_kind_v(i,J) = OBC_FLATHER_W - endif - if (G%mask2dCv(i,J-1) > 0.50) then - OBC%OBC_mask_v(i,J-1) = .true. - if (OBC%OBC_kind_v(i,J-1) == OBC_NONE) OBC%OBC_kind_v(i,J-1) = OBC_FLATHER_W - endif - endif - endif - enddo ; enddo - endif - - - if (apply_OBC_v_flather_north) then - ! Determine where v points are applied at north side - do J=JsdB,JedB ; do i=isd,ied - if ((J+G%jdg_offset) == north_boundary) then !northern side - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_kind_v(i,J) = OBC_FLATHER_N - if (G%mask2dCu(I,j+1) > 0.50) then - OBC%OBC_mask_u(I,j+1) = .true. - if (OBC%OBC_kind_u(I,j+1) == OBC_NONE) OBC%OBC_kind_u(I,j+1) = OBC_FLATHER_N - endif - if (G%mask2dCu(I-1,j+1) > 0.50) then - OBC%OBC_mask_u(I-1,j+1) = .true. - if (OBC%OBC_kind_u(I-1,j+1) == OBC_NONE) OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER_N - endif - endif - endif - enddo ; enddo - endif - - if (apply_OBC_v_flather_south) then - ! Determine where v points are applied at south side - do J=JsdB,JedB ; do i=isd,ied - if ((J+G%jdg_offset) == south_boundary) then !southern side - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_kind_v(i,J) = OBC_FLATHER_S - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - if (OBC%OBC_kind_u(I,j) == OBC_NONE) OBC%OBC_kind_u(I,j) = OBC_FLATHER_S - endif - if (G%mask2dCu(I-1,j) > 0.50) then - OBC%OBC_mask_u(I-1,j) = .true. - if (OBC%OBC_kind_u(I-1,j) == OBC_NONE) OBC%OBC_kind_u(I-1,j) = OBC_FLATHER_S - endif - endif - endif - enddo ; enddo - endif - - ! If there are no OBC points on this PE, there is no reason to keep the OBC - ! type, and it could be deallocated. - - - ! Define radiation coefficients r[xy]_old_[uvh] as needed. For now, there are - ! no radiation conditions applied to the thicknesses, since the thicknesses - ! might not be physically motivated. Instead, sponges should be used to - ! enforce the near-boundary layer structure. - if (apply_OBC_u_flather_west .or. apply_OBC_u_flather_east) then - allocate(OBC%rx_old_u(IsdB:IedB,jsd:jed,nz)) ; OBC%rx_old_u(:,:,:) = 0.0 - ! allocate(OBC%rx_old_h(Isd:Ied,jsd:jed,nz)) ; OBC%rx_old_h(:,:,:) = 0.0 - endif - if (apply_OBC_v_flather_south .or. apply_OBC_v_flather_north) then - allocate(OBC%ry_old_v(isd:ied,JsdB:JedB,nz)) ; OBC%ry_old_v(:,:,:) = 0.0 - ! 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 - allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) ; OBC_T_v(:,:,:) = 0.0 - allocate(OBC_S_v(isd:ied,JsdB:JedB,nz)) ; OBC_S_v(:,:,:) = 0.0 - - if (read_OBC_TS) then - call read_data(filename, 'OBC_T_u', OBC_T_u, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - call read_data(filename, 'OBC_S_u', OBC_S_u, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - - call read_data(filename, 'OBC_T_v', OBC_T_v, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - call read_data(filename, 'OBC_S_v', OBC_S_v, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - else - call pass_var(tv%T, G%Domain) - call pass_var(tv%S, G%Domain) - do k=1,nz ; do j=js,je ; do I=is-1,ie - if (OBC%OBC_mask_u(I,j)) then - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_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_kind_u(I,j) == OBC_FLATHER_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=js-1,je ; do i=is,ie - if (OBC%OBC_mask_v(i,J)) then - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_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_kind_v(i,J) == OBC_FLATHER_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 - endif - - 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) - do k=1,nz ; do j=js,je ; do I=is-1,ie - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) then - tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k) - elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then - tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k) - endif - enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je ; do i=is,ie - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) then - tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k) - elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then - tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k) - endif - enddo ; enddo ; enddo - endif - - do k=1,nz ; do j=js-1,je+1 ; do I=is-1,ie+1 - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) h(i+1,j,k) = h(i,j,k) - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) h(i,j,k) = h(i+1,j,k) - enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je+1 ; do i=is-1,ie+1 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) h(i,j+1,k) = h(i,j,k) - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) h(i,j,k) = h(i,j+1,k) - enddo ; enddo ; enddo - -end subroutine set_Flather_Bdry_Conds -! ----------------------------------------------------------------------------- - ! ----------------------------------------------------------------------------- subroutine set_velocity_depth_max(G) type(ocean_grid_type), intent(inout) :: G From 0248eb2db2348fafa84a616d7310ed155b53cb52 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 18:47:38 -0400 Subject: [PATCH 09/13] Separated set_Flather_... into fixed and state parts - The old routine set_Flather_Bdry_Conds used to set both the masks and data, all hard-corded to be along the edges of the domain. This separates the setting of masks and the construction of the boundary data. - This is part of a larger re-factor of OBC code. - No answer changes. --- src/core/MOM_open_boundary.F90 | 218 ++++++++---------- .../MOM_fixed_initialization.F90 | 3 + .../MOM_state_initialization.F90 | 13 +- 3 files changed, 113 insertions(+), 121 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 400c414698..027a5b3474 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -22,7 +22,8 @@ module MOM_open_boundary public open_boundary_query public open_boundary_end public Radiation_Open_Bdry_Conds -public set_Flather_Bdry_Conds +public set_Flather_positions +public set_Flather_data integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 integer, parameter, public :: OBC_FLATHER_E = 4, OBC_FLATHER_W = 5 @@ -312,9 +313,106 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & end subroutine Radiation_Open_Bdry_Conds +!> Sets the domain boundaries as Flather open boundaries using the original +!! Flather run-time logicals +subroutine set_Flather_positions(G, OBC) + type(ocean_grid_type), intent(inout) :: G + type(ocean_OBC_type), pointer :: OBC + ! Local variables + integer :: east_boundary, west_boundary, north_boundary, south_boundary + integer :: i, j + + if (.not.associated(OBC%OBC_mask_u)) then + allocate(OBC%OBC_mask_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_mask_u(:,:) = .false. + endif + if (.not.associated(OBC%OBC_kind_u)) then + allocate(OBC%OBC_kind_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE + endif + if (.not.associated(OBC%OBC_mask_v)) then + allocate(OBC%OBC_mask_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_mask_v(:,:) = .false. + endif + if (.not.associated(OBC%OBC_kind_v)) then + allocate(OBC%OBC_kind_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE + endif + + ! This code should be modified to allow OBCs to be applied anywhere. + + if (G%symmetric) then + east_boundary = G%ieg + west_boundary = G%isg-1 + north_boundary = G%jeg + south_boundary = G%jsg-1 + else + ! I am not entirely sure that this works properly. -RWH + east_boundary = G%ieg-1 + west_boundary = G%isg + north_boundary = G%jeg-1 + south_boundary = G%jsg + endif + + if (OBC%apply_OBC_u_flather_east) then + ! Determine where u points are applied at east side + do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB + if ((I+G%idg_offset) == east_boundary) then !eastern side + OBC%OBC_mask_u(I,j) = .true. + OBC%OBC_kind_u(I,j) = OBC_FLATHER_E + OBC%OBC_mask_v(i+1,J) = .true. + if (OBC%OBC_kind_v(i+1,J) == OBC_NONE) OBC%OBC_kind_v(i+1,J) = OBC_FLATHER_E + OBC%OBC_mask_v(i+1,J-1) = .true. + if (OBC%OBC_kind_v(i+1,J-1) == OBC_NONE) OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER_E + endif + enddo ; enddo + endif + + if (OBC%apply_OBC_u_flather_west) then + ! Determine where u points are applied at west side + do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB + if ((I+G%idg_offset) == west_boundary) then !western side + OBC%OBC_mask_u(I,j) = .true. + OBC%OBC_kind_u(I,j) = OBC_FLATHER_W + OBC%OBC_mask_v(i,J) = .true. + if (OBC%OBC_kind_v(i,J) == OBC_NONE) OBC%OBC_kind_v(i,J) = OBC_FLATHER_W + OBC%OBC_mask_v(i,J-1) = .true. + if (OBC%OBC_kind_v(i,J-1) == OBC_NONE) OBC%OBC_kind_v(i,J-1) = OBC_FLATHER_W + endif + enddo ; enddo + endif + + if (OBC%apply_OBC_v_flather_north) then + ! Determine where v points are applied at north side + do J=G%JsdB,G%JedB ; do i=G%isd,G%ied + if ((J+G%jdg_offset) == north_boundary) then !northern side + OBC%OBC_mask_v(i,J) = .true. + OBC%OBC_kind_v(i,J) = OBC_FLATHER_N + OBC%OBC_mask_u(I,j+1) = .true. + if (OBC%OBC_kind_u(I,j+1) == OBC_NONE) OBC%OBC_kind_u(I,j+1) = OBC_FLATHER_N + OBC%OBC_mask_u(I-1,j+1) = .true. + if (OBC%OBC_kind_u(I-1,j+1) == OBC_NONE) OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER_N + endif + enddo ; enddo + endif + + if (OBC%apply_OBC_v_flather_south) then + ! Determine where v points are applied at south side + do J=G%JsdB,G%JedB ; do i=G%isd,G%ied + if ((J+G%jdg_offset) == south_boundary) then !southern side + OBC%OBC_mask_v(i,J) = .true. + OBC%OBC_kind_v(i,J) = OBC_FLATHER_S + OBC%OBC_mask_u(I,j) = .true. + if (OBC%OBC_kind_u(I,j) == OBC_NONE) OBC%OBC_kind_u(I,j) = OBC_FLATHER_S + OBC%OBC_mask_u(I-1,j) = .true. + if (OBC%OBC_kind_u(I-1,j) == OBC_NONE) OBC%OBC_kind_u(I-1,j) = OBC_FLATHER_S + endif + enddo ; enddo + endif + + ! If there are no OBC points on this PE, there is no reason to keep the OBC + ! type, and it could be deallocated. +end subroutine set_Flather_positions + !> Sets the initial definitions of the characteristic open boundary conditions. !! \author Mehmet Ilicak -subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) +subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary structure type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure @@ -328,7 +426,6 @@ subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz integer :: isd_off, jsd_off integer :: IsdB, IedB, JsdB, JedB - integer :: east_boundary, west_boundary, north_boundary, south_boundary character(len=40) :: mod = "set_Flather_Bdry_Conds" ! This subroutine's name. character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path @@ -367,32 +464,6 @@ subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) call log_param(PF, mod, "INPUTDIR/OBC_FILE", filename) endif - if (G%symmetric) then - east_boundary = G%ieg - west_boundary = G%isg-1 - north_boundary = G%jeg - south_boundary = G%jsg-1 - else - ! I am not entirely sure that this works properly. -RWH - east_boundary = G%ieg-1 - west_boundary = G%isg - north_boundary = G%jeg-1 - south_boundary = G%jsg - endif - - if (.not.associated(OBC%OBC_mask_u)) then - allocate(OBC%OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_mask_u(:,:) = .false. - endif - if (.not.associated(OBC%OBC_kind_u)) then - allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE - endif - if (.not.associated(OBC%OBC_mask_v)) then - allocate(OBC%OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC%OBC_mask_v(:,:) = .false. - endif - if (.not.associated(OBC%OBC_kind_v)) then - allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE - endif - if (.not.associated(OBC%vbt_outer)) then allocate(OBC%vbt_outer(isd:ied,JsdB:JedB)) ; OBC%vbt_outer(:,:) = 0.0 endif @@ -426,93 +497,6 @@ subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) call pass_vector(OBC%eta_outer_u,OBC%eta_outer_v,G%Domain, To_All+SCALAR_PAIR, CGRID_NE) call pass_vector(OBC%ubt_outer,OBC%vbt_outer,G%Domain) - ! This code should be modified to allow OBCs to be applied anywhere. - - if (OBC%apply_OBC_u_flather_east) then - ! Determine where u points are applied at east side - do j=jsd,jed ; do I=IsdB,IedB - if ((I+G%idg_offset) == east_boundary) then !eastern side - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_kind_u(I,j) = OBC_FLATHER_E - if (G%mask2dCv(i+1,J) > 0.50) then - OBC%OBC_mask_v(i+1,J) = .true. - if (OBC%OBC_kind_v(i+1,J) == OBC_NONE) OBC%OBC_kind_v(i+1,J) = OBC_FLATHER_E - endif - if (G%mask2dCv(i+1,J-1) > 0.50) then - OBC%OBC_mask_v(i+1,J-1) = .true. - if (OBC%OBC_kind_v(i+1,J-1) == OBC_NONE) OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER_E - endif - endif - endif - enddo ; enddo - endif - - if (OBC%apply_OBC_u_flather_west) then - ! Determine where u points are applied at west side - do j=jsd,jed ; do I=IsdB,IedB - if ((I+G%idg_offset) == west_boundary) then !western side - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_kind_u(I,j) = OBC_FLATHER_W - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - if (OBC%OBC_kind_v(i,J) == OBC_NONE) OBC%OBC_kind_v(i,J) = OBC_FLATHER_W - endif - if (G%mask2dCv(i,J-1) > 0.50) then - OBC%OBC_mask_v(i,J-1) = .true. - if (OBC%OBC_kind_v(i,J-1) == OBC_NONE) OBC%OBC_kind_v(i,J-1) = OBC_FLATHER_W - endif - endif - endif - enddo ; enddo - endif - - - if (OBC%apply_OBC_v_flather_north) then - ! Determine where v points are applied at north side - do J=JsdB,JedB ; do i=isd,ied - if ((J+G%jdg_offset) == north_boundary) then !northern side - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_kind_v(i,J) = OBC_FLATHER_N - if (G%mask2dCu(I,j+1) > 0.50) then - OBC%OBC_mask_u(I,j+1) = .true. - if (OBC%OBC_kind_u(I,j+1) == OBC_NONE) OBC%OBC_kind_u(I,j+1) = OBC_FLATHER_N - endif - if (G%mask2dCu(I-1,j+1) > 0.50) then - OBC%OBC_mask_u(I-1,j+1) = .true. - if (OBC%OBC_kind_u(I-1,j+1) == OBC_NONE) OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER_N - endif - endif - endif - enddo ; enddo - endif - - if (OBC%apply_OBC_v_flather_south) then - ! Determine where v points are applied at south side - do J=JsdB,JedB ; do i=isd,ied - if ((J+G%jdg_offset) == south_boundary) then !southern side - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_kind_v(i,J) = OBC_FLATHER_S - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - if (OBC%OBC_kind_u(I,j) == OBC_NONE) OBC%OBC_kind_u(I,j) = OBC_FLATHER_S - endif - if (G%mask2dCu(I-1,j) > 0.50) then - OBC%OBC_mask_u(I-1,j) = .true. - if (OBC%OBC_kind_u(I-1,j) == OBC_NONE) OBC%OBC_kind_u(I-1,j) = OBC_FLATHER_S - endif - endif - endif - enddo ; enddo - endif - - ! If there are no OBC points on this PE, there is no reason to keep the OBC - ! type, and it could be deallocated. - - ! Define radiation coefficients r[xy]_old_[uvh] as needed. For now, there are ! no radiation conditions applied to the thicknesses, since the thicknesses ! might not be physically motivated. Instead, sponges should be used to @@ -625,7 +609,7 @@ subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) h(i,j,k) = h(i,j+1,k) enddo ; enddo ; enddo -end subroutine set_Flather_Bdry_Conds +end subroutine set_Flather_data !> \namespace mom_open_boundary !! This module implements some aspects of internal open boundary diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 175bc2d076..e9cd32e0ad 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -20,6 +20,7 @@ module MOM_fixed_initialization use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : open_boundary_config, open_boundary_query +use MOM_open_boundary, only : set_Flather_positions use MOM_string_functions, only : uppercase use user_initialization, only : user_initialize_topography, USER_set_OBC_positions use DOME_initialization, only : DOME_initialize_topography, DOME_set_OBC_positions @@ -99,6 +100,8 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) "The open boundary positions specified by OBC_CONFIG="//& trim(config)//" have not been fully implemented.") end select + elseif (open_boundary_query(OBC, apply_orig_Flather=.true.)) then + call set_Flather_positions(G, OBC) endif ! This call sets seamasks that prohibit flow over any point with ! diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index cf77120ee3..b1e7328f7c 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -23,7 +23,7 @@ module MOM_state_initialization use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE -use MOM_open_boundary, only : open_boundary_query, set_Flather_Bdry_Conds +use MOM_open_boundary, only : open_boundary_query, set_Flather_data, set_Flather_positions use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_restart, only : restore_state, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density @@ -432,10 +432,15 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & "OBC_CONFIG = "//trim(config)//" have not been fully implemented.") call set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg) endif + elseif (open_boundary_query(OBC, apply_orig_Flather=.true.)) then +! call set_Flather_positions(G, OBC) + call set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) endif - - if (open_boundary_query(OBC, apply_orig_Flather=.true.)) then - call set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) + if (debug.and.associated(OBC)) then + call hchksum(G%mask2dT, 'MOM_initialize_state: mask2dT ', G%HI) + call uchksum(G%mask2dCu, 'MOM_initialize_state: mask2dCu ', G%HI) + call vchksum(G%mask2dCv, 'MOM_initialize_state: mask2dCv ', G%HI) + call qchksum(G%mask2dBu, 'MOM_initialize_state: mask2dBu ', G%HI) endif call callTree_leave('MOM_initialize_state()') From 0b8dc7beaecf4b3cb33bc36b01e237a3889125d4 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 20:01:14 -0400 Subject: [PATCH 10/13] New function open_boundary_impose_normal_slope() - initialize_masks() in MOM_grid_initialization.F90 used to read the open-boundary parameters and then adjust the topography on the edge of the domain only. - There is no mention of "Flather" in MOM_grid_initialization.F90 any more. - MOM_fixed_initialization.F90 now calls this new routine prior to calling initialize_masks(). - I had to move a pass_var(G%bathyT) out of initialize_masks() and up to MOM_initialize_fixed() to do this. - This is part of a larger re-factor of OBC code. - No answer changes. --- src/core/MOM_open_boundary.F90 | 35 ++++++++++++ .../MOM_fixed_initialization.F90 | 11 +++- src/initialization/MOM_grid_initialize.F90 | 56 +------------------ 3 files changed, 44 insertions(+), 58 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 027a5b3474..5968bfd924 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -21,6 +21,7 @@ module MOM_open_boundary public open_boundary_init public open_boundary_query public open_boundary_end +public open_boundary_impose_normal_slope public Radiation_Open_Bdry_Conds public set_Flather_positions public set_Flather_data @@ -132,6 +133,14 @@ subroutine open_boundary_config(G, param_file, OBC) "Apply a Flather open boundary condition on the southern\n"//& "side of the global domain", & default=.false.) + + ! Safety check + if ((OBC%apply_OBC_u_flather_west .or. OBC%apply_OBC_v_flather_south) .and. & + .not.G%symmetric ) call MOM_error(FATAL, & + "MOM_open_boundary, open_boundary_config: "//& + "Symmetric memory must be used when APPLY_OBC_U_FLATHER_WEST "//& + "or APPLY_OBC_V_FLATHER_SOUTH is true.") + if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v .or. & OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south .or. & OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west)) then @@ -199,6 +208,32 @@ subroutine open_boundary_end(OBC) deallocate(OBC) end subroutine open_boundary_end +!> Sets the slope of bathymetry normal to an open bounndary to zero. +subroutine open_boundary_impose_normal_slope(OBC, G, depth) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points + ! Local variables + integer :: i, j + + if (.not.associated(OBC)) return + + if (associated(OBC%OBC_kind_u)) then + do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 + if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) depth(i+1,j) = depth(i,j) + if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) depth(i,j) = depth(i+1,j) + enddo ; enddo + endif + + if (associated(OBC%OBC_kind_v)) then + do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied + if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) depth(i,j+1) = depth(i,j) + if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) depth(i,j) = depth(i,j+1) + enddo ; enddo + endif + +end subroutine open_boundary_impose_normal_slope + !> Diagnose radiation conditions at open boundaries subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & h_new, h_old, G) diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index e9cd32e0ad..fd31128063 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -20,7 +20,7 @@ module MOM_fixed_initialization use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : open_boundary_config, open_boundary_query -use MOM_open_boundary, only : set_Flather_positions +use MOM_open_boundary, only : set_Flather_positions, open_boundary_impose_normal_slope use MOM_string_functions, only : uppercase use user_initialization, only : user_initialize_topography, USER_set_OBC_positions use DOME_initialization, only : DOME_initialize_topography, DOME_set_OBC_positions @@ -104,8 +104,13 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) call set_Flather_positions(G, OBC) endif -! This call sets seamasks that prohibit flow over any point with ! -! a bottom that is shallower than min_depth from PF. ! + ! To initialize masks, the bathymetry in halo regions must be filled in + call pass_var(G%bathyT, G%Domain) + + ! Make bathymetry consistent with open boundaries + call open_boundary_impose_normal_slope(OBC, G, G%bathyT) + + ! This call sets masks that prohibit flow over any point interpreted as land call initialize_masks(G, PF) if (debug) then call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 03332babf1..6bb268cd3a 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -1343,8 +1343,6 @@ subroutine initialize_masks(G, PF) real :: Dmin, min_depth, mask_depth integer :: i, j - logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west - logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south character(len=40) :: mod = "MOM_grid_init initialize_masks" call callTree_enter("initialize_masks(), MOM_grid_initialize.F90") @@ -1358,65 +1356,13 @@ subroutine initialize_masks(G, PF) "The depth below which to mask points as land points, for which all\n"//& "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", & units="m", default=-9999.0) - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,& - "Apply a Flather open boundary condition on the eastern \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,& - "Apply a Flather open boundary condition on the western \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,& - "Apply a Flather open boundary condition on the northern \n"//& - "side of the global domain", default=.false.) - call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,& - "Apply a Flather open boundary condition on the southern \n"//& - "side of the global domain", default=.false.) - - if ((apply_OBC_u_flather_west .or. apply_OBC_v_flather_south) .and. & - .not.G%symmetric ) & - call MOM_error(FATAL, "Symmetric memory must be used when "//& - "APPLY_OBC_U_FLATHER_WEST or APPLY_OBC_V_FLATHER_SOUTH is true.") Dmin = min_depth if (mask_depth>=0.) Dmin = mask_depth - call pass_var(G%bathyT, G%Domain) G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0 - ! Extrapolate the bottom depths at any points that are subject to Flather - ! open boundary conditions. This should be generalized for Flather OBCs - ! that are not necessarily at the edges of the domain. - if (apply_OBC_u_flather_west) then - do j=G%jsd,G%jed ; do I=G%isd+1,G%ied - if ((I+G%idg_offset) == G%isg) then - G%bathyT(i-1,j) = G%bathyT(i,j) - endif - enddo; enddo - endif - - if (apply_OBC_u_flather_east) then - do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 - if ((i+G%idg_offset) == G%ieg) then - G%bathyT(i+1,j) = G%bathyT(i,j) - endif - enddo; enddo - endif - - if (apply_OBC_v_flather_north) then - do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied - if ((j+G%jdg_offset) == G%jeg) then - G%bathyT(i,j+1) = G%bathyT(i,j) - endif - enddo; enddo - endif - - if (apply_OBC_v_flather_south) then - do J=G%jsd+1,G%jed ; do i=G%isd,G%ied - if ((J+G%jdg_offset) == G%jsg) then - G%bathyT(i,j-1) = G%bathyT(i,j) - endif - enddo; enddo - endif - + ! Construct the h-point or T-point mask do j=G%jsd,G%jed ; do i=G%isd,G%ied if (G%bathyT(i,j) <= Dmin) then G%mask2dT(i,j) = 0.0 From 643063c9a18c8011a6c4fd5eed2d6793d92d686c Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 20:15:10 -0400 Subject: [PATCH 11/13] Added license line to MOM_open_boundary.F90 - Forgot to add the license line when doxygenizing. - No answer changes. --- src/core/MOM_open_boundary.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 5968bfd924..be61791dc3 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1,6 +1,8 @@ !> Controls where open boundary conditions are applied module MOM_open_boundary +! This file is part of MOM6. See LICENSE.md for the license. + use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : pass_var, pass_vector From 2558954ec1366b2cb320974734840959f4b51cb8 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 21:04:29 -0400 Subject: [PATCH 12/13] Moved call to open_boundary_init from MOM to MOM_state_initialize - The call to open_boundary_init() was immediately after the call to MOM_initialize_state(). It is now immediately above the calls to DOME_set_OBC_data, etc. - No answer changes. --- src/core/MOM.F90 | 3 +-- src/core/MOM_open_boundary.F90 | 4 +--- src/initialization/MOM_state_initialization.F90 | 8 +++++--- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c2156d9cbf..4bdd7bf17d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -98,7 +98,7 @@ module MOM use MOM_mixed_layer_restrat, only : mixedlayer_restrat_register_restarts use MOM_neutral_diffusion, only : neutral_diffusion_CS, neutral_diffusion_diag_init use MOM_obsolete_diagnostics, only : register_obsolete_diagnostics -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_init use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS @@ -1799,7 +1799,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, param_file, & dirs, CS%restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) - call open_boundary_init(Time, G, param_file, diag, CS%OBC) call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index be61791dc3..c530519cc8 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -153,11 +153,9 @@ subroutine open_boundary_config(G, param_file, OBC) end subroutine open_boundary_config !> Initialize open boundary control structure -subroutine open_boundary_init(Time, G, param_file, diag, OBC) - type(time_type), target, intent(in) :: Time !< Current model time +subroutine open_boundary_init(G, param_file, OBC) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle - type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index b1e7328f7c..17290657bd 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -21,7 +21,7 @@ module MOM_state_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field use MOM_io, only : EAST_FACE, NORTH_FACE -use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE use MOM_open_boundary, only : open_boundary_query, set_Flather_data, set_Flather_positions use MOM_grid_initialize, only : initialize_masks, set_grid_metrics @@ -420,7 +420,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & end select endif -! This subroutine call sets optional open boundary conditions. + ! Reads OBC parameters not pertaining to the location of the boundaries + call open_boundary_init(G, PF, OBC) + + ! This is the legacy approach to turning on open boundaries if (open_boundary_query(OBC, apply_orig_OBCs=.true.)) then call get_param(PF, mod, "OBC_CONFIG", config, fail_if_missing=.true., do_not_log=.true.) if (trim(config) == "DOME") then @@ -433,7 +436,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & call set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg) endif elseif (open_boundary_query(OBC, apply_orig_Flather=.true.)) then -! call set_Flather_positions(G, OBC) call set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) endif if (debug.and.associated(OBC)) then From ee365bd6bc0544754a39cbec4f7832e2332025ac Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 21:26:16 -0400 Subject: [PATCH 13/13] Loop bound fix in set_Flather_Bdry_data() - This is a conflict resolution: @khedstrom implemented a loop bound fix in a routine that I had just moved to a different module. The resolution was to re-implemented it in the new location. - No answer changes. --- src/core/MOM_open_boundary.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c530519cc8..b153250248 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -619,14 +619,14 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) 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) - do k=1,nz ; do j=js,je ; do I=is-1,ie + do k=1,nz ; do j=jsd,jed ; do I=isd,ied-1 if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) then tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k) elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k) endif enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je ; do i=is,ie + do k=1,nz ; do J=jsd,jed-1 ; do i=isd,ied if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) then tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k) elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then @@ -635,11 +635,11 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) enddo ; enddo ; enddo endif - do k=1,nz ; do j=js-1,je+1 ; do I=is-1,ie+1 + do k=1,nz ; do j=jsd,jed ; do I=isd,ied-1 if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) h(i+1,j,k) = h(i,j,k) if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) h(i,j,k) = h(i+1,j,k) enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je+1 ; do i=is-1,ie+1 + do k=1,nz ; do J=jsd,jed-1 ; do i=isd,ied if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) h(i,j+1,k) = h(i,j,k) if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) h(i,j,k) = h(i,j+1,k) enddo ; enddo ; enddo