program gtopo30,43
C **********************************************************************
C **** Program to create topography data to be used for ****
C **** Turbulence. TURB topography data contains scales ****
C **** (0-10 km), and is the residual from: ****
C **** ****
C **** Raw (all scales) - Mean (scales > 100 km) ****
C **** - Grav (100 km > scales > 10 km) ****
C **** ****
C **********************************************************************
use dynamics_lattice_module
type ( dynamics_lattice_type ) lattice
integer imax, jmax
parameter ( imax = 120*360 )
parameter ( jmax = 120*180 )
real dum (imax,jmax)
real, allocatable :: hraw (:,:)
real, allocatable :: hmean(:,:)
c MPI stuff
c ---------
include 'mpif.h'
integer myid
integer npex,npey,ierror
character*4 cnpex
character*4 cnpey
call mpi_init (ierror)
call mpi_comm_rank ( mpi_comm_world,myid,ierror )
call getenv ('NPEX',cnpex)
call getenv ('NPEY',cnpey)
read(cnpex,101) npex
read(cnpey,101) npey
101 format(i4)
call create_dynamics_lattice ( lattice,npex,npey )
call init_dynamics_lattice ( lattice,mpi_comm_world,imax,jmax,1 )
myid = lattice%myid
im = lattice%im( lattice%pei )
jm = lattice%jm( lattice%pej )
allocate ( hraw (im,jm) )
allocate ( hmean(im,jm) )
C **********************************************************************
C **** Read GTOPO30 Raw & Mean Data ****
C **** Compute Residual and Distribute ****
C **********************************************************************
if( myid.eq.0 ) then
print *, 'Begin Reading Raw Data'
open (40,file='hraw.30x30sec.data2',form='unformatted',access='sequential')
do j=1,jmax
read (40) (dum(i,j),i=1,imax)
enddo
close(40)
print *, 'Finished Reading Raw Data'
endif
call scatter ( dum,hraw,lattice )
call printchar ('finished scatter of raw data',lattice)
if( myid.eq.0 ) then
print *, 'Begin Reading Mean Data'
open (40,file='hmean.30x30sec.data2',form='unformatted',access='sequential')
do j=1,jmax
read (40) (dum(i,j),i=1,imax)
enddo
close(40)
print *, 'Finished Reading Mean Data'
endif
call scatter ( dum,hmean,lattice )
call printchar ('finished scatter of mean data',lattice)
c Compute Residual
c ----------------
hraw(:,:) = hraw(:,:)-hmean(:,:)
if( myid.eq.0 ) then
print *, 'Begin Reading GRAV Data'
open (40,file='hgrav.30x30sec.data2',form='unformatted',access='sequential')
do j=1,jmax
read (40) (dum(i,j),i=1,imax)
enddo
close(40)
print *, 'Finished Reading GRAV Data'
endif
call scatter ( dum,hmean,lattice )
call printchar ('finished scatter of mean data',lattice)
c Compute Residual
c ----------------
hraw(:,:) = hraw(:,:)-hmean(:,:)
C **********************************************************************
C **** Gather Data and Write Output ****
C **********************************************************************
call gather ( dum,hraw,lattice )
call printchar ('finished gather of TURB data',lattice)
if( myid.eq.0 ) then
print *, 'Begin Writing TURB Data'
open (40,file='hturb.30x30sec.data2',form='unformatted',access='sequential')
do j=1,jmax
write(40) (dum(i,j),i=1,imax)
enddo
close(40)
print *, 'Finished Writing TURB Data'
endif
call my_finalize
stop
end
subroutine scatter ( qglobal,qlocal,lattice ) 12
use dynamics_lattice_module
implicit none
type ( dynamics_lattice_type ) lattice
real qglobal( lattice%imglobal,lattice%jmglobal )
real qlocal ( lattice%im(lattice%pei),lattice%jm(lattice%pej) )
include 'mpif.h'
integer status(mpi_status_size)
integer comm
integer myid,npe,ierror
integer nx,i,iloc,im,imx,imglobal,isum
integer ny,j,jloc,jm,jmy,jmglobal,jsum
real, allocatable :: buf(:,:)
comm = lattice%comm
myid = lattice%myid
iloc = lattice%pei
jloc = lattice%pej
im = lattice%im(iloc)
jm = lattice%jm(jloc)
imglobal = lattice%imglobal
jmglobal = lattice%jmglobal
if( myid.eq.0 ) then
jsum = 0
do ny=0,lattice%ny-1
jmy = lattice%jm(ny)
isum = 0
do nx=0,lattice%nx-1
imx = lattice%im(nx)
if( nx.eq.0 .and. ny.eq.0 ) then
do j=1,jmy
do i=1,imx
qlocal(i,j) = qglobal(i,j)
enddo
enddo
isum = isum + imx
else
allocate ( buf(imx,jmy) )
do j=1,jmy
do i=1,imx
buf(i,j) = qglobal(i+isum,j+jsum)
enddo
enddo
isum = isum + imx
npe = nx + ny*lattice%nx
call mpi_send ( buf,imx*jmy,mpi_real,npe,npe,comm,ierror )
deallocate ( buf )
endif
enddo
jsum = jsum + jmy
enddo
else
call mpi_recv ( qlocal,im*jm,mpi_real,0,myid,comm,status,ierror )
endif
call mpi_barrier ( comm,ierror )
return
end
subroutine gather ( qglobal,qlocal,lattice ) 9
use dynamics_lattice_module
implicit none
type ( dynamics_lattice_type ) lattice
real qglobal( lattice%imglobal,lattice%jmglobal )
real qlocal ( lattice%im(lattice%pei),lattice%jm(lattice%pej) )
include 'mpif.h'
integer status(mpi_status_size)
integer comm
integer myid,npe,ierror
integer nx,i,iloc,im,imx,imglobal,isum
integer ny,j,jloc,jm,jmy,jmglobal,jsum
real, allocatable :: buf(:,:)
comm = lattice%comm
myid = lattice%myid
iloc = lattice%pei
jloc = lattice%pej
im = lattice%im(iloc)
jm = lattice%jm(jloc)
imglobal = lattice%imglobal
jmglobal = lattice%jmglobal
if( myid.eq.0 ) then
jsum = 0
do ny=0,lattice%ny-1
jmy = lattice%jm(ny)
isum = 0
do nx=0,lattice%nx-1
imx = lattice%im(nx)
if( nx.eq.0 .and. ny.eq.0 ) then
do j=1,jmy
do i=1,imx
qglobal(i,j) = qlocal(i,j)
enddo
enddo
isum = isum + imx
else
allocate ( buf(imx,jmy) )
npe = nx + ny*lattice%nx
call mpi_recv ( buf,imx*jmy,mpi_real,npe,npe,comm,status,ierror )
do j=1,jmy
do i=1,imx
qglobal(i+isum,j+jsum) = buf(i,j)
enddo
enddo
isum = isum + imx
deallocate ( buf )
endif
enddo
jsum = jsum + jmy
enddo
else
call mpi_send ( qlocal,im*jm,mpi_real,0,myid,comm,ierror )
endif
call mpi_barrier ( comm,ierror )
return
end
subroutine ghostx (a,b,im,jm,lm,n,lattice,flag) 5,6
C **********************************************************************
C **** ****
C **** This program fills GHOST values in the Y-direction ****
C **** ****
C **** a ....... Input Array Non-Ghosted ****
C **** b ....... Output Array Ghosted ****
C **** im ...... Local Zonal Dimension ****
C **** jm ...... Local Meridional Dimension ****
C **** lm ...... Local Vertical Dimension ****
C **** n ....... Number of GHOST values ****
C **** lattice . Grid Lattice for MPI ****
C **** flag .... Character*(*) to do 'east' or 'west' only ****
C **** ****
C **********************************************************************
use dynamics_lattice_module
implicit none
type ( dynamics_lattice_type ) lattice
include 'mpif.h'
integer stat(mpi_status_size,4)
integer myid, nx, comm, ierror
integer send_reqeast, send_reqwest
integer recv_reqeast, recv_reqwest
character*(*) flag
integer im,jm,lm,n,i,j,L
integer east,west
real getcon,undef
real a(1 :im ,1:jm,lm)
real b(1-n:im+n,1:jm,lm)
real bufs(n,jm,lm,2)
real bufr(n,jm,lm,2)
undef = getcon('UNDEF')
comm = lattice%comm
myid = lattice%myid
nx = lattice%nx
if( n.gt.im ) then
print *
print *, 'Processor ',myid,' does not contain enough grid-points in X (',im,') to perform ',n,' point ghosting!'
call my_finalize
endif
east = mod( myid+nx+1,nx ) + (myid/nx)*nx
west = mod( myid+nx-1,nx ) + (myid/nx)*nx
do L=1,lm
do j=1,jm
do i=1,im
b(i,j,L) = a(i,j,L)
enddo
do i=1,n
bufs(i,j,L,1) = a(i,j,L)
bufs(i,j,L,2) = a(im-n+i,j,L)
bufr(i,j,L,1) = undef
bufr(i,j,L,2) = undef
b( 1-i,j,L) = undef ! Initialize ghost regions
b(im+i,j,L) = undef ! Initialize ghost regions
enddo
enddo
enddo
if( nx.gt.1 ) then
if( east.eq.west ) then
call mpi_sendrecv( bufs,2*n*jm*lm,mpi_real,east,east,
. bufr,2*n*jm*lm,mpi_real,west,myid,comm,stat,ierror )
else
if( flag.ne.'east' ) then
call mpi_isend ( bufs(1,1,1,2),n*jm*lm,mpi_real,east,east,comm,send_reqeast,ierror )
call mpi_irecv ( bufr(1,1,1,2),n*jm*lm,mpi_real,west,myid,comm,recv_reqwest,ierror )
endif
if( flag.ne.'west' ) then
call mpi_isend ( bufs(1,1,1,1),n*jm*lm,mpi_real,west,west,comm,send_reqwest,ierror )
call mpi_irecv ( bufr(1,1,1,1),n*jm*lm,mpi_real,east,myid,comm,recv_reqeast,ierror )
endif
if( flag.ne.'east' ) then
call mpi_wait ( send_reqeast,stat(1,1),ierror )
call mpi_wait ( recv_reqwest,stat(1,2),ierror )
endif
if( flag.ne.'west' ) then
call mpi_wait ( send_reqwest,stat(1,3),ierror )
call mpi_wait ( recv_reqeast,stat(1,4),ierror )
endif
endif
else
do L=1,lm
do j=1,jm
do i=1,n
bufr(i,j,L,1) = bufs(i,j,L,1)
bufr(i,j,L,2) = bufs(i,j,L,2)
enddo
enddo
enddo
endif
do L=1,lm
do j=1,jm
do i=-n+1,0
b(i,j,L) = bufr(i+n,j,L,2)
enddo
do i=im+1,im+n
b(i,j,L) = bufr(i-im,j,L,1)
enddo
enddo
enddo
return
end
subroutine ghosty (a,b,im,jm,lm,shift,msgn,n,lattice,flag) 5,6
C **********************************************************************
C **** ****
C **** This program fills GHOST values in the Y-direction ****
C **** ****
C **** a ....... Input Array Non-Ghosted ****
C **** b ....... Output Array Ghosted ****
C **** im ...... Local Zonal Dimension ****
C **** jm ...... Local Meridional Dimension ****
C **** lm ...... Local Vertical Dimension ****
C **** shift ... Integer Flag: 0 for A-Grid, 1 for C-Grid VWND ****
C **** msgn .... Integer Flag for Scaler (1) or Vector (-1) ****
C **** n ....... Number of GHOST values ****
C **** lattice . Grid Lattice for MPI ****
C **** flag .... Character*(*) to do 'north', 'south', or 'pole' ****
C **** ****
C **********************************************************************
use dynamics_lattice_module
implicit none
type ( dynamics_lattice_type ) lattice
include 'mpif.h'
integer status(mpi_status_size)
integer myid, nx, comm, ierror
character*(*) flag
integer im,jm,lm,shift,m,n,i,j,L,i0,req(im)
integer north,south,imx,ibeg,npes,request,msgn
real getcon,undef
real a(1:im, 1:jm ,lm)
real b(1:im,1-n:jm+n,lm)
real, allocatable :: apls(:,:,:)
real, allocatable :: amns(:,:,:)
real, allocatable :: buf(:,:)
undef = getcon('UNDEF')
comm = lattice%comm
myid = lattice%myid
nx = lattice%nx
if( (lattice%pej.eq.0 .and. n.gt.jm-1) .or. ! Pole values cannot be used for ghosting
. (lattice%pej.eq.lattice%ny-1 .and. n.gt.jm-1) .or. ! Pole values cannot be used for ghosting
. (n.gt.jm) ) then
print *
print *, 'Processor ',myid,' does not contain enough grid-points in Y (',jm,') to perform ',n,' point ghosting!'
call my_finalize
endif
do L=1,lm
do j=1,jm
do i=1,im
b(i,j,L) = a(i,j,L)
enddo
enddo
b(:, 1-n:0 ,L) = undef ! Initialize ghost regions
b(:,jm+1:jm+n,L) = undef ! Initialize ghost regions
enddo
c Ghost South Pole
c ----------------
if( lattice%pej.eq.0 .and. flag.ne.'north' ) then
do L=1,lm
do j=1,n
do i=1,im
b(i,1-j,L) = a(i,1+j-shift,L)
enddo
enddo
enddo
endif
c Ghost North Pole
c ----------------
if( lattice%pej.eq.lattice%ny-1 .and. flag.ne.'south' ) then
do L=1,lm
do j=1,n
do i=1,im
b(i,jm+j-shift,L) = a(i,jm-j,L)
enddo
enddo
enddo
endif
c Ghost Non-Pole Points North
c ---------------------------
if( flag.ne.'south' .and. flag.ne.'pole' ) then
if( lattice%pej.eq.lattice%ny-1 .and. lattice%pej.ne.0 ) then
south = myid - nx
allocate ( buf(im*n*lm,1) )
do L=1,lm
do j=1,n
do i=1,im
buf(i+(j-1)*im+(L-1)*im*n,1) = a(i,j,L)
enddo
enddo
enddo
call mpi_isend( buf,n*im*lm,mpi_real,south,south,comm,request,ierror )
call mpi_wait ( request,status,ierror )
deallocate ( buf )
endif
if( lattice%pej.eq.0 .and. lattice%pej.ne.lattice%ny-1 ) then
north = myid + nx
allocate ( apls(im,n,lm) )
call mpi_recv ( apls,n*im*lm,mpi_real,north,myid,comm,status,ierror )
do L=1,lm
do j=1,n
do i=1,im
b(i,jm+j,L) = apls(i,j,L)
enddo
enddo
enddo
deallocate ( apls )
endif
if( lattice%pej.ne.0 .and. lattice%pej.ne.lattice%ny-1 ) then
north = myid + nx
south = myid - nx
allocate ( apls(im,n,lm) )
allocate ( buf(im*n*lm,1) )
do L=1,lm
do j=1,n
do i=1,im
buf(i+(j-1)*im+(L-1)*im*n,1) = a(i,j,L)
enddo
enddo
enddo
call mpi_isend( buf,n*im*lm,mpi_real,south,south,comm,request,ierror )
call mpi_recv ( apls,n*im*lm,mpi_real,north,myid,comm,status,ierror )
call mpi_wait ( request,status,ierror )
do L=1,lm
do j=1,n
do i=1,im
b(i,jm+j,L) = apls(i,j,L)
enddo
enddo
enddo
deallocate ( apls )
deallocate ( buf )
endif
endif
c Ghost Non-Pole Points South
c ---------------------------
if( flag.ne.'north' .and. flag.ne.'pole' ) then
if( lattice%pej.eq.0 .and. lattice%pej.ne.lattice%ny-1 ) then
north = myid + nx
allocate ( buf(im*n*lm,1) )
do L=1,lm
do j=1,n
do i=1,im
buf(i+(j-1)*im+(L-1)*im*n,1) = a(i,jm-n+j,L)
enddo
enddo
enddo
call mpi_isend( buf,n*im*lm,mpi_real,north,north,comm,request,ierror )
call mpi_wait ( request,status,ierror )
deallocate ( buf )
endif
if( lattice%pej.eq.lattice%ny-1 .and. lattice%pej.ne.0 ) then
south = myid - nx
allocate ( amns(im,n,lm) )
call mpi_recv ( amns,n*im*lm,mpi_real,south,myid,comm,status,ierror )
do L=1,lm
do j=1,n
do i=1,im
b(i,j-n,L) = amns(i,j,L)
enddo
enddo
enddo
deallocate ( amns )
endif
if( lattice%pej.ne.0 .and. lattice%pej.ne.lattice%ny-1 ) then
north = myid + nx
south = myid - nx
allocate ( amns(im,n,lm) )
allocate ( buf(im*n*lm,1) )
do L=1,lm
do j=1,n
do i=1,im
buf(i+(j-1)*im+(L-1)*im*n,1) = a(i,jm-n+j,L)
enddo
enddo
enddo
call mpi_isend( buf,n*im*lm,mpi_real,north,north,comm,request,ierror )
call mpi_recv ( amns,n*im*lm,mpi_real,south,myid,comm,status,ierror )
call mpi_wait ( request,status,ierror )
do L=1,lm
do j=1,n
do i=1,im
b(i,j-n,L) = amns(i,j,L)
enddo
enddo
enddo
deallocate ( amns )
deallocate ( buf )
endif
endif
return
end
subroutine init_dynamics_lattice ( lattice,comm,imglobal,jmglobal,lm ) 6,12
use dynamics_lattice_module
implicit none
include 'mymalloc_interface'
type ( dynamics_lattice_type ) lattice
integer comm,imglobal,jmglobal,lm
include 'mpif.h'
integer status(mpi_status_size)
integer myid,ierror,im,jm,i,j,m,n,npes,rm
integer isum,jsum,imx,nx, img, ppe(imglobal)
call mpi_comm_rank ( comm,myid,ierror )
call mymalloc ( lattice%ppeg,imglobal )
lattice%comm = comm
lattice%myid = myid
lattice%nx = size( lattice%im )
lattice%ny = size( lattice%jm )
npes = lattice%nx*lattice%ny
lattice%imglobal = imglobal
lattice%jmglobal = jmglobal
lattice%lm = lm
c Initialize lattice%im
c ---------------------
im = imglobal/lattice%nx
rm = imglobal-lattice%nx*im
lattice%imax = im
do n=0,lattice%nx-1
lattice%im(n) = im
if( n.le.rm-1 ) lattice%im(n) = im+1
lattice%imax = max( lattice%imax,lattice%im(n) )
enddo
c Initialize lattice%jm
c ---------------------
jm = jmglobal/lattice%ny
rm = jmglobal-lattice%ny*jm
lattice%jmax = jm
do n=0,lattice%ny-1
lattice%jm(n) = jm
if( n.le.rm-1 ) lattice%jm(n) = jm+1
lattice%jmax = max( lattice%jmax,lattice%jm(n) )
enddo
c Initialize relative PE address
c ------------------------------
lattice%pei = mod(myid,lattice%nx)
lattice%pej = myid/lattice%nx
c Initialize global index for local locations
c -------------------------------------------
call mymalloc ( lattice%iglobal,lattice%im(lattice%pei) )
call mymalloc ( lattice%jglobal,lattice%jm(lattice%pej) )
isum = 0
do n=0,lattice%nx-1
if ( n.eq.lattice%pei ) then
do i=1,lattice%im(n)
lattice%iglobal(i) = i+isum
enddo
endif
isum = isum + lattice%im(n)
enddo
jsum = 0
do m=0,lattice%ny-1
if ( m.eq.lattice%pej ) then
do j=1,lattice%jm(m)
lattice%jglobal(j) = j+jsum
enddo
endif
jsum = jsum + lattice%jm(m)
enddo
c Initialize local index for global locations
c -------------------------------------------
call mymalloc ( lattice%ilocal,lattice%imglobal )
call mymalloc ( lattice%jlocal,lattice%jmglobal )
n = 0
isum = 0
do i = 1,imglobal
if(i.gt.isum+lattice%im(n) ) then
isum = isum+lattice%im(n)
n = n + 1
endif
lattice%ilocal(i) = i-isum
enddo
m = 0
jsum = 0
do j = 1,jmglobal
if(j.gt.jsum+lattice%jm(m) ) then
jsum = jsum+lattice%jm(m)
m = m + 1
endif
lattice%jlocal(j) = j-jsum
enddo
c Initialize relative PE address for global i-j locations
c -------------------------------------------------------
call mymalloc ( lattice%peiglobal,imglobal )
call mymalloc ( lattice%pejglobal,jmglobal )
isum = 0
do n=0,lattice%nx-1
do i=1,lattice%im(n)
lattice%peiglobal( i+isum ) = n
enddo
isum = isum + lattice%im(n)
enddo
jsum = 0
do m=0,lattice%ny-1
do j=1,lattice%jm(m)
lattice%pejglobal( j+jsum ) = m
enddo
jsum = jsum + lattice%jm(m)
enddo
c Initialize Pole PE ghosts for each iglobal
c ------------------------------------------
isum = 0
do nx=0,lattice%nx-1
imx = lattice%im(nx)
do i=1,imx
ppe(i+isum) = nx
enddo
isum = isum + imx
enddo
do i=1,imglobal/2
lattice%ppeg(i) = ppe(i+imglobal/2)
lattice%ppeg(i+imglobal/2) = ppe(i)
enddo
c Allocate Lattice%img
c --------------------
if(.not.associated(lattice%img)) then
allocate ( lattice%img(0:nx-1,imglobal) )
do m=1,imglobal
do n=0,nx-1
lattice%img(n,m) = 0
enddo
enddo
else
m=size(lattice%img)
if(m.ne.nx*imglobal) then
print *, 'Allocated Lattice Size (',m,') does not match request (',nx*imglobal,') for lattice%img!'
call my_finalize
endif
endif
c Allocate Lattice%im0
c --------------------
if(.not.associated(lattice%im0)) then
allocate ( lattice%im0(0:nx-1,imglobal) )
do m=1,imglobal
do n=0,nx-1
lattice%im0(n,m) = 0
enddo
enddo
else
m=size(lattice%im0)
if(m.ne.nx*imglobal) then
print *, 'Allocated Lattice Size (',m,') does not match request (',nx*imglobal,') for lattice%im0!'
call my_finalize
endif
endif
c Determine Pole PE ghosts for each Processor
c -------------------------------------------
isum = 0
do nx=0,lattice%nx-1
imx = lattice%im(nx)
lattice%npeg(nx) = 1
lattice% im0(nx,1) = 1
img = 1
lattice%img (nx,lattice%npeg(nx)) = img
do i=2,imx
if( lattice%ppeg(i+isum) .eq. lattice%ppeg(i-1+isum) ) then
img = img + 1
lattice%img (nx,lattice%npeg(nx)) = img
else
lattice%npeg(nx) = lattice%npeg(nx) + 1
lattice%im0 (nx,lattice%npeg(nx)) = i
img = 1
lattice%img (nx,lattice%npeg(nx)) = img
endif
enddo
isum = isum + imx
enddo
c Print Lattice Characteristics
c -----------------------------
if( myid.eq.0 ) then
print *, 'Number of Processors in X: ',lattice%nx
print *, 'Number of Processors in Y: ',lattice%ny
print *
endif
do n=0,npes-1
if( n.eq.myid ) then
write(6,1000) myid,lattice%pei,lattice%pej,lattice%im(lattice%pei),lattice%jm(lattice%pej)
endif
call mpi_barrier (lattice%comm,ierror)
enddo
if( myid.eq.npes-1 ) print *
1000 format(1x,'absolute PE id: ',i3,' relative PE (i,j): (',i2,',',i2,') (im,jm) = (',i6,',',i6,')')
return
end
subroutine my_barrier (comm)
implicit none
integer comm,ierror
call mpi_barrier ( comm,ierror )
return
end
subroutine my_finalize 30
implicit none
integer ierror
call mpi_finalize (ierror )
return
end