copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
***s/r fwrd3d -- Compute so-called Q terms
*
subroutine fwrd3d ( xyzd, t1, dtm, dtp ) 1,1
implicit none
*
real t1(*)
real*8 xyzd(*),dtm,dtp
*
*
*AUTHORs M. Desgagne & C. Girard
*
*OBJECT
*
* *******************************************************************
* * *
* * ANDRE ROBERT SEMI-LAGRANGIAN SCHEME *
* * *
* * *
* * CALCULATES *
* * *
* * Q = P (t-dt,r-dr)+dtm*R (t,r-dr)+dtp*R (t,r) *
* * u u u u *
* * *
* * *
* * Q = P (t-dt,r-dr)+dtm*R (t,r-dr)+dtp*R (t,r) *
* * v v v v *
* * *
* * *
* * Q = P (t-dt,r-dr)+dtm*R (t,r-dr)+dtp*R (t,r) *
* * w w w w *
* * *
* * *
* * Q = P (t-dt,r-dr)+dtm*R (t,r-dr)+dtp*R (t,r) *
* * t t t t *
* * *
* * *
* * Q = P (t-dt,r-dr)+dtm*R (t,r-dr)+dtp*R (t,r) *
* * q q q q *
* * *
* * *
* * upstream interpolations are done by *
* * slag3d *
* * *
* *******************************************************************
*
* _________________________________________________________________
* | |
* | |
* | Example: Solving dQ/dt = R + S as follows |
* | |
* | |
* | |
* | Q(t+dt,r) = Q(t-dt,r-dr) + dtm*R(t,r-dr) + dtm*S(t-dt,r-dr) |
* | |
* | + dtp*R(t,r) + dtp*S(t-dt,r-dr) |
* | |
* | |
* | |
* | Grouping terms at (t-dt,r-dr) under the letter P : |
* | |
* | |
* | P(t-dt,r-dr) = Q(t-dt,r-dr) + (dtm+dtp)*S(t-dt,r-dr) |
* | |
* | |
* | we end up with: |
* | |
* | |
* | Q(t+dt,r) = P(t-dt,r-dr) + dtm*R(t,r-dr) + dtp*R(t,r) |
* | |
* | |
* | Grouping terms at r-dr under the letter Q" : |
* | |
* | |
* | Q"(r-dr) = P(t-dt,r-dr) + dtm*R(t,r-dr) |
* | |
* | |
* | we end up with: |
* | |
* | |
* | Q(t+dt,r) = Q"(r-dr) + dtp*R(t,r) |
* | |
* |_________________________________________________________________|
*
*
*ARGUMENTS
* _________________________________________________________________
* | | | |
* | NAME | DESCRIPTION | I/O |
* |---------|-------------------------------------------------|-----|
* | | | |
* | xyzd | displacements xd, yd, zd | i |
* | | | |
* | pptmp | Q term for variable pp (work) | |
* | uutmp | Q term for variable uu (work) | |
* | vvtmp | Q term for variable vv (work) | |
* | wwtmp | Q term for variable ww (work) | |
* | bbtmp | Q term for variable bb (work) | |
* | hutmp | Q term for variable hu (work) | |
* | trtmp | Q term for variable tr (work) | |
* | | | |
* | dtm | timestep lenght (-) | i |
* | dtp | timestep lenght (+) | i |
* | | | |
* |_________|_________________________________________________|_____|
*
*
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "dynmem.cdk"
#include "wrnmem.cdk"
#include "nbcpu.cdk"
#include "partopo.cdk"
*
**
integer i,j,k,n,id,jd,iff,jf,err
real pptmp,uutmp,vvtmp,wwtmp,bbtmp,hutmp,trtmp
pointer (papptmp, pptmp(minx:maxx,miny:maxy,*)),
$ (pauutmp, uutmp(minx:maxx,miny:maxy,*)),
$ (pavvtmp, vvtmp(minx:maxx,miny:maxy,*)),
$ (pawwtmp, wwtmp(minx:maxx,miny:maxy,*)),
$ (pabbtmp, bbtmp(minx:maxx,miny:maxy,*)),
$ (pahutmp, hutmp(minx:maxx,miny:maxy,*)),
$ (patrtmp, trtmp(minx:maxx,miny:maxy,gnk,*))
real*8 xdu, ydu, zdu, fxx, fyy, fzz
pointer (paxdu, xdu(minx:maxx,miny:maxy,*)),
$ (paydu, ydu(minx:maxx,miny:maxy,*)),
$ (pazdu, zdu(minx:maxx,miny:maxy,*)),
$ (pafxx, fxx(minx:maxx,miny:maxy,*)),
$ (pafyy, fyy(minx:maxx,miny:maxy,*)),
$ (pafzz, fzz(minx:maxx,miny:maxy,*))
include 'mpif.h'
real *8 time_i,time_f,time_it,time_ft
*----------------------------------------------------------------------
*
papptmp = loc(t1( 1))
pauutmp = loc(t1( dim3d+1))
pavvtmp = loc(t1(2*dim3d+1))
pawwtmp = loc(t1(3*dim3d+1))
pabbtmp = loc(t1(4*dim3d+1))
pahutmp = loc(t1(5*dim3d+1))
if (ntr.ge.1) patrtmp = loc(t1(6*dim3d+1))
paxdu = loc(xyzd(3*dim3d+1))
paydu = loc(xyzd(4*dim3d+1))
pazdu = loc(xyzd(5*dim3d+1))
pafxx = loc(xyzd(6*dim3d+1))
pafyy = loc(xyzd(7*dim3d+1))
pafzz = loc(xyzd(8*dim3d+1))
*
* A) Compute the part of Q terms to be evaluated at upwind positions
*
* Q' = P (t-dt,r)+dtm*R (t,r)
* x x x
*
* upwind-point range
*
id =1-hx*west
jd =1-hy*south
iff=ldni+(hx-1)*east
jf =ldnj+(hy-1)*north
*
!$omp do
do k=1,gnk
if(k.ne.gnk) then
do j=jd,jf
do i=id+west,iff
uutmp(i,j,k) = uup(i,j,k) + dtm*uur(i,j,k)
end do
end do
do j=jd+south,jf
do i=id,iff
vvtmp(i,j,k) = vvp(i,j,k) + dtm*vvr(i,j,k)
end do
end do
endif
do j=jd,jf
do i=id,iff
wwtmp(i,j,k) = wwp(i,j,k) + dtm*wwr(i,j,k)
bbtmp(i,j,k) = bbp(i,j,k) + dtm*bbr(i,j,k)
pptmp(i,j,k) = ppp(i,j,k) + dtm*ppr(i,j,k)
hutmp(i,j,k) = hup(i,j,k)
end do
end do
do n=1,ntr
do j=jd,jf
do i=id,iff
trtmp(i,j,k,n) = trp(i,j,k,n)
end do
end do
end do
end do
!$omp enddo
*
* B) Interpolate this part at upwind positions
* and perform its Lagrangian displacement
*
* from Q'(t-dt,r) to Q"(t-dt,r-dr)
* x x
*
!$omp single
call rpn_comm_xch_halon (xyzd,minx,maxx,miny,maxy,ldni,ldnj,
$ 3*gnk,hx,hy,period_x,period_y,ldni,0,2)
call rpn_comm_xch_halo (pptmp,minx,maxx,miny,maxy,ldni,ldnj,
$ 6*gnk,hx,hy,period_x,period_y,ldni,0)
if (ntr.gt.0)
$ call rpn_comm_xch_halo (trtmp,minx,maxx,miny,maxy,ldni,ldnj,
$ ntr*gnk,hx,hy,period_x,period_y,ldni,0)
!$omp end single
*
call slag3d ( ppp(minx,miny,1),uup,vvp,wwp,bbp,hup,trp,
$ pptmp,uutmp,vvtmp,wwtmp,bbtmp,hutmp,trtmp,
$ xyzd, xyzd(dim3d+1),xyzd(2*dim3d+1),
$ xdu,ydu,zdu,fxx,fyy,fzz,
$ minx,maxx,miny,maxy )
*
*
* C) Add the final contribution to the Q terms
*
* Q = Q" + dtp*R
* x x x
*
* at destination points
*
id =1
jd =1
iff=ldni-east
jf =ldnj-north
*
!$omp do
do k=1,gnk
if(k.ne.gnk) then
do j=jd,jf
do i=id+west,iff
uup(i,j,k) = uup(i,j,k) + dtp*uur(i,j,k)
end do
end do
do j=jd+south,jf
do i=id,iff
vvp(i,j,k) = vvp(i,j,k) + dtp*vvr(i,j,k)
end do
end do
endif
do j=jd,jf
do i=id,iff
wwp(i,j,k) = wwp(i,j,k) + dtp*wwr(i,j,k)
bbp(i,j,k) = bbp(i,j,k) + dtp*bbr(i,j,k)
ppp(i,j,k) = ppp(i,j,k) + dtp*ppr(i,j,k)
end do
end do
end do
!$omp enddo
*
*----------------------------------------------------------------------
return
end