copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
***s/r adipc2_3d
*
subroutine adipc2_3d (g,h,s,czz,c2,dge,c5,nlim,nx,ny,nz,niter) 2,3
implicit none
*
integer nlim,nx,ny,nz,niter
real*8 g(nx,ny,nz), h(nx,ny,nz), dge, c2, c5,
$ czz(1-(niter-1):nx+(niter-1),1-(niter-1):ny+(niter-1),nz,3)
real s(1-(niter-1):nx+(niter-1),1-(niter-1):ny+(niter-1))
*
*AUTHOR Andre Robert. Sep 1987
*
*REVISION
* Evhen Yakimiw May 1988
* - version nested
* M. Tanguay May 1989
* - gal-chen sans montagnes
* Yves Chartier/Michel Desgagne Oct/Nov 1992
* - implicit none
* - structural documentation
* - working vectors memory allocation
* - in lining
* Michel Giguere/Michel Desgagne Dec 1992
* - solution now stable for small dt
* Guy Bergeron Sept. 1994
* - version 2D dans le plan XZ
* - automatisation du solveur
* Michel Desgagne Feb 2001
* - Transpose strategy for DM version
*
*OBJECT
* Effectue une iteration avec une variante tridimensionnelle
* du scheme adi (Peaceman-Rachford) utilise frequemment pour
* resoudre l'equation de poisson en deux dimensions.
*
* Seules les valeurs a l'exterieur a la zone frontiere sont evaluees
* (c.f. figure 3.2.3). Les valeurs des variables situees sur des
* points contenuent dans la zone de pilotage sont entierement dictees
* par le modele pilote.
*
*FILES
*
*ARGUMENTS
* NAMES I/O TYPE A/S DESCRIPTION
*
* g O R A valeur de qp a (t+dt)
* h I R A (-(A2/RT*)(dz/dt1)**2)
* s I R A Sbarxy*C1*(dz/dx)**2
* dz I R S distance entre deux niveaux du modele
* dt I R S (dt) modifie 0.5*(1+grepsi)*frtss
* dge I R S coeff. d'iteration de la methode ADI
* nx I I S dimension de la grille selon X
* ny I I S dimension de la grille selon Y
* nz I I S dimension de la grille selon Z
*
*IMPLICIT
#include "lcldim.cdk"
#include "nbcpu.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "transpose.cdk"
#include "partopo.cdk"
include 'mpif.h'
*
**
integer i,j,k,kp1,km1,offj,iter,rpn_comm_topo,err
real*8 dgel,dgex,con,conx,conz,one,two,four,
$ q(minx:maxx,miny:maxy,nz),
$ w1(gni,gnj),sg(gni,gnj),r(nh_maxx,nh_maxy,nz),
$ g1(nh_maxy,t1maxk,gni+npex),ax(nh_maxy,t1maxk,gni+npex),
$ bx(nh_maxy,t1maxk,gni+npex),cx(nh_maxy,t1maxk,gni+npex),
$ g2(t1maxk ,t2maxx,gnj+npey),ay(t1maxk ,t2maxx,gnj+npex),
$ by(t1maxk ,t2maxx,gnj+npex),cy(t1maxk ,t2maxx,gnj+npex),
$ az(nh_maxx,nh_maxy,nz),bz(nh_maxx,nh_maxy,nz),
$ cz(nh_maxx,nh_maxy,nz)
parameter (one = 1.0d0, two = 2.0d0, four = 4.0d0)
*
*----------------------------------------------------------------------
*
!$omp single
conx = one / ( dble(grdx) * dble(grdx) )
*
* Globalizing s into sg
*
do j=1,gnj
do i=1,gni
w1(i,j) = 0.
end do
end do
do j=1,ldnj
offj = gc_ld(3,myproc)-1
do i=1,ldni
w1(gc_ld(1,myproc)+i-1,offj+j) = s(i,j)*conx
end do
end do
if (numproc.gt.1) then
call MPI_ALLREDUCE (w1,sg,gni*gnj,MPI_DOUBLE_PRECISION,MPI_SUM,
$ MPI_COMM_WORLD,err)
else
do j=1,gnj
do i=1,gni
sg(i,j) = w1(i,j)
end do
end do
endif
*
dgel=dge
*
* g is pushed into q for halo exchange purposes
*
do k = 1, nz
do j = 1, ny
do i = 1, nx
q(i,j,k) = g(i,j,k)
end do
end do
end do
*
* Establishing ax and cx and filling out borders of bx
*
do j = 1, ldnj-north
do i = 1, gni-1
ax(j,1,i) = sg(i,gc_ld(3,myproc)+j-1)
cx(j,1,i) = sg(i,gc_ld(3,myproc)+j-1)
end do
ax(j,1,gni-1) = 0.0
cx(j,1,1 ) = 0.0
end do
do i = 1, gni-1
do j = ldnj-north+1, nh_maxy
ax(j,1,i) = 0.
bx(j,1,i) = 1.
cx(j,1,i) = 0.
end do
do k = 2, t1maxk
do j = 1, nh_maxy
ax(j,k,i) = ax(j,1,i)
cx(j,k,i) = cx(j,1,i)
end do
end do
end do
*
* Establishing ay and cy and filling out borders of by
*
do i = 1, t2n-teast
do j = 1, gnj-1
ay(1,i,j) = sg(t2n0+i-1,j)
cy(1,i,j) = sg(t2n0+i-1,j)
end do
ay(1,i,gnj-1) = 0.0
cy(1,i, 1) = 0.0
end do
do j = 1, gnj-1
do i = t2n-teast+1, t2maxx
ay(1,i,j) = 0.
by(1,i,j) = 1.
cy(1,i,j) = 0.
end do
do k = 2, t1maxk
do i = 1, t2maxx
ay(k,i,j) = ay(1,i,j)
cy(k,i,j) = cy(1,i,j)
end do
end do
end do
*
* Establishing az and cz and filling out borders of bz
*
do k = 1, nz
do j = 1, ldnj-north
do i = 1, ldni-east
az(i,j,k) = czz(i,j,k,3)
cz(i,j,k) = czz(i,j,k,1)
end do
end do
do j = ldnj-north+1, nh_maxy
do i = 1, ldni-east
az(i,j,k) = 0.
bz(i,j,k) = 1.
cz(i,j,k) = 0.
end do
end do
do i = ldni-east+1, nh_maxx
do j = 1, nh_maxy
az(i,j,k) = 0.
bz(i,j,k) = 1.
cz(i,j,k) = 0.
end do
end do
end do
*
* Perform nlim iterations
*
do 100 iter = 1, nlim
*
dgel = dgel*c5**two
dgex = dgel-c2
con = (two*dgel-c2)*dgex
*
do k = 1, nz
if (west.gt.0) then
do j = 1, ldnj-north
q(0,j,k) = q(1,j,k)
end do
endif
if (east.gt.0) then
do j = 1, ldnj-north
q(ldni,j,k) = q(ldni-1,j,k)
end do
endif
if (south.gt.0) then
do i = 1, ldni-east
q(i,0,k) = q(i,1,k)
end do
endif
if (north.gt.0) then
do i = 1, ldni-east
q(i,ldnj,k) = q(i,ldnj-1,k)
end do
endif
end do
*
call rpn_comm_xch_halon (q,minx,maxx,miny,maxy,ldni,ldnj,
$ nz,hx,hy,period_x,period_y,ldni,0,2)
*
do k = 1, nz
conz = conx
if(k.eq.gnk) conz = 0.
kp1 = min( nz,k+1 )
km1 = max( k-1, 1 )
do j = 1, ldnj-north
do i = 1, ldni-east
r(i,j,k) = conz * s(i,j) *
$ ( q(i+1,j ,k) + q(i-1,j ,k) +
$ q(i ,j+1,k) + q(i ,j-1,k) - four*q(i,j,k) ) +
$ czz(i,j,k,1) * q(i,j,km1) +
$ czz(i,j,k,2) * q(i,j,k ) +
$ czz(i,j,k,3) * q(i,j,kp1) - h(i,j,k)
end do
end do
end do
*
* Solving 3-diag matrix in x direction
*
call rpn_comm_transpose ( r,1,nh_maxx,gni,nh_maxy,1,t1maxk,nz,
$ g1, 1, 2 )
do j=ldnj-north+1,nh_maxy
do k=1,t1n
do i=1,gni-1
g1(j,k,i) = 0.
end do
end do
end do
*
do i = 1, gni-1
do j = 1, ldnj-north
bx(j,1,i) = dgex +ax(j,1,i) +cx(j,1,i)
end do
end do
do i = 2, gni-1
do j = 1, ldnj-north
bx(j,1,i) = bx(j,1,i)-cx(j,1,i)*ax(j,1,i-1)/bx(j,1,i-1)
end do
end do
do k = 2, t1maxk
do j = 1, nh_maxy
do i = 1, gni-1
bx(j,k,i) = bx(j,1,i)
end do
end do
end do
*
call tridiag (g1,ax,bx,cx,nh_maxy*t1maxk,nh_maxy*t1n,gni)
*
* Solving 3-diag matrix in y direction
*
call rpn_comm_transpose ( g1,1,nh_maxy,gnj,t1maxk,1,t2maxx,Gni,
$ g2, 2, 2 )
do k=t1n+1,t1maxk
do i=1,t2n-teast
do j=1,gnj-1
g2(k,i,j) = 0.
end do
end do
end do
*
do j = 1, gnj-1
do i = 1, t2n-teast
by(1,i,j) = dgex +ay(1,i,j) +cy(1,i,j)
end do
end do
do j = 2, gnj-1
do i = 1, t2n-teast
by(1,i,j) = by(1,i,j) - cy(1,i,j)*ay(1,i,j-1)/by(1,i,j-1)
end do
end do
do k = 2, t1maxk
do i = 1, t2maxx
do j = 1, gnj-1
by(k,i,j) = by(1,i,j)
end do
end do
end do
*
call tridiag (g2,ay,by,cy,t1maxk*t2maxx,t1maxk*(t2n-teast),gnj)
*
call rpn_comm_transpose (g1,1,nh_maxy,gnj,t1maxk ,1,t2maxx,Gni,
$ g2, -2, 2 )
call rpn_comm_transpose (r ,1,nh_maxx,gni,nh_maxy,1,t1maxk,nz,
$ g1, -1, 2 )
*
* Solving 3-diag matrix in z direction
*
do k = 1, nz
do j = 1, ldnj-north
do i = ldni-east+1,nh_maxx
r(i,j,k) = 0.
end do
end do
end do
do k = 1, nz
do j = 1, ldnj-north
do i = 1, ldni-east
bz(i,j,k) = -czz(i,j,k,2) + dgel
end do
end do
end do
do k = 2, nz
do j = 1, ldnj-north
do i = 1, ldni-east
bz(i,j,k)= bz(i,j,k)-cz(i,j,k)*az(i,j,k-1)/bz(i,j,k-1)
end do
end do
end do
*
call tridiag (r,az,bz,cz,nh_maxx*nh_maxy,nh_maxx*(ldnj-north),
$ nz+1)
*
do k = 1, nz
do j = 1, ldnj-north
do i = 1, ldni-east
q(i,j,k) = q(i,j,k) + con*r(i,j,k)
end do
end do
end do
*
100 continue
*
do k = 1, nz
do j = 1, ny
do i = 1, nx
g(i,j,k) = q(i,j,k)
end do
end do
end do
*
!$omp end single
*----------------------------------------------------------------------
return
end
*
subroutine tridiag (g,a,b,c,ns,n,gn) 3
implicit none
*
integer ns,n,gn
real*8 g(ns,*),a(ns,*),b(ns,*),c(ns,*)
*
integer i,j
*
do i=2,gn-1
do j=1,n
g(j,i) = g(j,i) + c(j,i) * g(j,i-1)/b(j,i-1)
end do
end do
do j=1,n
g(j,gn-1) = g(j,gn-1)/b(j,gn-1)
end do
do i=gn-2,1,-1
do j=1,n
g(j,i) = (g(j,i) + a(j,i) * g(j,i+1))/b(j,i)
end do
end do
*
return
end