block_control/block_control.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup block_control_mod block_control_mod
  20 !> @ingroup block_control
  21 !> @brief Routines for "blocks" used for  OpenMP threading of column-based
  22 !!        calculations
  23
  24 module block_control_mod
  25
  26 use mpp_mod,         only: mpp_error, NOTE, WARNING, FATAL
  27 use mpp_domains_mod, only: mpp_compute_extent
  28 implicit none
  29
  30 public block_control_type
  31
  32 !> Type to dereference packed index from global index.
  33 !> @ingroup block_control_mod
  34 type :: ix_type
  35   integer, dimension(:,:), allocatable :: ix
  36 end type ix_type
  37
  38 !> Type to dereference packed index from global indices.
  39 !> @ingroup block_control_mod
  40 type :: pk_type
  41   integer, dimension(:), allocatable :: ii
  42   integer, dimension(:), allocatable :: jj
  43 end type pk_type
  44
  45 !> @brief Block data and extents for OpenMP threading of column-based calculations
  46 !> @ingroup block_control_mod
  47 type :: block_control_type
  48   integer :: nx_block, ny_block  !< blocking factor using mpp-style decomposition
  49   integer :: nblks               !< number of blocks cover MPI domain
  50   integer :: isc, iec, jsc, jec  !< MPI domain global extents
  51   integer :: npz                 !< vertical extent
  52   integer, dimension(:),        allocatable :: ibs  , &  !< block extents for mpp-style
  53                                                ibe  , &  !! decompositions
  54                                                jbs  , &
  55                                                jbe
  56   type(ix_type), dimension(:),  allocatable :: ix    !< dereference packed index from global index
  57   !--- packed blocking fields
  58   integer, dimension(:),        allocatable :: blksz !< number of points in each individual block
  59                                                             !! blocks are not required to be uniforom in size
  60   integer, dimension(:,:),      allocatable :: blkno !< dereference block number using global indices
  61   integer, dimension(:,:),      allocatable :: ixp   !< dereference packed index from global indices
  62                                                             !! must be used in conjuction with blkno
  63   type(pk_type), dimension(:),  allocatable :: index !< dereference global indices from
  64                                                             !! block/ixp combo
  65 end type block_control_type
  66
  67 !> @addtogroup block_control_mod
  68 !> @{
  69
  70 public :: define_blocks, define_blocks_packed
  71
  72 contains
  73
  74 !###############################################################################
  75 !> @brief Sets up "blocks" used for OpenMP threading of column-based
  76 !!        calculations using rad_n[x/y]xblock from coupler_nml
  77 !!
  78   subroutine define_blocks (component, Block, isc, iec, jsc, jec, kpts, &
  79                             nx_block, ny_block, message)
  80     character(len=*),         intent(in)    :: component !< Component name string
  81     type(block_control_type), intent(inout) :: Block !< Returns instantiated @ref block_control_type
  82     integer,                  intent(in)    :: isc, iec, jsc, jec, kpts
  83     integer,                  intent(in)    :: nx_block, ny_block
  84     logical,                  intent(inout) :: message !< flag for outputting debug message
  85
  86 !-------------------------------------------------------------------------------
  87 ! Local variables:
  88 !       blocks
  89 !       i1
  90 !       i2
  91 !       j1
  92 !       j2
  93 !       text
  94 !       i
  95 !       j
  96 !       nblks
  97 !       ix
  98 !       ii
  99 !       jj
 100 !-------------------------------------------------------------------------------
 101
 102     integer :: blocks
 103     integer, dimension(nx_block) :: i1, i2
 104     integer, dimension(ny_block) :: j1, j2
 105     character(len=256) :: text
 106     integer :: i, j, nblks, ix, ii, jj
 107
 108     if (message) then
 109       if ((mod(iec-isc+1,nx_block) .ne. 0) .or. (mod(jec-jsc+1,ny_block) .ne. 0)) then
 110         write( text,'(a,a,2i4,a,2i4,a)' ) trim(component),'define_blocks: domain (',&
 111              (iec-isc+1), (jec-jsc+1),') is not an even divisor with definition (',&
 112              nx_block, ny_block,') - blocks will not be uniform'
 113         call mpp_error (WARNING, trim(text))
 114       endif
 115       message = .false.
 116     endif
 117
 118 !--- set up blocks
 119     if (iec-isc+1 .lt. nx_block) &
 120         call mpp_error(FATAL, 'block_control: number of '//trim(component)//' nxblocks .gt. &
 121                              &number of elements in MPI-domain size')
 122     if (jec-jsc+1 .lt. ny_block) &
 123         call mpp_error(FATAL, 'block_control: number of '//trim(component)//' nyblocks .gt. &
 124                              &number of elements in MPI-domain size')
 125     call mpp_compute_extent(isc,iec,nx_block,i1,i2)
 126     call mpp_compute_extent(jsc,jec,ny_block,j1,j2)
 127
 128     nblks = nx_block*ny_block
 129     Block%isc = isc
 130     Block%iec = iec
 131     Block%jsc = jsc
 132     Block%jec = jec
 133     Block%npz = kpts
 134     Block%nx_block = nx_block
 135     Block%ny_block = ny_block
 136     Block%nblks = nblks
 137
 138     if (.not.allocated(Block%ibs)) &
 139          allocate (Block%ibs(nblks), &
 140                    Block%ibe(nblks), &
 141                    Block%jbs(nblks), &
 142                    Block%jbe(nblks), &
 143                    Block%ix(nblks) )
 144
 145     blocks=0
 146     do j = 1, ny_block
 147       do i = 1, nx_block
 148         blocks = blocks + 1
 149         Block%ibs(blocks) = i1(i)
 150         Block%jbs(blocks) = j1(j)
 151         Block%ibe(blocks) = i2(i)
 152         Block%jbe(blocks) = j2(j)
 153         allocate(Block%ix(blocks)%ix(i1(i):i2(i),j1(j):j2(j)) )
 154         ix = 0
 155         do jj = j1(j), j2(j)
 156           do ii = i1(i), i2(i)
 157             ix = ix+1
 158             Block%ix(blocks)%ix(ii,jj) = ix
 159           enddo
 160         enddo
 161       enddo
 162     enddo
 163
 164   end subroutine define_blocks
 165
 166
 167
 168 !###############################################################################
 169 !> @brief Creates and populates a data type which is used for defining the
 170 !!        sub-blocks of the MPI-domain to enhance OpenMP and memory performance.
 171 !!        Uses a packed concept.
 172 !!
 173   subroutine define_blocks_packed (component, Block, isc, iec, jsc, jec, &
 174                                    kpts, blksz, message)
 175     character(len=*),         intent(in)    :: component !< Component name string
 176     type(block_control_type), intent(inout) :: Block !< Returns instantiated @ref block_control_type
 177     integer,                  intent(in)    :: isc, iec, jsc, jec, kpts
 178     integer,                  intent(inout) :: blksz !< block size
 179     logical,                  intent(inout) :: message !< flag for outputting debug message
 180
 181 !-------------------------------------------------------------------------------
 182 ! Local variables:
 183 !       nblks
 184 !       lblksz
 185 !       tot_pts
 186 !       nb
 187 !       ix
 188 !       ii
 189 !       jj
 190 !       text
 191 !-------------------------------------------------------------------------------
 192
 193     integer :: nblks, lblksz, tot_pts, nb, ix, ii, jj
 194     character(len=256) :: text
 195
 196     tot_pts = (iec - isc + 1) * (jec - jsc + 1)
 197     if (blksz < 0) then
 198       nblks = 1
 199       blksz = tot_pts
 200     else
 201       nblks = tot_pts/blksz
 202       if (mod(tot_pts,blksz) .gt. 0) then
 203         nblks = nblks + 1
 204       endif
 205     endif
 206
 207     if (message) then
 208       if (mod(tot_pts,blksz) .ne. 0) then
 209         write( text,'(a,a,2i4,a,i4,a,i4)' ) trim(component),'define_blocks_packed: domain (',&
 210              (iec-isc+1), (jec-jsc+1),') is not an even divisor with definition (',&
 211              blksz,') - blocks will not be uniform with a remainder of ',mod(tot_pts,blksz)
 212         call mpp_error (WARNING, trim(text))
 213       endif
 214       message = .false.
 215     endif
 216
 217     Block%isc   = isc
 218     Block%iec   = iec
 219     Block%jsc   = jsc
 220     Block%jec   = jec
 221     Block%npz   = kpts
 222     Block%nblks = nblks
 223     if (.not. allocated(Block%blksz)) &
 224       allocate (Block%blksz(nblks), &
 225                 Block%index(nblks), &
 226                 Block%blkno(isc:iec,jsc:jec), &
 227                 Block%ixp(isc:iec,jsc:jec))
 228
 229 !--- set up blocks
 230     do nb = 1, nblks
 231       lblksz = blksz
 232       if (nb .EQ. nblks) lblksz = tot_pts - (nb-1) * blksz
 233       Block%blksz(nb) = lblksz
 234       allocate (Block%index(nb)%ii(lblksz), &
 235                 Block%index(nb)%jj(lblksz))
 236     enddo
 237
 238 !--- set up packed indices
 239     nb = 1
 240     ix = 0
 241     do jj = jsc, jec
 242       do ii = isc, iec
 243         ix = ix + 1
 244         if (ix .GT. blksz) then
 245           ix = 1
 246           nb = nb + 1
 247         endif
 248         Block%ixp(ii,jj) = ix
 249         Block%blkno(ii,jj) = nb
 250         Block%index(nb)%ii(ix) = ii
 251         Block%index(nb)%jj(ix) = jj
 252       enddo
 253     enddo
 254
 255   end subroutine define_blocks_packed
 256
 257 end module block_control_mod
 258 !> @}
 259 ! close documentation grouping