Skip to content

Commit

Permalink
Some index fixes.
Browse files Browse the repository at this point in the history
  • Loading branch information
@kshedstrom
kshedstrom committed Jan 12, 2017
1 parent a083e65 commit ec30811
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 54 deletions.
70 changes: 33 additions & 37 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,24 +85,21 @@ module MOM_open_boundary
!<wave speed (m -s) at OBC-points.
real, pointer, dimension(:,:) :: Htot=>NULL() !<The total column thickness (m) at OBC-points.
real, pointer, dimension(:,:,:) :: h=>NULL() !<The cell thickness (m) at OBC-points.
real, pointer, dimension(:,:,:) :: un=>NULL() !<The grid-oriented normal layer velocity (m s-1).
real, pointer, dimension(:,:,:) :: unh=>NULL() !<The grid-oriented normal layer transports (m3 s-1).
real, pointer, dimension(:,:) :: unhbt=>NULL() !<The vertically summed normal layer transports (m3 s-1).
real, pointer, dimension(:,:) :: unbt=>NULL() !<The vertically-averaged normal velocity (m s-1).
real, pointer, dimension(:,:,:) :: tangent_vel=>NULL() !<The layer velocity tangential to the OB segment (m s-1).
real, pointer, dimension(:,:,:) :: tangent_trans=>NULL() !<The layer transport tangential to the OB segment (m3 s-1).
real, pointer, dimension(:,:,:) :: tangent_trans=>NULL() !<The layer transport tangential to the OB segment (m3 s-1).
real, pointer, dimension(:,:) :: tangent_vel_bt=>NULL() !<The barotropic velocity tangential to the OB segment (m s-1).
real, pointer, dimension(:,:,:) :: normal_vel=>NULL() !<The layer velocity normal to the OB segment (m s-1).
real, pointer, dimension(:,:,:) :: normal_trans=>NULL() !<The layer transport normal to the OB segment (m3 s-1).
real, pointer, dimension(:,:) :: normal_vel_bt=>NULL() !<The barotropic velocity normal to the OB segment (m s-1).
real, pointer, dimension(:,:,:) :: normal_vel=>NULL() !<The layer velocity normal to the OB segment (m s-1).
real, pointer, dimension(:,:,:) :: normal_trans=>NULL() !<The layer transport normal to the OB segment (m3 s-1).
real, pointer, dimension(:,:) :: normal_vel_bt=>NULL() !<The barotropic velocity normal to the OB segment (m s-1).
real, pointer, dimension(:,:) :: normal_trans_bt=>NULL()!<The barotropic transport normal to the OB segment (m3 s-1).
real, pointer, dimension(:,:) :: eta=>NULL() !<The sea-surface elevation along the segment (m).
type(hor_index_type) :: HI !< Horizontal index ranges
end type OBC_segment_type

!> Open-boundary data
type, public :: ocean_OBC_type
integer :: number_of_segments = 0 !< The number of open-boundary segments.
integer :: ke = 0 !< The number of model layers
integer :: ke = 0 !< The number of model layers
logical :: Flather_u_BCs_exist_globally = .false. !< True if any zonal velocity points
!! in the global domain use Flather BCs.
logical :: Flather_v_BCs_exist_globally = .false. !< True if any meridional velocity points
Expand Down Expand Up @@ -430,12 +427,14 @@ subroutine setup_segment_indices(G, seg, Is_obc, Ie_obc, Js_obc, Je_obc)
! Global space I*_obc but sorted
seg%HI%IsgB = Isg ; seg%HI%IegB = Ieg
seg%HI%isg = Isg+1 ; seg%HI%ieg = Ieg
seg%HI%JsgB = Isg ; seg%HI%JegB = Jeg
seg%HI%jsg = Jsg+1 ; seg%HI%jeg = Jeg
seg%HI%JsgB = Jsg ; seg%HI%JegB = Jeg
seg%HI%jsg = Jsg+1 ; seg%HI%Jeg = Jeg

! Move into local index space
Isg = Isg - G%idg_offset
Jsg = Jsg - G%jdg_offset
Ieg = Ieg - G%idg_offset
Jeg = Jeg - G%jdg_offset

! This is the i-extent of the segment on this PE.
! The values are nonsence if the segment is not on this PE.
Expand Down Expand Up @@ -476,7 +475,7 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg)
call parse_segment_str(G%ieg, G%jeg, segment_str, I_obc, Js_obc, Je_obc, action_str )

call setup_segment_indices(G, OBC%OBC_segment_number(l_seg),I_obc,I_obc,Js_obc,Je_obc)
call allocate_OBC_segment_data(G, OBC%OBC_segment_number(l_seg))
call allocate_OBC_segment_data(OBC, OBC%OBC_segment_number(l_seg))

I_obc = I_obc - G%idg_offset ! Convert to local tile indices on this tile
Js_obc = Js_obc - G%jdg_offset ! Convert to local tile indices on this tile
Expand Down Expand Up @@ -602,7 +601,7 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg)
call parse_segment_str(G%ieg, G%jeg, segment_str, J_obc, Is_obc, Ie_obc, action_str )

call setup_segment_indices(G, OBC%OBC_segment_number(l_seg),Is_obc,Ie_obc,J_obc,J_obc)
call allocate_OBC_segment_data(G, OBC%OBC_segment_number(l_seg))
call allocate_OBC_segment_data(OBC, OBC%OBC_segment_number(l_seg))

J_obc = J_obc - G%jdg_offset ! Convert to local tile indices on this tile
Is_obc = Is_obc - G%idg_offset ! Convert to local tile indices on this tile
Expand Down Expand Up @@ -1832,8 +1831,8 @@ subroutine fill_OBC_halos(G, GV, OBC, tv, h, tracer_Reg)
end subroutine fill_OBC_halos

!> Allocate segment data fields
subroutine allocate_OBC_segment_data(G, segment)
type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure
subroutine allocate_OBC_segment_data(OBC, segment)
type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
type(OBC_segment_type), intent(inout) :: segment !< Open boundary segment
! Local variables
integer :: isd, ied, jsd, jed
Expand All @@ -1845,23 +1844,20 @@ subroutine allocate_OBC_segment_data(G, segment)
IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB
JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB

if (.not. associated(seg)) return

if (.not. segment%on_pe) return ! continue to next segment if not in computational domain

if (.not. ASSOCIATED(segment%Cg)) then ! finishing allocating storage for segments
allocate(segment%Cg(isd:ied,jsd:jed)); segment%Cg(:,:)=0.
allocate(segment%Htot(isd:ied,jsd:jed)); segment%Htot(:,:)=0.0
allocate(segment%h(isd:ied,jsd:jed,G%ke)); segment%h(:,:,:)=0.0
allocate(segment%un(isd:ied,jsd:jed,G%ke)); segment%un(:,:,:)=0.0
allocate(segment%unh(isd:ied,jsd:jed,G%ke)); segment%unh(:,:,:)=0.0
allocate(segment%unhbt(isd:ied,jsd:jed)); segment%unhbt(:,:)=0.0
allocate(segment%unbt(isd:ied,jsd:jed)); segment%unbt(:,:)=0.0
allocate(segment%tangent_vel(IsdB:IedB,JsdB:JedB,G%ke)); segment%tangent_vel(:,:,:)=0.0
allocate(segment%tangent_vel_bt(IsdB:IedB,JsdB:JedB)); segment%tangent_vel_bt(:,:)=0.0
allocate(segment%eta(isd:ied,jsd:jed)); segment%eta(:,:)=0.0
endif
enddo
if (.not. segment%on_pe) return

if (.not. ASSOCIATED(segment%Cg)) then ! finishing allocating storage for segments
allocate(segment%Cg(isd:ied,jsd:jed)); segment%Cg(:,:)=0.
allocate(segment%Htot(isd:ied,jsd:jed)); segment%Htot(:,:)=0.0
allocate(segment%h(isd:ied,jsd:jed,OBC%ke)); segment%h(:,:,:)=0.0
allocate(segment%normal_vel(isd:ied,jsd:jed,OBC%ke)); segment%normal_vel(:,:,:)=0.0
allocate(segment%normal_trans(isd:ied,jsd:jed,OBC%ke)); segment%normal_trans(:,:,:)=0.0
allocate(segment%normal_vel_bt(isd:ied,jsd:jed)); segment%normal_vel_bt(:,:)=0.0
allocate(segment%normal_trans_bt(isd:ied,jsd:jed)); segment%normal_trans_bt(:,:)=0.0
allocate(segment%tangent_vel(IsdB:IedB,JsdB:JedB,OBC%ke));segment%tangent_vel(:,:,:)=0.0
allocate(segment%tangent_vel_bt(IsdB:IedB,JsdB:JedB)); segment%tangent_vel_bt(:,:)=0.0
allocate(segment%eta(isd:ied,jsd:jed)); segment%eta(:,:)=0.0
endif
end subroutine allocate_OBC_segment_data

subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time)
Expand Down Expand Up @@ -2060,13 +2056,13 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time)
segment%direction == OBC_DIRECTION_S))) then
do j=js_obc,je_obc
do i=is_obc,ie_obc
segment%unhbt(i,j) = 0.0
segment%normal_trans_bt(i,j) = 0.0
do k=1,G%ke
segment%un(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
segment%unh(i,j,k) = segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k)
segment%unhbt(i,j)= segment%unhbt(i,j)+segment%unh(i,j,k)
segment%normal_vel(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
segment%normal_trans(i,j,k) = segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k)
segment%normal_trans_bt(i,j)= segment%normal_trans_bt(i,j)+segment%normal_trans(i,j,k)
enddo
segment%unbt(i,j) = segment%unhbt(i,j)/max(segment%Htot(i,j),1.e-12)
segment%normal_vel_bt(i,j) = segment%normal_trans_bt(i,j)/max(segment%Htot(i,j),1.e-12)
enddo
enddo
endif
Expand Down
19 changes: 9 additions & 10 deletions src/user/DOME_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ module DOME_initialization
use MOM_get_input, only : directories
use MOM_grid, only : ocean_grid_type
use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE
use MOM_open_boundary, only : OBC_segment_type
use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
Expand Down Expand Up @@ -270,9 +271,6 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg)
integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz
integer :: IsdB, IedB, JsdB, JedB
type(OBC_segment_type), pointer :: segment
integer :: ni_seg, nj_seg ! number of src gridpoints along the segments
integer :: i2, j2 ! indices for referencing local domain array
integer :: ishift, jshift ! offsets for staggered locations

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
Expand All @@ -297,7 +295,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg)
return !!! Need a better error message here
endif
segment => OBC%OBC_segment_number(1)
if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain
if (.not. segment%on_pe) return

do k=1,nz
rst = -1.0
Expand Down Expand Up @@ -332,12 +330,13 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg)
! New way
isd = segment%HI%isd ; ied = segment%HI%ied
JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB
do J=JsdB,JedB ; do i=isd,ied
lon_im1 = 2.0*G%geoLonCv(i,J) - G%geoLonBu(I,J)
segment%unh(i,J,k) = tr_k * (exp(-2.0*(lon_im1 - 1000.0)/Def_Rad) -&
exp(-2.0*(G%geoLonBu(I,J) - 1000.0)/Def_Rad))
segment%un(i,J,k) = v_k * exp(-2.0*(G%geoLonCv(i,J) - 1000.0)/Def_Rad)
enddo ; enddo
print *, 'DOME segment indices', isd, ied, JsdB, JedB
! do J=JsdB,JedB ; do i=isd,ied
! lon_im1 = 2.0*G%geoLonCv(i,J) - G%geoLonBu(I,J)
! segment%normal_trans(i,J,k) = tr_k * (exp(-2.0*(lon_im1 - 1000.0)/Def_Rad) -&
! exp(-2.0*(G%geoLonBu(I,J) - 1000.0)/Def_Rad))
! segment%normal_vel(i,J,k) = v_k * exp(-2.0*(G%geoLonCv(i,J) - 1000.0)/Def_Rad)
! enddo ; enddo
enddo

! The inflow values of temperature and salinity also need to be set here if
Expand Down
16 changes: 9 additions & 7 deletions src/user/tidal_bay_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module tidal_bay_initialization
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE
use MOM_open_boundary, only : OBC_segment_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_time_manager, only : time_type, set_time, time_type_to_real

Expand All @@ -49,8 +50,9 @@ subroutine tidal_bay_set_OBC_data(OBC, G, h, Time)
real :: my_area, my_flux
real :: PI
character(len=40) :: mod = "tidal_bay_set_OBC_data" ! This subroutine's name.
integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz
integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n
integer :: IsdB, IedB, JsdB, JedB
type(OBC_segment_type), pointer :: segment

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Expand Down Expand Up @@ -90,15 +92,15 @@ subroutine tidal_bay_set_OBC_data(OBC, G, h, Time)
enddo ; enddo

! New way
do n = 1, OBC%number_of_segments
segment => OBC%OBC_segment_number(n)
! do n = 1, OBC%number_of_segments
! segment => OBC%OBC_segment_number(n)

if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain
! if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain

segment%unbt(:,:) = my_flux/my_area
segment%eta(:,:) = cff
! segment%normal_vel_bt(:,:) = my_flux/my_area
! segment%eta(:,:) = cff

enddo ! end segment loop
! enddo ! end segment loop

end subroutine tidal_bay_set_OBC_data

Expand Down

0 comments on commit ec30811

Please sign in to comment.