module general_sub2grid_mod
!$$$ module documentation block
! . . . .
! module: generic_sub2grid_mod generalized sub2grid and grid2sub style routines
! prgmmr: parrish org: np22 date: 2010-02-05
!
! abstract: This module contains generalized sub2grid and grid2sub like routines
! which are largely independent of other gsi routines/modules.
! This has been created first to allow easier introduction of dual
! resolution capability for the hybrid ensemble option. But
! it may be used eventually to replace most of the specialized
! sub2grid/grid2sub routines.
! NOTE: Since this will initially be used only for ensemble space
! where it is not necessary to have haloes on subdomains, the
! routine general_grid2sub_ strips off the haloes from the output,
! while still including them internally. They are included internally
! since this version still uses the same computation of mpi_alltoallv
! arguments used for existing grid2sub routine. Also, the input for
! general_sub2grid_ has no halo.
! The initial use for this module will be to manipulate fields of ensemble
! perturbations as required when running GSI with the hybrid ensemble mode
! turned on. In this case, the haloes are not required.
! To make this a more useful generalized code, later add option of halo size,
! from 0 to any value.
!
! program history log:
! 2010-02-05 parrish, initial documentation
! 2010-03-02 parrish - restore halo to size 1 and duplicate what is in current sub2grid/grid2sub.
! also, make sure that periodic and periodic_s are properly specified.
! there is a bug in existing use of periodic when running in regional mode.
! periodic should always be false when regional=.true. this has been
! corrected in this version.
! 2010-03-11 parrish - add parameter kend_alloc to type sub2grid_info. this is to fix problem
! of allocating arrays when kend_loc = kbegin_loc-1, which happens for
! processors not involved with sub2grid/grid2sub.
! to fix this, use kend_alloc = max(kend_loc,kbegin_loc) for allocation.
!
! subroutines included:
! sub general_sub2grid_r_single - convert from subdomains to grid for real single precision (4 byte)
! sub general_grid2sub_r_single - convert from grid to subdomains for real single precision (4 byte)
! Variable Definitions:
! def sub2grid_info - contains all information needed for general_sub2grid and general_grid2sub
!
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
implicit none
! set default to private
private
! set subroutines to public
public :: general_sub2grid_r_single
public :: general_grid2sub_r_single
public :: general_sub2grid_r_double
public :: general_grid2sub_r_double
public :: general_sub2grid_create_info
public :: general_sube2suba_r_double
public :: general_sube2suba_r_double_ad
public :: general_suba2sube_r_double
! set passed variables to public
public :: sub2grid_info
type sub2grid_info
integer(i_kind) inner_vars ! number of inner-most loop variables
integer(i_kind) lat1 ! no. of lats on subdomain (no buffer)
integer(i_kind) lon1 ! no. of lons on subdomain (no buffer)
integer(i_kind) lat2 ! no. of lats on subdomain (buffer)
integer(i_kind) lon2 ! no. of lons on subdomain (buffer)
integer(i_kind) latlon11 ! no. of points on subdomain (including buffer)
integer(i_kind) latlon1n ! latlon11*nsig
integer(i_kind) nlat ! no. of latitudes
integer(i_kind) nlon ! no. of longitudes
integer(i_kind) nsig ! no. of vertical levels
integer(i_kind) num_fields ! total number of fields/levels
integer(i_kind) iglobal ! number of horizontal points on global grid
integer(i_kind) itotsub ! number of horizontal points of all subdomains combined
integer(i_kind) kbegin_loc ! starting slab index for local processor
integer(i_kind) kend_loc ! ending slab index for local processor
integer(i_kind) kend_alloc ! kend_loc can = kbegin_loc - 1, for a processor not involved.
! this causes problems with array allocation:
! to correct this, use kend_alloc=max(kend_loc,kbegin_loc)
integer(i_kind) npe ! total number of processors
integer(i_kind) mype ! local processor
logical periodic ! logical flag for periodic e/w domains
logical,pointer :: periodic_s(:) => NULL() ! logical flag for periodic e/w subdomain (all tasks)
logical,pointer :: vector(:) => NULL() ! logical flag, true for vector variables
integer(i_kind),pointer :: ilat1(:) => NULL() ! no. of lats for each subdomain (no buffer)
integer(i_kind),pointer :: jlon1(:) => NULL() ! no. of lons for each subdomain (no buffer)
integer(i_kind),pointer :: istart(:) => NULL() ! start lat of the whole array on each pe
integer(i_kind),pointer :: jstart(:) => NULL() ! start lon of the whole array on each pe
integer(i_kind),pointer :: recvcounts(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: displs_g(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: rdispls(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: sendcounts(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: sdispls(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: ijn(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: ltosj(:) => NULL() ! lat index for reordering slab
integer(i_kind),pointer :: ltosi(:) => NULL() ! lon index for reordering slab
integer(i_kind),pointer :: recvcounts_s(:)=> NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: irc_s(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: ird_s(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: displs_s(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: rdispls_s(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: sendcounts_s(:)=> NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: sdispls_s(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: ijn_s(:) => NULL() ! for mpi_alltoallv (sub2grid)
integer(i_kind),pointer :: ltosj_s(:) => NULL() ! lat index for reordering slab
integer(i_kind),pointer :: ltosi_s(:) => NULL() ! lon index for reordering slab
integer(i_kind),pointer :: kbegin(:) => NULL() ! starting slab index for each processor
integer(i_kind),pointer :: kend(:) => NULL() ! ending slab index for each processor
logical:: lallocated = .false.
end type sub2grid_info
! other declarations ...
contains
subroutine general_sub2grid_create_info(s,inner_vars,nlat,nlon,nsig,num_fields,regional,vector)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sub2grid_create_info populate info variable s
! prgmmr: parrish org: np22 date: 2010-02-12
!
! abstract: given dimensions of horizontal domain and various other
! information, obtain all required constants to allow
! use of general_sub2grid_ and general_grid2sub_ and store them
! in structure variable s.
!
! program history log:
! 2010-02-11 parrish, initial documentation
! 2010-03-02 parrish - add regional flag to input. if regional=.true.,
! then periodic, periodic_s=.false. always. this corrects a bug
! in existing code. (never a problem, except when npe=1).
!
! input argument list:
! s - structure variable, waiting for all necessary information for
! use with general_sub2grid and general_grid2sub.
! inner_vars - inner index, reserved for eventually putting all ensemble members
! on 1st (most rapidly varying) array index.
! nlat - number of horizontal grid points in "latitude" direction
! nlon - number of horizontal grid points in "longitude"
! nsig - number of vertical levels for 1 3d variable.
! num_fields - total number of 2d fields to be processed.
! vector - optional logical array of length num_fields, set to true for
! each field which will be a vector component.
! if not present, s%vector = .false.
! regional - if true, then no periodicity in "longitude" direction
!
! output argument list:
! s - structure variable, contains all necessary information for
! moving this set of subdomain variables sub_vars to
! the corresponding set of full horizontal grid variables.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_single
use constants, only: izero,ione
use mpimod, only: mpi_comm_world
implicit none
type(sub2grid_info),intent(inout) :: s
integer(i_kind), intent(in ) :: inner_vars,nlat,nlon,nsig,num_fields
logical, intent(in ) :: regional
logical,optional, intent(in ) :: vector(num_fields)
integer(i_kind) i,ierror,j,k,num_loc_groups,nextra,mm1,n,ns
integer(i_kind),allocatable:: isc_g(:),isd_g(:)
call mpi_comm_size(mpi_comm_world,s%npe,ierror)
call mpi_comm_rank(mpi_comm_world,s%mype,ierror)
s%inner_vars=inner_vars
s%nlat=nlat
s%nlon=nlon
s%iglobal=nlat*nlon
s%nsig=nsig
s%num_fields=num_fields
if(s%lallocated) then
deallocate(s%periodic_s,s%ilat1,s%istart,s%jlon1,s%jstart,s%kbegin,s%kend,s%ijn)
deallocate(s%sendcounts,s%sdispls,s%recvcounts,s%rdispls)
deallocate(s%sendcounts_s,s%sdispls_s,s%recvcounts_s,s%rdispls_s)
deallocate(s%ltosi,s%ltosj,s%ltosi_s,s%ltosj_s)
deallocate(s%displs_s,s%irc_s,s%ird_s)
deallocate(s%displs_g)
deallocate(s%vector)
s%lallocated=.false.
end if
allocate(s%periodic_s(s%npe),s%jstart(s%npe),s%istart(s%npe),s%ilat1(s%npe),s%jlon1(s%npe))
allocate(s%ijn(s%npe),s%ijn_s(s%npe))
allocate(s%vector(num_fields))
if(present(vector)) then
s%vector=vector
else
s%vector=.false.
end if
! first determine subdomains
call general_deter_subdomain_nolayout(s%npe,s%mype,s%nlat,s%nlon,regional, &
s%periodic,s%periodic_s,s%lon1,s%lon2,s%lat1,s%lat2,s%ilat1,s%istart,s%jlon1,s%jstart)
s%latlon11=s%lat2*s%lon2
s%latlon1n=s%latlon11*s%nsig
allocate(isc_g(s%npe),isd_g(s%npe),s%displs_g(s%npe),s%displs_s(s%npe),s%ird_s(s%npe),s%irc_s(s%npe))
s%ijn=s%ilat1*s%jlon1
s%ijn_s=(s%ilat1+2_i_kind)*(s%jlon1+2_i_kind)
mm1=s%mype+ione
do i=1,s%npe
s%irc_s(i)=s%ijn_s(mm1)
isc_g(i)=s%ijn(mm1)
end do
! obtain ltosi,ltosj
allocate(s%ltosi(s%nlat*s%nlon),s%ltosj(s%nlat*s%nlon))
do i=1,s%nlat*s%nlon
s%ltosi(i)=izero
s%ltosj(i)=izero
end do
! load arrays dealing with global grids
isd_g(1)=izero
s%displs_g(1)=izero
do n=1,s%npe
if(n/=ione) then
isd_g(n)=isd_g(n-ione)+isc_g(n-ione)
s%displs_g(n)=s%displs_g(n-ione)+s%ijn(n-ione)
end if
do j=1,s%jlon1(n)
ns=s%displs_g(n)+(j-ione)*s%ilat1(n)
do i=1,s%ilat1(n)
ns=ns+ione
s%ltosi(ns)=s%istart(n)+i-ione
s%ltosj(ns)=s%jstart(n)+j-ione
end do
end do
end do
! Load arrays dealing with subdomain grids
s%ird_s(1)=izero
s%displs_s(1)=izero
do n=1,s%npe
if(n/=ione) then
s%ird_s(n)=s%ird_s(n-ione)+s%irc_s(n-ione)
s%displs_s(n)=s%displs_s(n-ione)+s%ijn_s(n-ione)
end if
end do
! set total number of points from all subdomain grids
s%itotsub=s%displs_s(s%npe)+s%ijn_s(s%npe)
! obtain ltosi_s,ltosj_s
allocate(s%ltosi_s(s%itotsub),s%ltosj_s(s%itotsub))
do i=1,s%itotsub
s%ltosi_s(i)=izero
s%ltosj_s(i)=izero
end do
if(regional)then
do n=1,s%npe
do j=1,s%jlon1(n)+2_i_kind
ns=s%displs_s(n)+(j-ione)*(s%ilat1(n)+2_i_kind)
do i=1,s%ilat1(n)+2_i_kind
ns=ns+ione
s%ltosi_s(ns)=s%istart(n)+i-2_i_kind
s%ltosj_s(ns)=s%jstart(n)+j-2_i_kind
if(s%ltosi_s(ns)==izero) s%ltosi_s(ns)=ione
if(s%ltosi_s(ns)==nlat+ione) s%ltosi_s(ns)=s%nlat
if(s%ltosj_s(ns)==izero) s%ltosj_s(ns)=ione
if(s%ltosj_s(ns)==nlon+ione) s%ltosj_s(ns)=s%nlon
end do
end do
end do ! end do over npe
else
do n=1,s%npe
do j=1,s%jlon1(n)+2_i_kind
ns=s%displs_s(n)+(j-ione)*(s%ilat1(n)+2_i_kind)
do i=1,s%ilat1(n)+2_i_kind
ns=ns+ione
s%ltosi_s(ns)=s%istart(n)+i-2_i_kind
s%ltosj_s(ns)=s%jstart(n)+j-2_i_kind
if(s%ltosi_s(ns)==izero) s%ltosi_s(ns)=ione
if(s%ltosi_s(ns)==nlat+ione) s%ltosi_s(ns)=nlat
if(s%ltosj_s(ns)==izero) s%ltosj_s(ns)=nlon
if(s%ltosj_s(ns)==nlon+ione) s%ltosj_s(ns)=ione
end do
end do
end do ! end do over npe
endif
deallocate(isc_g,isd_g)
! next, determine vertical layout:
allocate(s%kbegin(0:s%npe),s%kend(0:s%npe-ione))
num_loc_groups=s%num_fields/s%npe
nextra=s%num_fields-num_loc_groups*s%npe
s%kbegin(0)=ione
if(nextra > izero) then
do k=1,nextra
s%kbegin(k)=s%kbegin(k-ione)+ione+num_loc_groups
end do
end if
do k=nextra+ione,s%npe
s%kbegin(k)=s%kbegin(k-ione)+num_loc_groups
end do
do k=0,s%npe-ione
s%kend(k)=s%kbegin(k+ione)-ione
end do
if(s%mype == izero) then
write(6,*)' in general_sub2grid_create_info, kbegin=',s%kbegin
write(6,*)' in general_sub2grid_create_info, kend= ',s%kend
end if
s%kbegin_loc=s%kbegin(s%mype)
s%kend_loc=s%kend(s%mype)
s%kend_alloc=max(s%kend_loc,s%kbegin_loc)
! get alltoallv indices for sub2grid
allocate(s%sendcounts(0:s%npe-ione),s%sdispls(0:s%npe))
allocate(s%recvcounts(0:s%npe-ione),s%rdispls(0:s%npe))
s%sdispls(0)=izero
do k=0,s%npe-ione
s%sendcounts(k)=s%ijn(k+ione)*(s%kend_loc-s%kbegin_loc+ione)
s%sdispls(k+ione)=s%sdispls(k)+s%sendcounts(k)
end do
s%rdispls(0)=izero
do k=0,s%npe-ione
s%recvcounts(k)=s%ijn(s%mype+ione)*(s%kend(k)-s%kbegin(k)+ione)
s%rdispls(k+ione)=s%rdispls(k)+s%recvcounts(k)
end do
! get alltoallv indices for grid2sub
allocate(s%sendcounts_s(0:s%npe-ione),s%sdispls_s(0:s%npe))
allocate(s%recvcounts_s(0:s%npe-ione),s%rdispls_s(0:s%npe))
s%sdispls_s(0)=izero
do k=0,s%npe-ione
s%sendcounts_s(k)=s%ijn_s(k+ione)*(s%kend_loc-s%kbegin_loc+ione)
s%sdispls_s(k+ione)=s%sdispls_s(k)+s%sendcounts_s(k)
end do
s%rdispls_s(0)=izero
do k=0,s%npe-ione
s%recvcounts_s(k)=s%ijn_s(s%mype+ione)*(s%kend(k)-s%kbegin(k)+ione)
s%rdispls_s(k+ione)=s%rdispls_s(k)+s%recvcounts_s(k)
end do
s%lallocated=.true.
end subroutine general_sub2grid_create_info
subroutine general_deter_subdomain_nolayout(npe,mype,nlat,nlon,regional, &
periodic,periodic_s,lon1,lon2,lat1,lat2,ilat1,istart,jlon1,jstart)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_deter_subdomain_nolayout perform domain decomposition
! prgmmr: weiyu yang org: np20 date: 1998-05-14
!
! abstract: Given an array of the observation computation load and
! the number of available mpi tasks (npe), this routine
! decomposes the total analysis grid into npe subdomains
!
! program history log:
! 1998-05-14 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-06-01 treadon - simplify algorithm
! 2004-07-28 treadon - add only to module use, add intent in/out
! 2005-10-17 derber - rewrite routine using simpler algorithm
! 2005-10-26 treadon - correct error in 100 format text
! 2008-06-04 safford - rm unused vars
! 2008-09-05 lueken - merged ed's changes into q1fy09 code
! 2010-02-12 parrish - make copy for use in general_sub2grid_mod
! 2010-03-02 parrish - add regional flag to input. if regional=.true.,
! then periodic, periodic_s=.false. always. this corrects a bug
! in existing code. (never a problem, except when npe=1).
!
! input argument list:
! mype - mpi task number
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: izero,ione
implicit none
! Declare passed variables
integer(i_kind),intent(in ) :: npe,mype,nlat,nlon
logical ,intent(in ) :: regional
logical ,intent( out) :: periodic,periodic_s(npe)
integer(i_kind),intent( out) :: lon1,lon2,lat1,lat2
integer(i_kind),intent( out) :: ilat1(npe),istart(npe),jlon1(npe),jstart(npe)
! Declare local variables
integer(i_kind) npts,nrnc,iinum,iileft,jrows,jleft,k,i,jjnum
integer(i_kind) j,mm1,iicnt,ipts,jjleft
integer(i_kind),dimension(npe+ione):: iiend,jjend,iistart
real(r_kind):: anperpe
!************************************************************************
periodic=.false.
periodic_s=.false.
! Compute number of points on full grid and target number of
! point per mpi task (pe)
npts=nlat*nlon
anperpe=float(npts)/float(npe)
! Start with square subdomains
nrnc=sqrt(anperpe)
iinum=nlon/nrnc
if(iinum==izero) iinum=ione
iicnt=nlon/iinum
iileft=nlon-iicnt*iinum
jrows=npe/iinum
jleft=npe-jrows*iinum
! Adjust subdomain boundaries
k=izero
istart=ione
jstart=ione
iistart(1)=ione
do i=1,iinum
ipts = iicnt
if(i <= iileft)ipts=ipts+ione
iiend(i)=iistart(i)+ipts-ione
iistart(i+ione)=iiend(i)+ione
jjnum=jrows
if(i <= jleft)jjnum=jrows+ione
do j=1,jjnum
k=k+ione
jlon1(k)=ipts
jstart(k)= iistart(i)
ilat1(k)=nlat/jjnum
jjleft=nlat-ilat1(k)*jjnum
if(j <= jjleft)ilat1(k)=ilat1(k)+ione
if(j > ione)istart(k)=jjend(j-1)+ione
jjend(j)=istart(k)+ilat1(k)-ione
if (jlon1(k)==nlon.and..not.regional) then
periodic=.true.
periodic_s(k)=.true.
endif
if(mype == izero) &
write(6,100) k-ione,istart(k),jstart(k),ilat1(k),jlon1(k)
end do
end do
100 format('general_DETER_SUBDOMAIN: task,istart,jstart,ilat1,jlon1=',6(i6,1x))
! Set number of latitude and longitude for given subdomain
mm1=mype+ione
lat1=ilat1(mm1)
lon1=jlon1(mm1)
lat2=lat1+2_i_kind
lon2=lon1+2_i_kind
return
end subroutine general_deter_subdomain_nolayout
subroutine general_sub2grid_r_single(s,sub_vars,grid_vars)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sub2grid_r_single convert from subdomains to full horizontal grid
! prgmmr: parrish org: np22 date: 2010-02-11
!
! abstract: generalized version of sub2grid--uses only gsi module kinds.
! All information needed is contained in the structure variable
! "s", instead of various modules. This allows
! for easy adaptation for any collection/ordering of variables
! defined on subdomains, which need to be made available on
! full horizontal grid for horizontal operations.
! The structure variable is specified by subroutine general_sub2grid_setup.
! This version works with single precision (4-byte) real variables.
! Input sub_vars, the desired arrays on horizontal subdomains, has one
! halo row, for now, which is filled with zero, since for ensemble use,
! there is no need for a halo, but is easiest for now to keep it.
! A later version will have variable number of halo rows, filled with proper values.
!
! program history log:
! 2010-02-11 parrish, initial documentation
!
! input argument list:
! s - structure variable, contains all necessary information for
! moving this set of subdomain variables sub_vars to
! the corresponding set of full horizontal grid variables.
! sub_vars - input grid values in vertical subdomain mode (contains one halo row)
!
! output argument list:
! grid_vars - output grid values in horizontal slab mode.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_single,i_kind,i_long
use constants, only: izero,ione
use mpimod, only: mpi_comm_world,mpi_real4
implicit none
type(sub2grid_info),intent(in ) :: s
real(r_single), intent(in ) :: sub_vars(s%inner_vars,s%lat2,s%lon2,s%num_fields)
real(r_single), intent( out) :: grid_vars(s%inner_vars,s%nlat,s%nlon,s%kbegin_loc:s%kend_alloc)
real(r_single) :: sub_vars0(s%inner_vars,s%lat1,s%lon1,s%num_fields)
real(r_single) :: work(s%inner_vars,s%itotsub*(s%kend_alloc-s%kbegin_loc+ione))
integer(i_kind) iloc,iskip,i,i0,ii,j,j0,k,n,k_in,ilat,jlon,ierror
integer(i_long) mpi_string
! remove halo row
do k=1,s%num_fields
do j=2,s%lon2-1
j0=j-1
do i=2,s%lat2-1
i0=i-1
do ii=1,s%inner_vars
sub_vars0(ii,i0,j0,k)=sub_vars(ii,i,j,k)
end do
end do
end do
end do
call mpi_type_contiguous(s%inner_vars,mpi_real4,mpi_string,ierror)
call mpi_type_commit(mpi_string,ierror)
call mpi_alltoallv(sub_vars0,s%recvcounts,s%rdispls,mpi_string, &
work,s%sendcounts,s%sdispls,mpi_string,mpi_comm_world,ierror)
call mpi_type_free(mpi_string,ierror)
k_in=s%kend_loc-s%kbegin_loc+ione
! Load temp array in desired order
do k=s%kbegin_loc,s%kend_loc
iskip=izero
iloc=izero
do n=1,s%npe
if (n/=ione) then
iskip=iskip+s%ijn(n-ione)*k_in
end if
do i=1,s%ijn(n)
iloc=iloc+ione
ilat=s%ltosi(iloc)
jlon=s%ltosj(iloc)
do ii=1,s%inner_vars
grid_vars(ii,ilat,jlon,k)=work(ii,i + iskip + (k-s%kbegin_loc)*s%ijn(n))
end do
end do
end do
end do
end subroutine general_sub2grid_r_single
subroutine general_grid2sub_r_single(s,grid_vars,sub_vars)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sub2grid convert from subdomains to full horizontal grid
! prgmmr: parrish org: np22 date: 2010-02-11
!
! abstract: generalized version of grid2sub--uses only gsi module kinds.
! All information needed is contained in the structure variable
! "s", instead of various modules. This allows
! for easy adaptation for any collection/ordering of variables
! defined on subdomains, which need to be made available on
! full horizontal grid for horizontal operations.
! The structure variable is specified by subroutine general_sub2grid_setup.
! This version works with single precision (4-byte) real variables.
! Output sub_vars, the desired arrays on horizontal subdomains, has one
! halo row, for now, which is filled with zero, since for ensemble use,
! there is no need for a halo, but is easiest for now to keep it.
! A later version will have variable number of halo rows, filled with proper values.
!
! program history log:
! 2010-02-11 parrish, initial documentation
! 2010-03-02 parrish - remove setting halo to zero in output
!
! input argument list:
! s - structure variable, contains all necessary information for
! moving this set of subdomain variables sub_vars to
! the corresponding set of full horizontal grid variables.
! grid_vars - input grid values in horizontal slab mode.
!
! output argument list:
! sub_vars - output grid values in vertical subdomain mode
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_single,i_kind,i_long
use constants, only: izero,ione,zero
use mpimod, only: mpi_comm_world,mpi_real4
implicit none
type(sub2grid_info),intent(in ) :: s
real(r_single), intent(in ) :: grid_vars(s%inner_vars,s%nlat,s%nlon,s%kbegin_loc:s%kend_alloc)
real(r_single), intent( out) :: sub_vars(s%inner_vars,s%lat2,s%lon2,s%num_fields)
real(r_single),allocatable :: temp(:,:),work(:,:,:)
integer(i_kind) iloc,iskip,i,ii,j,k,n,k_in,ilat,jlon,ierror
integer(i_long) mpi_string
allocate(temp(s%inner_vars,s%itotsub*(s%kend_alloc-s%kbegin_loc+ione)))
allocate(work(s%inner_vars,s%itotsub,s%kbegin_loc:s%kend_alloc))
! reorganize for eventual distribution to local domains
work=zero !????????????not needed??
do k=s%kbegin_loc,s%kend_loc
do i=1,s%itotsub
ilat=s%ltosi_s(i)
jlon=s%ltosj_s(i)
do ii=1,s%inner_vars
work(ii,i,k)=grid_vars(ii,ilat,jlon,k)
end do
end do
end do
! load temp array in order of subdomains
temp=zero
iloc=izero
iskip=izero
do n=1,s%npe
if (n/=ione) then
iskip=iskip+s%ijn_s(n-ione)
end if
do k=s%kbegin_loc,s%kend_loc
do i=1,s%ijn_s(n)
iloc=iloc+ione
do ii=1,s%inner_vars
temp(ii,iloc)=work(ii,iskip+i,k)
end do
end do
end do
end do
! Now load the temp array back into work
iloc=izero
do k=s%kbegin_loc,s%kend_loc
do i=1,s%itotsub
iloc=iloc+ione
do ii=1,s%inner_vars
work(ii,i,k)=temp(ii,iloc)
end do
end do
end do
deallocate(temp)
call mpi_type_contiguous(s%inner_vars,mpi_real4,mpi_string,ierror)
call mpi_type_commit(mpi_string,ierror)
call mpi_alltoallv(work,s%sendcounts_s,s%sdispls_s,mpi_string, &
sub_vars,s%recvcounts_s,s%rdispls_s,mpi_string,mpi_comm_world,ierror)
deallocate(work)
call mpi_type_free(mpi_string,ierror)
end subroutine general_grid2sub_r_single
subroutine general_sub2grid_r_double(s,sub_vars,grid_vars)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sub2grid_r_double convert from subdomains to full horizontal grid
! prgmmr: parrish org: np22 date: 2010-02-11
!
! abstract: generalized version of sub2grid--uses only gsi module kinds.
! All information needed is contained in the structure variable
! "s", instead of various modules. This allows
! for easy adaptation for any collection/ordering of variables
! defined on subdomains, which need to be made available on
! full horizontal grid for horizontal operations.
! The structure variable is specified by subroutine general_sub2grid_setup.
! This version works with double precision (8-byte) real variables.
! Input sub_vars, the desired arrays on horizontal subdomains, has one
! halo row, for now, which is filled with zero, since for ensemble use,
! there is no need for a halo, but is easiest for now to keep it.
! A later version will have variable number of halo rows, filled with proper values.
!
! program history log:
! 2010-02-11 parrish, initial documentation
!
! input argument list:
! s - structure variable, contains all necessary information for
! moving this set of subdomain variables sub_vars to
! the corresponding set of full horizontal grid variables.
! sub_vars - input grid values in vertical subdomain mode
!
! output argument list:
! grid_vars - output grid values in horizontal slab mode.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_double,i_kind,i_long
use constants, only: izero,ione
use mpimod, only: mpi_comm_world,mpi_real8
implicit none
type(sub2grid_info),intent(in ) :: s
real(r_double), intent(in ) :: sub_vars(s%inner_vars,s%lat2,s%lon2,s%num_fields)
real(r_double), intent( out) :: grid_vars(s%inner_vars,s%nlat,s%nlon,s%kbegin_loc:s%kend_alloc)
real(r_double) :: sub_vars0(s%inner_vars,s%lat1,s%lon1,s%num_fields)
real(r_double) :: work(s%inner_vars,s%itotsub*(s%kend_alloc-s%kbegin_loc+ione))
integer(i_kind) iloc,iskip,i,i0,ii,j,j0,k,n,k_in,ilat,jlon,ierror
integer(i_long) mpi_string
! remove halo row
do k=1,s%num_fields
do j=2,s%lon2-1
j0=j-1
do i=2,s%lat2-1
i0=i-1
do ii=1,s%inner_vars
sub_vars0(ii,i0,j0,k)=sub_vars(ii,i,j,k)
end do
end do
end do
end do
call mpi_type_contiguous(s%inner_vars,mpi_real8,mpi_string,ierror)
call mpi_type_commit(mpi_string,ierror)
call mpi_alltoallv(sub_vars0,s%recvcounts,s%rdispls,mpi_string, &
work,s%sendcounts,s%sdispls,mpi_string,mpi_comm_world,ierror)
call mpi_type_free(mpi_string,ierror)
k_in=s%kend_loc-s%kbegin_loc+ione
! Load temp array in desired order
do k=s%kbegin_loc,s%kend_loc
iskip=izero
iloc=izero
do n=1,s%npe
if (n/=ione) then
iskip=iskip+s%ijn(n-ione)*k_in
end if
do i=1,s%ijn(n)
iloc=iloc+ione
ilat=s%ltosi(iloc)
jlon=s%ltosj(iloc)
do ii=1,s%inner_vars
grid_vars(ii,ilat,jlon,k)=work(ii,i + iskip + (k-s%kbegin_loc)*s%ijn(n))
end do
end do
end do
end do
end subroutine general_sub2grid_r_double
subroutine general_grid2sub_r_double(s,grid_vars,sub_vars)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sub2grid convert from subdomains to full horizontal grid
! prgmmr: parrish org: np22 date: 2010-02-11
!
! abstract: generalized version of grid2sub--uses only gsi module kinds.
! All information needed is contained in the structure variable
! "s", instead of various modules. This allows
! for easy adaptation for any collection/ordering of variables
! defined on subdomains, which need to be made available on
! full horizontal grid for horizontal operations.
! The structure variable is specified by subroutine general_sub2grid_setup.
! This version works with double precision (8-byte) real variables.
! Output sub_vars, the desired arrays on horizontal subdomains, has one
! halo row, for now, which is filled with zero, since for ensemble use,
! there is no need for a halo, but is easiest for now to keep it.
! A later version will have variable number of halo rows, filled with proper values.
!
! program history log:
! 2010-02-11 parrish, initial documentation
! 2010-03-02 parrish - remove setting halo to zero in output
!
! input argument list:
! s - structure variable, contains all necessary information for
! moving this set of subdomain variables sub_vars to
! the corresponding set of full horizontal grid variables.
! grid_vars - input grid values in horizontal slab mode.
!
! output argument list:
! sub_vars - output grid values in vertical subdomain mode.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_double,i_kind,i_long
use constants, only: izero,ione,zero
use mpimod, only: mpi_comm_world,mpi_real8
implicit none
type(sub2grid_info),intent(in ) :: s
real(r_double), intent(in ) :: grid_vars(s%inner_vars,s%nlat,s%nlon,s%kbegin_loc:s%kend_alloc)
real(r_double), intent( out) :: sub_vars(s%inner_vars,s%lat2,s%lon2,s%num_fields)
real(r_double),allocatable :: temp(:,:),work(:,:,:)
integer(i_kind) iloc,iskip,i,ii,j,k,n,k_in,ilat,jlon,ierror
integer(i_long) mpi_string
allocate(temp(s%inner_vars,s%itotsub*(s%kend_alloc-s%kbegin_loc+ione)))
allocate(work(s%inner_vars,s%itotsub,s%kbegin_loc:s%kend_alloc))
! reorganize for eventual distribution to local domains
work=zero
do k=s%kbegin_loc,s%kend_loc
do i=1,s%itotsub
ilat=s%ltosi_s(i)
jlon=s%ltosj_s(i)
do ii=1,s%inner_vars
work(ii,i,k)=grid_vars(ii,ilat,jlon,k)
end do
end do
end do
! load temp array in order of subdomains
temp=zero
iloc=izero
iskip=izero
do n=1,s%npe
if (n/=ione) then
iskip=iskip+s%ijn_s(n-ione)
end if
do k=s%kbegin_loc,s%kend_loc
do i=1,s%ijn_s(n)
iloc=iloc+ione
do ii=1,s%inner_vars
temp(ii,iloc)=work(ii,iskip+i,k)
end do
end do
end do
end do
! Now load the temp array back into work
iloc=izero
do k=s%kbegin_loc,s%kend_loc
do i=1,s%itotsub
iloc=iloc+ione
do ii=1,s%inner_vars
work(ii,i,k)=temp(ii,iloc)
end do
end do
end do
deallocate(temp)
call mpi_type_contiguous(s%inner_vars,mpi_real8,mpi_string,ierror)
call mpi_type_commit(mpi_string,ierror)
call mpi_alltoallv(work,s%sendcounts_s,s%sdispls_s,mpi_string, &
sub_vars,s%recvcounts_s,s%rdispls_s,mpi_string,mpi_comm_world,ierror)
deallocate(work)
call mpi_type_free(mpi_string,ierror)
end subroutine general_grid2sub_r_double
subroutine general_sube2suba_r_double(se,sa,p_e2a,sube_vars,suba_vars,regional)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sube2suba_r_double interpolate ens grid to anl grid
! prgmmr: parrish org: np22 date: 2010-02-27
!
! abstract: interpolate ensemble grid variables to analysis grid variables,
! where input and output are in the respective subdomains as defined
! by the structure variables se and sa.
!
! program history log:
! 2010-02-27 parrish, initial documentation
!
! input argument list:
! se - ensemble grid structure variable
! sa - analysis grid structure variable
! p_e2a - interpolation from ensemble to grid to analysis grid structure variable
! sube_vars - input ensemble grid values in ensemble subdomain mode (as defined by se)
! regional - true for regional grids--this code currently works only with global grids
! need to fix this.
!
! output argument list:
! suba_vars - output analysis grid values in analysis subdomain mode (as defined by sa)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_double,i_kind,i_long
use constants, only: izero,ione
use egrid2agrid_mod, only: g_egrid2agrid,egrid2agrid_parm
use mpimod, only: mpi_comm_world,mpi_real8
implicit none
type(sub2grid_info), intent(in ) :: se,sa
type(egrid2agrid_parm),intent(in ) :: p_e2a
real(r_double), intent(in ) :: sube_vars(se%inner_vars,se%lat2,se%lon2,se%num_fields)
real(r_double), intent( out) :: suba_vars(sa%inner_vars,sa%lat2,sa%lon2,sa%num_fields)
logical, intent(in ) :: regional
real(r_double),allocatable:: gride_vars(:,:),grida_vars(:,:)
integer(i_kind) k
allocate(gride_vars(se%inner_vars*se%nlat*se%nlon,se%kbegin_loc:se%kend_alloc))
call general_sub2grid_r_double(se,sube_vars,gride_vars)
allocate(grida_vars(sa%inner_vars*sa%nlat*sa%nlon,sa%kbegin_loc:sa%kend_alloc))
if(regional) then
write(6,*)' not ready for regional dual_res yet'
call mpi_finalize(k)
stop
else
do k=se%kbegin_loc,se%kend_loc
call g_egrid2agrid(p_e2a,gride_vars(:,k),grida_vars(:,k),se%vector(k))
end do
end if
deallocate(gride_vars)
call general_grid2sub_r_double(sa,grida_vars,suba_vars)
deallocate(grida_vars)
end subroutine general_sube2suba_r_double
subroutine general_sube2suba_r_double_ad(se,sa,p_e2a,sube_vars,suba_vars,regional)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_sube2suba_r_double_ad adjoint of interpolate ens grid to anl grid
! prgmmr: parrish org: np22 date: 2010-02-28
!
! abstract: adjoint of general_sube2suba_r_double.
!
! program history log:
! 2010-02-28 parrish, initial documentation
!
! input argument list:
! se - ensemble grid structure variable
! sa - analysis grid structure variable
! p_e2a - interpolation from ensemble to grid to analysis grid structure variable
! sube_vars - input ensemble grid values in ensemble subdomain mode (as defined by se)
! regional - true for regional grids--this code currently works only with global grids
! need to fix this.
!
! output argument list:
! suba_vars - output analysis grid values in analysis subdomain mode (as defined by sa)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_double,i_kind,i_long
use constants, only: izero,ione
use egrid2agrid_mod, only: g_egrid2agrid_ad,egrid2agrid_parm
use mpimod, only: mpi_comm_world,mpi_real8
implicit none
type(sub2grid_info), intent(in ) :: se,sa
type(egrid2agrid_parm),intent(in ) :: p_e2a
real(r_double), intent( out) :: sube_vars(se%inner_vars,se%lat2,se%lon2,se%num_fields)
real(r_double), intent(in ) :: suba_vars(sa%inner_vars,sa%lat2,sa%lon2,sa%num_fields)
logical, intent(in ) :: regional
real(r_double),allocatable:: gride_vars(:,:),grida_vars(:,:)
integer(i_kind) k
allocate(grida_vars(sa%inner_vars*sa%nlat*sa%nlon,sa%kbegin_loc:sa%kend_alloc))
call general_sub2grid_r_double(sa,suba_vars,grida_vars)
allocate(gride_vars(se%inner_vars*se%nlat*se%nlon,se%kbegin_loc:se%kend_alloc))
if(regional) then
write(6,*)' not ready for regional dual_res yet'
call mpi_finalize(k)
stop
else
do k=se%kbegin_loc,se%kend_loc
call g_egrid2agrid_ad(p_e2a,gride_vars(:,k),grida_vars(:,k),se%vector(k))
end do
end if
deallocate(grida_vars)
call general_grid2sub_r_double(se,gride_vars,sube_vars)
deallocate(gride_vars)
end subroutine general_sube2suba_r_double_ad
subroutine general_suba2sube_r_double(sa,se,p_e2a,suba_vars,sube_vars,regional)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_suba2sube_r_double smoothing interpolate anl grid to ens grid
! prgmmr: parrish org: np22 date: 2010-03-01
!
! abstract: smoothing interpolation from analysis grid to ensemble grid (analysis subdomain
! input, ensemble subdomain output).
!
! program history log:
! 2010-03-01 parrish, initial documentation
!
! input argument list:
! sa - ensemble grid structure variable
! se - analysis grid structure variable
! p_e2a - interpolation from ensemble to grid to analysis grid structure variable
! suba_vars - input analysis grid values in analysis subdomain mode (as defined by sa)
! regional - true for regional grids--this code currently works only with global grids
! need to fix this.
!
! output argument list:
! sube_vars - output ensemble grid values in ensemble subdomain mode (as defined by se)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_double,i_kind,i_long
use constants, only: izero,ione
use egrid2agrid_mod, only: g_agrid2egrid,egrid2agrid_parm
use mpimod, only: mpi_comm_world,mpi_real8
implicit none
type(sub2grid_info), intent(in ) :: sa,se
type(egrid2agrid_parm),intent(in ) :: p_e2a
real(r_double), intent(in ) :: suba_vars(sa%inner_vars,sa%lat2,sa%lon2,sa%num_fields)
real(r_double), intent( out) :: sube_vars(se%inner_vars,se%lat2,se%lon2,se%num_fields)
logical, intent(in ) :: regional
real(r_double),allocatable:: gride_vars(:,:),grida_vars(:,:)
integer(i_kind) k
allocate(grida_vars(sa%inner_vars*sa%nlat*sa%nlon,sa%kbegin_loc:sa%kend_alloc))
call general_sub2grid_r_double(sa,suba_vars,grida_vars)
allocate(gride_vars(se%inner_vars*se%nlat*se%nlon,se%kbegin_loc:se%kend_alloc))
if(regional) then
write(6,*)' not ready for regional dual_res yet'
call mpi_finalize(k)
stop
else
do k=se%kbegin_loc,se%kend_loc
call g_agrid2egrid(p_e2a,grida_vars(:,k),gride_vars(:,k),se%vector(k))
end do
end if
deallocate(grida_vars)
call general_grid2sub_r_double(se,gride_vars,sube_vars)
deallocate(gride_vars)
end subroutine general_suba2sube_r_double
end module general_sub2grid_mod