[cig-commits] r17265 - in seismo/3D/ADJOINT_TOMO/iterate_adj/pangu: model_vp_vs/src smooth/src
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Oct 13 07:12:09 PDT 2010
Author: danielpeter
Date: 2010-10-13 07:12:09 -0700 (Wed, 13 Oct 2010)
New Revision: 17265
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/exit_mpi.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/gll_library.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sub.f90
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe
Log:
adds missing files in directory iterate_adj/pangu/smooth/src/
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe 2010-10-13 13:14:21 UTC (rev 17264)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe 2010-10-13 14:12:09 UTC (rev 17265)
@@ -1,3 +1,4 @@
+FC=mpif90
FFLAGS=-O3 -assume byterecl
all: add_model_globe_iso \
@@ -5,13 +6,13 @@
add_model_globe_tiso_iso
add_model_globe_iso: add_model_globe_iso.f90
- mpif90 -o add_model_globe_iso $(FFLAGS) add_model_globe_iso.f90 gll_library.f90 exit_mpi.f90
+ $(FC) -o add_model_globe_iso $(FFLAGS) add_model_globe_iso.f90 gll_library.f90 exit_mpi.f90
add_model_globe_tiso: add_model_globe_tiso.f90
- mpif90 -o add_model_globe_tiso $(FFLAGS) add_model_globe_tiso.f90 gll_library.f90 exit_mpi.f90
+ $(FC) -o add_model_globe_tiso $(FFLAGS) add_model_globe_tiso.f90 gll_library.f90 exit_mpi.f90
add_model_globe_tiso_iso: add_model_globe_tiso_iso.f90
- mpif90 -o add_model_globe_tiso_iso $(FFLAGS) add_model_globe_tiso_iso.f90 gll_library.f90 exit_mpi.f90
+ $(FC) -o add_model_globe_tiso_iso $(FFLAGS) add_model_globe_tiso_iso.f90 gll_library.f90 exit_mpi.f90
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe 2010-10-13 13:14:21 UTC (rev 17264)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe 2010-10-13 14:12:09 UTC (rev 17265)
@@ -1,3 +1,4 @@
+FC=mpif90
FFLAGS=-O3 -assume byterecl #-g -traceback
all: smooth_sem_globe
@@ -3,5 +4,5 @@
smooth_sem_globe: smooth_sem_globe.f90
- mpif90 -o smooth_sem_globe $(FFLAGS) smooth_sem_globe.f90 smooth_sub.f90 exit_mpi.f90 gll_library_globe.f90
+ $(FC) -o smooth_sem_globe $(FFLAGS) smooth_sem_globe.f90 smooth_sub.f90 exit_mpi.f90 gll_library.f90
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/exit_mpi.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/exit_mpi.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/exit_mpi.f90 2010-10-13 14:12:09 UTC (rev 17265)
@@ -0,0 +1,61 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 3 . 6
+! --------------------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+! end the simulation and exit MPI
+
+! version with rank number printed in the error message
+ subroutine exit_MPI(myrank,error_msg)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+
+! identifier for error message file
+ integer, parameter :: IERROR = 30
+
+ integer myrank
+ character(len=*) error_msg
+
+ integer ier
+ character(len=80) outputname
+
+! write error message to screen
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI... proc ',myrank
+
+! write error message to file
+ write(outputname,"('OUTPUT_FILES/error_message',i6.6,'.txt')") myrank
+ open(unit=IERROR,file=trim(outputname),status='unknown')
+ write(IERROR,*) error_msg(1:len(error_msg))
+ write(IERROR,*) 'Error detected, aborting MPI... proc ',myrank
+ close(IERROR)
+
+
+! stop all the MPI processes, and exit
+! on some machines, MPI_FINALIZE needs to be called before MPI_ABORT
+ call MPI_FINALIZE(ier)
+ call MPI_ABORT(ier)
+ stop 'error, program ended in exit_MPI'
+
+ end subroutine exit_MPI
+
+!
+!----
+!
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/gll_library.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/gll_library.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/gll_library.f90 2010-10-13 14:12:09 UTC (rev 17265)
@@ -0,0 +1,529 @@
+
+!=======================================================================
+!
+! Library to compute the Gauss-Lobatto-Legendre points and weights
+! Based on Gauss-Lobatto routines from M.I.T.
+! Department of Mechanical Engineering
+!
+!=======================================================================
+
+ double precision function endw1(n,alpha,beta)
+
+ implicit none
+
+ integer n
+ double precision alpha,beta
+
+ double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0
+ double precision apb,f1,fint1,fint2,f2,di,abn,abnn,a1,a2,a3,f3
+ double precision, external :: gammaf
+ integer i
+
+ f3 = zero
+ apb = alpha+beta
+ if (n == 0) then
+ endw1 = zero
+ return
+ endif
+ f1 = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
+ f1 = f1*(apb+two)*two**(apb+two)/two
+ if (n == 1) then
+ endw1 = f1
+ return
+ endif
+ fint1 = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
+ fint1 = fint1*two**(apb+two)
+ fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
+ fint2 = fint2*two**(apb+three)
+ f2 = (-two*(beta+two)*fint1 + (apb+four)*fint2) * (apb+three)/four
+ if (n == 2) then
+ endw1 = f2
+ return
+ endif
+ do i=3,n
+ di = dble(i-1)
+ abn = alpha+beta+di
+ abnn = abn+di
+ a1 = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
+ a2 = (two*(alpha-beta))/(abnn*(abnn+two))
+ a3 = (two*(abn+one))/((abnn+two)*(abnn+one))
+ f3 = -(a2*f2+a1*f1)/a3
+ f1 = f2
+ f2 = f3
+ enddo
+ endw1 = f3
+
+ end function endw1
+
+!
+!=======================================================================
+!
+
+ double precision function endw2(n,alpha,beta)
+
+ implicit none
+
+ integer n
+ double precision alpha,beta
+
+ double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0
+ double precision apb,f1,fint1,fint2,f2,di,abn,abnn,a1,a2,a3,f3
+ double precision, external :: gammaf
+ integer i
+
+ apb = alpha+beta
+ f3 = zero
+ if (n == 0) then
+ endw2 = zero
+ return
+ endif
+ f1 = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
+ f1 = f1*(apb+two)*two**(apb+two)/two
+ if (n == 1) then
+ endw2 = f1
+ return
+ endif
+ fint1 = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
+ fint1 = fint1*two**(apb+two)
+ fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
+ fint2 = fint2*two**(apb+three)
+ f2 = (two*(alpha+two)*fint1 - (apb+four)*fint2) * (apb+three)/four
+ if (n == 2) then
+ endw2 = f2
+ return
+ endif
+ do i=3,n
+ di = dble(i-1)
+ abn = alpha+beta+di
+ abnn = abn+di
+ a1 = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
+ a2 = (two*(alpha-beta))/(abnn*(abnn+two))
+ a3 = (two*(abn+one))/((abnn+two)*(abnn+one))
+ f3 = -(a2*f2+a1*f1)/a3
+ f1 = f2
+ f2 = f3
+ enddo
+ endw2 = f3
+
+ end function endw2
+
+!
+!=======================================================================
+!
+
+ double precision function gammaf (x)
+
+ implicit none
+
+ double precision, parameter :: pi = 3.141592653589793d0
+
+ double precision x
+
+ double precision, parameter :: half=0.5d0,one=1.d0,two=2.d0
+
+ gammaf = one
+
+ if (x == -half) gammaf = -two*dsqrt(pi)
+ if (x == half) gammaf = dsqrt(pi)
+ if (x == one ) gammaf = one
+ if (x == two ) gammaf = one
+ if (x == 1.5d0) gammaf = dsqrt(pi)/2.d0
+ if (x == 2.5d0) gammaf = 1.5d0*dsqrt(pi)/2.d0
+ if (x == 3.5d0) gammaf = 2.5d0*1.5d0*dsqrt(pi)/2.d0
+ if (x == 3.d0 ) gammaf = 2.d0
+ if (x == 4.d0 ) gammaf = 6.d0
+ if (x == 5.d0 ) gammaf = 24.d0
+ if (x == 6.d0 ) gammaf = 120.d0
+
+ end function gammaf
+
+!
+!=====================================================================
+!
+
+ subroutine jacg (xjac,np,alpha,beta)
+
+!=======================================================================
+!
+! computes np Gauss points, which are the zeros of the
+! Jacobi polynomial with parameters alpha and beta
+!
+! .alpha = beta = 0.0 -> Legendre points
+! .alpha = beta = -0.5 -> Chebyshev points
+!
+!=======================================================================
+
+ implicit none
+
+ integer np
+ double precision alpha,beta
+ double precision xjac(np)
+
+ integer k,j,i,jmin,jm,n
+ double precision xlast,dth,x,x1,x2,recsum,delx,xmin,swap
+ double precision p,pd,pm1,pdm1,pm2,pdm2
+
+ integer, parameter :: K_MAX_ITER = 10
+ double precision, parameter :: zero = 0.d0, eps = 1.0d-12
+
+ pm1 = zero
+ pm2 = zero
+ pdm1 = zero
+ pdm2 = zero
+
+ xlast = 0.d0
+ n = np-1
+ dth = 4.d0*datan(1.d0)/(2.d0*dble(n)+2.d0)
+ p = 0.d0
+ pd = 0.d0
+ jmin = 0
+ do j=1,np
+ if(j == 1) then
+ x = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
+ else
+ x1 = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
+ x2 = xlast
+ x = (x1+x2)/2.d0
+ endif
+ do k=1,K_MAX_ITER
+ call jacobf (p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
+ recsum = 0.d0
+ jm = j-1
+ do i=1,jm
+ recsum = recsum+1.d0/(x-xjac(np-i+1))
+ enddo
+ delx = -p/(pd-recsum*p)
+ x = x+delx
+ if(abs(delx) < eps) goto 31
+ enddo
+ 31 continue
+ xjac(np-j+1) = x
+ xlast = x
+ enddo
+ do i=1,np
+ xmin = 2.d0
+ do j=i,np
+ if(xjac(j) < xmin) then
+ xmin = xjac(j)
+ jmin = j
+ endif
+ enddo
+ if(jmin /= i) then
+ swap = xjac(i)
+ xjac(i) = xjac(jmin)
+ xjac(jmin) = swap
+ endif
+ enddo
+
+ end subroutine jacg
+
+!
+!=====================================================================
+!
+
+ subroutine jacobf (poly,pder,polym1,pderm1,polym2,pderm2,n,alp,bet,x)
+
+!=======================================================================
+!
+! Computes the Jacobi polynomial of degree n and its derivative at x
+!
+!=======================================================================
+
+ implicit none
+
+ double precision poly,pder,polym1,pderm1,polym2,pderm2,alp,bet,x
+ integer n
+
+ double precision apb,polyl,pderl,dk,a1,a2,b3,a3,a4,polyn,pdern,psave,pdsave
+ integer k
+
+ apb = alp+bet
+ poly = 1.d0
+ pder = 0.d0
+ psave = 0.d0
+ pdsave = 0.d0
+
+ if (n == 0) return
+
+ polyl = poly
+ pderl = pder
+ poly = (alp-bet+(apb+2.d0)*x)/2.d0
+ pder = (apb+2.d0)/2.d0
+ if (n == 1) return
+
+ do k=2,n
+ dk = dble(k)
+ a1 = 2.d0*dk*(dk+apb)*(2.d0*dk+apb-2.d0)
+ a2 = (2.d0*dk+apb-1.d0)*(alp**2-bet**2)
+ b3 = (2.d0*dk+apb-2.d0)
+ a3 = b3*(b3+1.d0)*(b3+2.d0)
+ a4 = 2.d0*(dk+alp-1.d0)*(dk+bet-1.d0)*(2.d0*dk+apb)
+ polyn = ((a2+a3*x)*poly-a4*polyl)/a1
+ pdern = ((a2+a3*x)*pder-a4*pderl+a3*poly)/a1
+ psave = polyl
+ pdsave = pderl
+ polyl = poly
+ poly = polyn
+ pderl = pder
+ pder = pdern
+ enddo
+
+ polym1 = polyl
+ pderm1 = pderl
+ polym2 = psave
+ pderm2 = pdsave
+
+ end subroutine jacobf
+
+!
+!------------------------------------------------------------------------
+!
+
+ double precision FUNCTION PNDLEG (Z,N)
+
+!------------------------------------------------------------------------
+!
+! Compute the derivative of the Nth order Legendre polynomial at Z.
+! Based on the recursion formula for the Legendre polynomials.
+!
+!------------------------------------------------------------------------
+ implicit none
+
+ double precision z
+ integer n
+
+ double precision P1,P2,P1D,P2D,P3D,FK,P3
+ integer k
+
+ P1 = 1.d0
+ P2 = Z
+ P1D = 0.d0
+ P2D = 1.d0
+ P3D = 1.d0
+
+ do K = 1, N-1
+ FK = dble(K)
+ P3 = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
+ P3D = ((2.d0*FK+1.d0)*P2 + (2.d0*FK+1.d0)*Z*P2D - FK*P1D) / (FK+1.d0)
+ P1 = P2
+ P2 = P3
+ P1D = P2D
+ P2D = P3D
+ enddo
+
+ PNDLEG = P3D
+
+ end function pndleg
+
+!
+!------------------------------------------------------------------------
+!
+
+ double precision FUNCTION PNLEG (Z,N)
+
+!------------------------------------------------------------------------
+!
+! Compute the value of the Nth order Legendre polynomial at Z.
+! Based on the recursion formula for the Legendre polynomials.
+!
+!------------------------------------------------------------------------
+ implicit none
+
+ double precision z
+ integer n
+
+ double precision P1,P2,P3,FK
+ integer k
+
+ P1 = 1.d0
+ P2 = Z
+ P3 = P2
+
+ do K = 1, N-1
+ FK = dble(K)
+ P3 = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
+ P1 = P2
+ P2 = P3
+ enddo
+
+ PNLEG = P3
+
+ end function pnleg
+
+!
+!------------------------------------------------------------------------
+!
+
+ double precision function pnormj (n,alpha,beta)
+
+ implicit none
+
+ double precision alpha,beta
+ integer n
+
+ double precision one,two,dn,const,prod,dindx,frac
+ double precision, external :: gammaf
+ integer i
+
+ one = 1.d0
+ two = 2.d0
+ dn = dble(n)
+ const = alpha+beta+one
+
+ if (n <= 1) then
+ prod = gammaf(dn+alpha)*gammaf(dn+beta)
+ prod = prod/(gammaf(dn)*gammaf(dn+alpha+beta))
+ pnormj = prod * two**const/(two*dn+const)
+ return
+ endif
+
+ prod = gammaf(alpha+one)*gammaf(beta+one)
+ prod = prod/(two*(one+const)*gammaf(const+one))
+ prod = prod*(one+alpha)*(two+alpha)
+ prod = prod*(one+beta)*(two+beta)
+
+ do i=3,n
+ dindx = dble(i)
+ frac = (dindx+alpha)*(dindx+beta)/(dindx*(dindx+alpha+beta))
+ prod = prod*frac
+ enddo
+
+ pnormj = prod * two**const/(two*dn+const)
+
+ end function pnormj
+
+!
+!------------------------------------------------------------------------
+!
+
+ subroutine zwgjd(z,w,np,alpha,beta)
+
+!=======================================================================
+!
+! Z w g j d : Generate np Gauss-Jacobi points and weights
+! associated with Jacobi polynomial of degree n = np-1
+!
+! Note : Coefficients alpha and beta must be greater than -1.
+! ----
+!=======================================================================
+
+ implicit none
+
+ double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0
+
+ integer np
+ double precision z(np),w(np)
+ double precision alpha,beta
+
+ integer n,np1,np2,i
+ double precision p,pd,pm1,pdm1,pm2,pdm2
+ double precision apb,dnp1,dnp2,fac1,fac2,fac3,fnorm,rcoef
+ double precision, external :: gammaf,pnormj
+
+ pd = zero
+ pm1 = zero
+ pm2 = zero
+ pdm1 = zero
+ pdm2 = zero
+
+ n = np-1
+ apb = alpha+beta
+ p = zero
+ pdm1 = zero
+
+ if (np <= 0) stop 'minimum number of Gauss points is 1'
+
+ if ((alpha <= -one) .or. (beta <= -one)) stop 'alpha and beta must be greater than -1'
+
+ if (np == 1) then
+ z(1) = (beta-alpha)/(apb+two)
+ w(1) = gammaf(alpha+one)*gammaf(beta+one)/gammaf(apb+two) * two**(apb+one)
+ return
+ endif
+
+ call jacg(z,np,alpha,beta)
+
+ np1 = n+1
+ np2 = n+2
+ dnp1 = dble(np1)
+ dnp2 = dble(np2)
+ fac1 = dnp1+alpha+beta+one
+ fac2 = fac1+dnp1
+ fac3 = fac2+one
+ fnorm = pnormj(np1,alpha,beta)
+ rcoef = (fnorm*fac2*fac3)/(two*fac1*dnp2)
+ do i=1,np
+ call jacobf(p,pd,pm1,pdm1,pm2,pdm2,np2,alpha,beta,z(i))
+ w(i) = -rcoef/(p*pdm1)
+ enddo
+
+ end subroutine zwgjd
+
+!
+!------------------------------------------------------------------------
+!
+
+ subroutine zwgljd(z,w,np,alpha,beta)
+
+!=======================================================================
+!
+! Z w g l j d : Generate np Gauss-Lobatto-Jacobi points and the
+! ----------- weights associated with Jacobi polynomials of degree
+! n = np-1.
+!
+! Note : alpha and beta coefficients must be greater than -1.
+! Legendre polynomials are special case of Jacobi polynomials
+! just by setting alpha and beta to 0.
+!
+!=======================================================================
+
+ implicit none
+
+ double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0
+
+ integer np
+ double precision alpha,beta
+ double precision z(np), w(np)
+
+ integer n,nm1,i
+ double precision p,pd,pm1,pdm1,pm2,pdm2
+ double precision alpg,betg
+ double precision, external :: endw1,endw2
+
+ p = zero
+ pm1 = zero
+ pm2 = zero
+ pdm1 = zero
+ pdm2 = zero
+
+ n = np-1
+ nm1 = n-1
+ pd = zero
+
+ if (np <= 1) stop 'minimum number of Gauss-Lobatto points is 2'
+
+! with spectral elements, use at least 3 points
+ if (np <= 2) stop 'minimum number of Gauss-Lobatto points for the SEM is 3'
+
+ if ((alpha <= -one) .or. (beta <= -one)) stop 'alpha and beta must be greater than -1'
+
+ if (nm1 > 0) then
+ alpg = alpha+one
+ betg = beta+one
+ call zwgjd(z(2),w(2),nm1,alpg,betg)
+ endif
+
+ z(1) = - one
+ z(np) = one
+
+ do i=2,np-1
+ w(i) = w(i)/(one-z(i)**2)
+ enddo
+
+ call jacobf(p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(1))
+ w(1) = endw1(n,alpha,beta)/(two*pd)
+ call jacobf(p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(np))
+ w(np) = endw2(n,alpha,beta)/(two*pd)
+
+ end subroutine zwgljd
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sub.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sub.f90 2010-10-13 14:12:09 UTC (rev 17265)
@@ -0,0 +1,221 @@
+
+! --------------------------------------------------------------------------------
+
+subroutine get_all_eight_slices(ichunk,ixi,ieta,&
+ ileft,iright,ibot,itop, ilb,ilt,irb,irt,&
+ nproc_xi,nproc_eta)
+
+ implicit none
+
+ integer, intent(IN) :: ichunk,ixi,ieta,nproc_xi,nproc_eta
+
+ integer, intent(OUT) :: ileft,iright,ibot,itop,ilb,ilt,irb,irt
+ integer :: get_slice_number
+
+
+ integer :: ichunk_left, islice_xi_left, islice_eta_left, &
+ ichunk_right, islice_xi_right, islice_eta_right, &
+ ichunk_bot, islice_xi_bot, islice_eta_bot, &
+ ichunk_top, islice_xi_top, islice_eta_top, &
+ ileft0,iright0,ibot0,itop0, &
+ ichunk_left0, islice_xi_left0, islice_eta_left0, &
+ ichunk_right0, islice_xi_right0, islice_eta_right0, &
+ ichunk_bot0, islice_xi_bot0, islice_eta_bot0, &
+ ichunk_top0, islice_xi_top0, islice_eta_top0
+
+
+! get the first 4 immediate slices
+ call get_lrbt_slices(ichunk,ixi,ieta, &
+ ileft, ichunk_left, islice_xi_left, islice_eta_left, &
+ iright, ichunk_right, islice_xi_right, islice_eta_right, &
+ ibot, ichunk_bot, islice_xi_bot, islice_eta_bot, &
+ itop, ichunk_top, islice_xi_top, islice_eta_top, &
+ nproc_xi,nproc_eta)
+
+! get the 4 diagonal neighboring slices (actually 3 diagonal slices at the corners)
+ ilb = get_slice_number(ichunk,ixi-1,ieta-1,nproc_xi,nproc_eta)
+ ilt = get_slice_number(ichunk,ixi-1,ieta+1,nproc_xi,nproc_eta)
+ irb = get_slice_number(ichunk,ixi+1,ieta-1,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk,ixi+1,ieta+1,nproc_xi,nproc_eta)
+
+ if (ixi==0) then
+ call get_lrbt_slices(ichunk_left,islice_xi_left,islice_eta_left, &
+ ileft0, ichunk_left0, islice_xi_left0, islice_eta_left0, &
+ iright0, ichunk_right0, islice_xi_right0, islice_eta_right0, &
+ ibot0, ichunk_bot0, islice_xi_bot0, islice_eta_bot0, &
+ itop0, ichunk_top0, islice_xi_top0, islice_eta_top0, &
+ nproc_xi,nproc_eta)
+
+ if (ichunk == 0 .or. ichunk == 1 .or. ichunk == 3 .or. ichunk == 5) then
+ ilb = get_slice_number(ichunk_bot0,islice_xi_bot0,islice_eta_bot0,nproc_xi,nproc_eta)
+ ilt = get_slice_number(ichunk_top0,islice_xi_top0,islice_eta_top0,nproc_xi,nproc_eta)
+ else if (ichunk == 2) then
+ ilb = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ ilt = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ else
+ ilb = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ ilt = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ endif
+ endif
+
+ if (ixi==nproc_xi-1) then
+ call get_lrbt_slices(ichunk_right,islice_xi_right,islice_eta_right, &
+ ileft0, ichunk_left0, islice_xi_left0, islice_eta_left0, &
+ iright0, ichunk_right0, islice_xi_right0, islice_eta_right0, &
+ ibot0, ichunk_bot0, islice_xi_bot0, islice_eta_bot0, &
+ itop0, ichunk_top0, islice_xi_top0, islice_eta_top0, &
+ nproc_xi,nproc_eta)
+ if (ichunk == 0 .or. ichunk == 1 .or. ichunk == 3 .or. ichunk == 5) then
+ irb = get_slice_number(ichunk_bot0,islice_xi_bot0,islice_eta_bot0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_top0,islice_xi_top0,islice_eta_top0,nproc_xi,nproc_eta)
+ else if (ichunk == 2) then
+ irb = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ else
+ irb = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ endif
+ endif
+
+ if (ieta==0) then
+ call get_lrbt_slices(ichunk_bot,islice_xi_bot,islice_eta_bot, &
+ ileft0, ichunk_left0, islice_xi_left0, islice_eta_left0, &
+ iright0, ichunk_right0, islice_xi_right0, islice_eta_right0, &
+ ibot0, ichunk_bot0, islice_xi_bot0, islice_eta_bot0, &
+ itop0, ichunk_top0, islice_xi_top0, islice_eta_top0, &
+ nproc_xi,nproc_eta)
+ if (ichunk == 1 .or. ichunk == 2) then
+ ilb = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ irb = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ else if (ichunk == 3 .or. ichunk == 4) then
+ ilb = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ irb = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ else if (ichunk == 0) then
+ ilb = get_slice_number(ichunk_top0,islice_xi_top0,islice_eta_top0,nproc_xi,nproc_eta)
+ irb = get_slice_number(ichunk_bot0,islice_xi_bot0,islice_eta_bot0,nproc_xi,nproc_eta)
+ else
+ ilb = get_slice_number(ichunk_bot0,islice_xi_bot0,islice_eta_bot0,nproc_xi,nproc_eta)
+ irb = get_slice_number(ichunk_top0,islice_xi_top0,islice_eta_top0,nproc_xi,nproc_eta)
+ endif
+ endif
+
+ if (ieta==nproc_eta-1) then
+ call get_lrbt_slices(ichunk_top,islice_xi_top,islice_eta_top, &
+ ileft0, ichunk_left0, islice_xi_left0, islice_eta_left0, &
+ iright0, ichunk_right0, islice_xi_right0, islice_eta_right0, &
+ ibot0, ichunk_bot0, islice_xi_bot0, islice_eta_bot0, &
+ itop0, ichunk_top0, islice_xi_top0, islice_eta_top0, &
+ nproc_xi,nproc_eta)
+
+ if (ichunk == 1 .or. ichunk == 4) then
+ ilt = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ else if (ichunk == 2 .or. ichunk == 3) then
+ ilt = get_slice_number(ichunk_right0,islice_xi_right0,islice_eta_right0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_left0,islice_xi_left0,islice_eta_left0,nproc_xi,nproc_eta)
+ else if (ichunk == 0) then
+ ilt = get_slice_number(ichunk_bot0,islice_xi_bot0,islice_eta_bot0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_top0,islice_xi_top0,islice_eta_top0,nproc_xi,nproc_eta)
+ else
+ ilt = get_slice_number(ichunk_top0,islice_xi_top0,islice_eta_top0,nproc_xi,nproc_eta)
+ irt = get_slice_number(ichunk_bot0,islice_xi_bot0,islice_eta_bot0,nproc_xi,nproc_eta)
+ endif
+
+ endif
+
+end subroutine get_all_eight_slices
+
+!--------------------------------------------------------------------------------------------
+
+subroutine get_lrbt_slices(ichunk,ixi,ieta, &
+ ileft, ichunk_left, islice_xi_left, islice_eta_left, &
+ iright, ichunk_right, islice_xi_right, islice_eta_right, &
+ ibot, ichunk_bot, islice_xi_bot, islice_eta_bot, &
+ itop, ichunk_top, islice_xi_top, islice_eta_top, &
+ nproc_xi,nproc_eta)
+
+ implicit none
+
+ integer, intent(IN) :: ichunk, ixi, ieta, nproc_xi, nproc_eta
+ integer, intent(OUT) :: ileft, ichunk_left, islice_xi_left, islice_eta_left, &
+ iright, ichunk_right, islice_xi_right, islice_eta_right, &
+ ibot, ichunk_bot, islice_xi_bot, islice_eta_bot, &
+ itop, ichunk_top, islice_xi_top, islice_eta_top
+
+ integer, parameter :: NCHUNKS = 6
+
+ integer, dimension(NCHUNKS) :: chunk_left,chunk_right,chunk_bot,chunk_top, &
+ slice_xi_left,slice_eta_left,slice_xi_right,slice_eta_right, &
+ slice_xi_bot,slice_eta_bot,slice_xi_top,slice_eta_top
+ integer :: get_slice_number
+
+! set up mapping arrays -- assume chunk/slice number starts from 0
+ chunk_left(:) = (/2,6,6,1,6,4/) - 1
+ chunk_right(:) = (/4,1,1,6,1,2/) - 1
+ chunk_bot(:) = (/5,5,2,5,4,5/) - 1
+ chunk_top(:) = (/3,3,4,3,2,3/) - 1
+
+ slice_xi_left(:) = (/nproc_xi-1,nproc_xi-1,nproc_xi-1-ieta,nproc_xi-1,ieta,nproc_xi-1/)
+ slice_eta_left(:) = (/ieta,ieta,nproc_eta-1,ieta,0,ieta/)
+ slice_xi_right(:) = (/0,0,ieta,0,nproc_xi-1-ieta,0/)
+ slice_eta_right(:) = (/ieta,ieta,nproc_eta-1,ieta,0,ieta/)
+
+ slice_xi_bot(:) = (/nproc_xi-1,ixi,ixi,nproc_xi-1-ixi,nproc_xi-1-ixi,0/)
+ slice_eta_bot(:) = (/nproc_eta-1-ixi,nproc_eta-1,nproc_eta-1,0,0,ixi/)
+ slice_xi_top(:) = (/nproc_xi-1,ixi,nproc_xi-1-ixi,nproc_xi-1-ixi,ixi,0/)
+ slice_eta_top(:) = (/ixi,0,nproc_eta-1,nproc_eta-1,0,nproc_eta-1-ixi /)
+
+ ichunk_left = ichunk
+ ichunk_right = ichunk
+ ichunk_bot = ichunk
+ ichunk_top = ichunk
+
+ islice_xi_left = ixi-1
+ islice_eta_left = ieta
+ islice_xi_right = ixi+1
+ islice_eta_right = ieta
+
+ islice_xi_bot = ixi
+ islice_eta_bot = ieta-1
+ islice_xi_top = ixi
+ islice_eta_top = ieta+1
+
+ if (ixi == 0) then
+ ichunk_left=chunk_left(ichunk+1)
+ islice_xi_left=slice_xi_left(ichunk+1)
+ islice_eta_left=slice_eta_left(ichunk+1)
+ endif
+ if (ixi == nproc_xi - 1) then
+ ichunk_right=chunk_right(ichunk+1)
+ islice_xi_right=slice_xi_right(ichunk+1)
+ islice_eta_right=slice_eta_right(ichunk+1)
+ endif
+ if (ieta == 0) then
+ ichunk_bot=chunk_bot(ichunk+1)
+ islice_xi_bot=slice_xi_bot(ichunk+1)
+ islice_eta_bot=slice_eta_bot(ichunk+1)
+ endif
+ if (ieta == nproc_eta - 1) then
+ ichunk_top=chunk_top(ichunk+1)
+ islice_xi_top=slice_xi_top(ichunk+1)
+ islice_eta_top=slice_eta_top(ichunk+1)
+ endif
+
+ ileft = get_slice_number(ichunk_left,islice_xi_left,islice_eta_left,nproc_xi,nproc_eta)
+ iright = get_slice_number(ichunk_right,islice_xi_right,islice_eta_right,nproc_xi,nproc_eta)
+ ibot = get_slice_number(ichunk_bot,islice_xi_bot,islice_eta_bot,nproc_xi,nproc_eta)
+ itop = get_slice_number(ichunk_top,islice_xi_top,islice_eta_top,nproc_xi,nproc_eta)
+
+end subroutine get_lrbt_slices
+
+! ---------------------------------------------------------------------------------------
+
+integer function get_slice_number(ichunk,ixi,ieta,nproc_xi,nproc_eta)
+
+ implicit none
+
+ integer :: ichunk, ixi, ieta, nproc_xi, nproc_eta
+
+ get_slice_number = ichunk*nproc_xi*nproc_eta+ieta*nproc_xi+ixi
+
+ end function get_slice_number
More information about the CIG-COMMITS
mailing list