From 0248eb2db2348fafa84a616d7310ed155b53cb52 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 14 Jun 2016 18:47:38 -0400 Subject: [PATCH] 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()')