program gtopo30,43 C ********************************************************************** C **** Program to create topography variance data to be used for **** C **** Gravity Wave Drag. GWD topography data is a smoothed **** C **** version of the residual between: **** C **** **** C **** Raw data ( hraw.30x30.data) **** C **** & GCM Model Resolved data (hmean.30x30.data) **** C **** **** C **** nbox .... Number of grid-points to use for average **** C **** **** C ********************************************************************** use dynamics_lattice_module implicit none type ( dynamics_lattice_type ) lattice integer imax, jmax, nbox parameter ( imax = 120*360 ) parameter ( jmax = 120*180 ) parameter ( nbox = 100 ) real, allocatable :: dum (:,:) real, allocatable :: halx (:,:) real, allocatable :: halo (:,:) real, allocatable :: hmean(:,:) real, allocatable :: hgrav(:,:) real, allocatable :: varimj(:,:) ! Variance in SW-NE Direction real, allocatable :: varipj(:,:) ! Variance in NW-SE Direction real, allocatable :: vimj(:) real, allocatable :: himj(:) integer, allocatable :: Nimj(:) real, allocatable :: vipj(:) real, allocatable :: hipj(:) integer, allocatable :: Nipj(:) real sumimj,sumipj integer i,j,ii,jj,k,imj,ipj integer im,jm 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 ( hmean(im,jm) ) allocate ( hgrav(im,jm) ) allocate ( varimj(im,jm) ) allocate ( varipj(im,jm) ) allocate ( vimj (-nbox+1:nbox-1) ) allocate ( himj (-nbox+1:nbox-1) ) allocate ( Nimj (-nbox+1:nbox-1) ) allocate ( vipj (-nbox+1:nbox-1) ) allocate ( hipj (-nbox+1:nbox-1) ) allocate ( Nipj (-nbox+1:nbox-1) ) C ********************************************************************** C **** Read and Distribute GTOPO30 Mean & GRAV Data **** C ********************************************************************** if( myid.eq.0 ) then allocate ( dum(imax,jmax) ) else allocate ( dum(1,1) ) endif if( myid.eq.0 ) then print *, 'Begin Reading MEAN Data' open (40,file='hmean.30x30sec.data',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) if( myid.eq.0 ) then print *, 'Begin Reading GRAV Data' open (40,file='hgrav.30x30sec.data',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,hgrav,lattice ) call printchar ('finished scatter of GRAV data',lattice) c Ghost GRAV Data c --------------- allocate ( halx(1-nbox/2:im+nbox/2,1:jm) ) call ghostx ( hgrav(1,1), halx(1-nbox/2,1), im,jm,1,nbox/2,lattice,'both' ) call printchar ('finished halo in x-direction',lattice) allocate ( halo(1-nbox/2:im+nbox/2,1-nbox/2:jm+nbox/2) ) call ghosty ( halx(1-nbox/2,1), halo(1-nbox/2,1-nbox/2), im+nbox,jm,1,0,1,nbox/2,lattice,'both' ) call printchar ('finished halo in y-direction',lattice) deallocate ( halx ) C ********************************************************************** C **** Compute (nbox X nbox) Variances **** C ********************************************************************** do j=1,jm do i=1,im c Water Point c ----------- if( hmean(i,j).eq.0.0 ) then varimj(i,j) = 0.0 varipj(i,j) = 0.0 else c Land Point: c ----------- c Compute Diagonal Means and Counters c ----------------------------------- himj(:) = 0.0 hipj(:) = 0.0 Nimj(:) = 0 Nipj(:) = 0 do jj = j-nbox/2,j+nbox/2 do ii = i-nbox/2,i+nbox/2 imj = (ii-i)-(jj-j) ipj = (ii-i)+(jj-j) himj(imj) = himj(imj) + halo(ii,jj) Nimj(imj) = Nimj(imj) + 1 hipj(ipj) = hipj(ipj) + halo(ii,jj) Nipj(ipj) = Nipj(ipj) + 1 enddo enddo himj(:) = himj(:) / Nimj(:) ! Diagonal Mean in SW-NE Direction hipj(:) = hipj(:) / Nipj(:) ! Diagonal Mean in NW-SE Direction c Compute Diagonal Variances c -------------------------- vimj(:) = 0.0 vipj(:) = 0.0 do jj = j-nbox/2,j+nbox/2 do ii = i-nbox/2,i+nbox/2 imj = (ii-i)-(jj-j) ipj = (ii-i)+(jj-j) vimj(imj) = vimj(imj) + ( halo(ii,jj)-himj(imj) )**2 vipj(ipj) = vipj(ipj) + ( halo(ii,jj)-hipj(ipj) )**2 enddo enddo vimj(:) = vimj(:) / ( Nimj(:)-1 ) ! Diagonal Variance in SW-NE Direction vipj(:) = vipj(:) / ( Nipj(:)-1 ) ! Diagonal Variance in NW-SE Direction c Average Diagonal Variances c -------------------------- sumimj = 0.0 sumipj = 0.0 do k = -nbox+1,nbox-1 sumimj = sumimj + Nimj(k)*vimj(k) sumipj = sumipj + Nipj(k)*vipj(k) enddo varimj(i,j) = sumimj / ( (nbox+1)*(nbox+1) ) varipj(i,j) = sumipj / ( (nbox+1)*(nbox+1) ) endif enddo if( myid.eq.0 ) print *, 'Finished J=',j enddo print *, 'Variance Calculation Finished, myid = ',myid C ********************************************************************** C **** Gather Data and Write Output **** C ********************************************************************** call gather ( dum,varimj,lattice ) call printchar ('finished gather of VARxy data',lattice) if( myid.eq.0 ) then print *, 'Begin Writing VARxy Data' open (40,file='hgrav_varxy.30x30sec.data',form='unformatted',access='sequential') do j=1,jmax write(40) (dum(i,j),i=1,imax) enddo close(40) print *, 'Finished Writing VARxy Data' endif call gather ( dum,varipj,lattice ) call printchar ('finished gather of VARyx data',lattice) if( myid.eq.0 ) then print *, 'Begin Writing VARyx Data' open (40,file='hgrav_varyx.30x30sec.data',form='unformatted',access='sequential') do j=1,jmax write(40) (dum(i,j),i=1,imax) enddo close(40) print *, 'Finished Writing VARyx 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