!
!--------------------------------------------------------------------------
!
SUBROUTINE updbc( mptr,samlvl,tmwght ) 1,10
!
! Provide boundary values (conditons) for grid MPTR by
! interpolating from course grid or from overlapping fine grids
! at the same level depending on the value of samlvl..
!
! MPTR is the grid needing boundary values by interpolation
! samlvl is a logical, if true - interp from same level
! (i.e., from overlapping grid)
! false - interp from next level
! TMWGHT is the weight for the interpolated value, with
!
! phi(t)=TMWGHT*phi(interp)+(1-TMWGHT)*phi(t-dt(receive)) (1)
!
! where phi(interp) is the interpolated value from source grid
! at t+dt(source).
!
! This routine resets phi(t-dt(receive)) at the boundaries so that
! formula (1) is satisfied, then linear extrapolation of phi
! from its value at t and t-dt(receive) will set phi(t+n*dt(receive))
! equal to the interpolated value, where n=dt(source)/dt(receive)..
!
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'agrialloc.inc'
INCLUDE 'grddsc.inc'
INCLUDE 'agricpu.inc'
INCLUDE 'agricst.inc'
DIMENSION msrcs(30),ibcvar(nxy2d+nxyz3d),ibcsta(nxy2d+nxyz3d)
LOGICAL :: samlvl,smstg2
!
! here we assume the number of source grids never exceeds 30
!
IF( intrpodr == 2 ) THEN
nwhts=9
ELSE
nwhts=4
END IF
level = node(4,mptr)
IF(level > lfine.OR.level == 0) RETURN
cpu0 = f_cputime
()
!
! Get the possible source grids for the b.c.'s of grid mptr.
!
numscs=0
lchk = level
IF(.NOT.samlvl) lchk = lchk-1
!
! do 20 lchk=max(1,level-1),level
!
mgrid=lstart(lchk)
15 IF(mgrid == 0) GO TO 20
numscs=numscs+1
msrcs(numscs)=mgrid
mgrid=node(10,mgrid)
GO TO 15
20 CONTINUE
5 CONTINUE
nx=node(5,mptr)
ny=node(6,mptr)
!
! --------------------------------------------------------------------
! calculate b.c.'s for 2-d scalars
! --------------------------------------------------------------------
!
ndim=2
!
! initialize b.c. calculating status with value 1 for each scalar
!
DO i=1,nxy2d
ibcsta(i)=1
END DO
100 CONTINUE
!
! search for variables that need b.c.'s and group them according to
! their staggering
!
nbcvar=0
DO ibc=1,nxy2d
!
! look for 1st 2-d scalar that needs b.c. calculation
!
IF(ibdxy(1,ibc) == 1.AND.ibcsta(ibc) == 1.AND. ibdxy(2,ibc) == 0)THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ibc
ibcsta(ibc)=0
DO ivar=ibc+1,nxy2d
!
! list additonal variables with matching staggering to that of ibc.
!
IF(ibdxy(1,ivar) == 1.AND.ibcsta(ivar) == 1.AND. &
ibdxy(2,ivar) == 0.AND.smstg2(ivar,ibc,2))THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ivar
ibcsta(ivar)=0
END IF
END DO
END IF
END DO
!
! if the b.c.'s for all 2-d scalars are calculated goto next step
!
IF(nbcvar == 0)GO TO 200
!
! perform interpolations for 2-d scalar variables listed in ibcvar(nbcvar).
!
mxpts=2*(nx+ny)
iptr1=igetsp(2*mxpts)
iptr2=igetsp(2*mxpts)
iptr3=igetsp(2*mxpts)
iptr4=igetsp( mxpts)
iptr5=igetsp(nwhts*mxpts)
iptr6=igetsp(2*mxpts)
iptr7=igetsp(2*mxpts)
!
CALL itpbsl
(mptr,msrcs,numscs,a(iptr1),a(iptr2),a(iptr3), &
a(iptr4),a(iptr5),a(iptr6),a(iptr7), &
mxpts,nwhts,ibcvar,nbcvar,ndim,samlvl,tmwght)
!
! restore tem space
CALL resett
!
! go back to check if any more variable needs b.c. calculation
!
GO TO 100
200 CONTINUE
!
! --------------------------------------------------------------------
! calculate b.c.'s for 3-d scalars
! --------------------------------------------------------------------
!
ndim=3
!
! initialize b.c. calculating status with value 1 for each vari.
!
DO i=1,nxyz3d
ibcsta(i)=1
END DO
!
! search for variables that need b.c.'s and group them according to
! their staggering
!
300 CONTINUE
nbcvar=0
DO ibc=1,nxyz3d
!
! look for 1st 3-d scalar that needs b.c. calculation
!
IF(ibdxyz(1,ibc) == 1.AND.ibcsta(ibc) == 1.AND. ibdxyz(2,ibc) == 0)THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ibc
ibcsta(ibc)=0
DO ivar=ibc+1,nxyz3d
!
! list additonal variables with matching staggering to that of ibc.
!
IF(ibdxyz(1,ivar) == 1.AND.ibcsta(ivar) == 1.AND. &
ibdxyz(2,ivar) == 0.AND.smstg2(ivar,ibc,3))THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ivar
ibcsta(ivar)=0
END IF
END DO
END IF
END DO
!
! if the b.c.'s for all 3-d scalars are calculated goto next step
!
IF(nbcvar == 0) GO TO 400
!
! perform interpolations for 3-d scalars listed in ibcvar(nbcvar).
!
mxpts=2*(nx+ny)
iptr1=igetsp(2*mxpts)
iptr2=igetsp(2*mxpts)
iptr3=igetsp(2*mxpts)
iptr4=igetsp( mxpts)
iptr5=igetsp(nwhts*mxpts)
iptr6=igetsp(2*mxpts)
iptr7=igetsp(2*mxpts)
!
CALL itpbsl
(mptr,msrcs,numscs,a(iptr1),a(iptr2),a(iptr3), &
a(iptr4),a(iptr5),a(iptr6),a(iptr7), &
mxpts,nwhts,ibcvar,nbcvar,ndim,samlvl,tmwght)
!
! restore tem space
!
CALL resett
!
! go back to check if any more 3-d scalars needs b.c. calculation
!
GO TO 300
400 CONTINUE
!
!-------------------------------------------------------------------
! calculate b.c.'s for 2-d vectors
! --------------------------------------------------------------------
!
ndim=2
!
! initialize b.c. calculating status with value 1 for each vari.
!
DO i=1,nxy2d
ibcsta(i)=1
END DO
500 CONTINUE
!
! search for 2-d vectors that need b.c.'s
!
nbcvar=0
DO ibc=1,nxy2d
!
! look for 1st 2-d vector that needs b.c. calculation
!
IF(ibdxy(1,ibc) == 1.AND.ibdxy(2,ibc) > 0.AND. ibcsta(ibc) == 1) THEN
ibc2=ibdxy(2,ibc)
ibcvar(nbcvar+1)=ibc
ibcvar(nbcvar+2)=ibc2
nbcvar=nbcvar+2
ibcsta(ibc )=0
ibcsta(ibc2)=0
EXIT
END IF
END DO
67 CONTINUE
!
! if the b.c.'s for all 2-d vectors are calculated goto next step
!
IF(nbcvar == 0)GO TO 600
!
! perform interpolations for 2-d vectors listed in ibcvar(nbcvar).
!
mxpts=2*(nx+ny)
iptr1=igetsp(4*mxpts)
iptr2=igetsp(4*mxpts)
iptr3=igetsp(4*mxpts)
iptr4=igetsp(2*mxpts)
iptr5=igetsp(2*nwhts*mxpts)
iptr6=igetsp(4*mxpts)
iptr7=igetsp(4*mxpts)
!
CALL itpbvt
(mptr,msrcs,numscs,a(iptr1),a(iptr2),a(iptr3), &
a(iptr4),a(iptr5),a(iptr6),a(iptr7), &
mxpts,nwhts,ibcvar,nbcvar,ndim,samlvl,tmwght)
!
! restore tem space
!
CALL resett
!
! go back to check if any more 2-d vector needs b.c. calculation
!
GO TO 500
600 CONTINUE
!
!------------------------------------------------------------------
! calculate b.c.'s for 3-d vectors
!------------------------------------------------------------------
!
ndim=3
!
! initialize b.c. calculating status with value 1 for each vari.
!
DO i=1,nxyz3d
ibcsta(i)=1
END DO
!
700 CONTINUE
!
! search for variables that need b.c.'s and group them according to
! their staggering
!
nbcvar=0
DO ibc=1,nxyz3d
!
! look for 1st 3-d vector that needs b.c. calculation
!
IF(ibdxyz(1,ibc) == 1.AND.ibdxyz(2,ibc) > 0.AND. ibcsta(ibc) == 1) THEN
ibc2=ibdxyz(2,ibc)
ibcvar(nbcvar+1)=ibc
ibcvar(nbcvar+2)=ibc2
nbcvar=nbcvar+2
ibcsta(ibc )=0
ibcsta(ibc2)=0
EXIT
END IF
END DO
92 CONTINUE
!
! if the b.c.'s for all 3-d vectors are calculated goto next step
!
IF(nbcvar == 0)GO TO 800
!
! perform interpolations for 3-d vectors listed in ibcvar(nbcvar).
!
mxpts=2*(nx+ny)
iptr1=igetsp(4*mxpts)
iptr2=igetsp(4*mxpts)
iptr3=igetsp(4*mxpts)
iptr4=igetsp(2*mxpts)
iptr5=igetsp(2*nwhts*mxpts)
iptr6=igetsp(4*mxpts)
iptr7=igetsp(4*mxpts)
!
IF( verbose6 ) WRITE(6,'('' ITPBVT CALLED FOR 3-D VECTORS'')')
!
CALL itpbvt
(mptr,msrcs,numscs,a(iptr1),a(iptr2),a(iptr3), &
a(iptr4),a(iptr5),a(iptr6),a(iptr7), &
mxpts,nwhts,ibcvar,nbcvar,ndim,samlvl,tmwght)
!
IF( verbose6 ) THEN
DO ib=1,nbcvar,2
WRITE(6,'(2i5)') ibcvar(ib),ibcvar(ib+1)
END DO
END IF
!
! restore tem space
!
CALL resett
!
! go back to check if any more 3-d vector needs b.c. calculation
!
GO TO 700
800 CONTINUE
cpu_bndcint = cpu_bndcint + f_cputime
() - cpu0
RETURN
END SUBROUTINE updbc
!
!--------------------------------------------------------------------------
!
SUBROUTINE itpbsl(mrec,msrcs,numscs,ptsrec, & 2,10
irec,isrc,igsrc,whts,pts,irecp,mxpts,nwhts, &
ibcvar,nbcvar,ndim,samlvl,tmwght)
!
DIMENSION msrcs(numscs),ptsrec(2,mxpts),irec(2,mxpts), &
isrc(2,mxpts),igsrc(mxpts),whts(nwhts,mxpts), &
pts(2,mxpts),irecp(2,mxpts), &
ibcvar(nbcvar)
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'agrialloc.inc'
INCLUDE 'grddsc.inc'
LOGICAL :: samlvl
nw1=nint( SQRT(FLOAT(nwhts)) )
nw2=nw1
!
! Calculate and store the coor's of p-points of the receive
! grid in array ptsrec with their indecies stored in irec.
!
IF(nbcvar == 0) RETURN
CALL chkdim
(ndim,'ITPBSL')
nx=node(5,mrec)
ny=node(6,mrec)
cosr=rnode(21,mrec)
sinr=rnode(22,mrec)
dxr=rnode(9 ,mrec)
dyr=rnode(10,mrec)
xor=rnode(1,mrec)
yor=rnode(2,mrec)
ivar=ibcvar(1)
IF(ndim == 2) THEN
xshift=stgxy(1,ivar)
yshift=stgxy(2,ivar)
nx1=nx-idmxy(1,ivar)
ny1=ny-idmxy(2,ivar)
ELSE
xshift=stgxyz(1,ivar)
yshift=stgxyz(2,ivar)
nx1=nx-idmxyz(1,ivar)
ny1=ny-idmxyz(2,ivar)
END IF
xc=xshift*dxr*cosr
ys=yshift*dyr*sinr
xs=xshift*dxr*sinr
yc=yshift*dyr*cosr
xorp=xor+xc-ys
yorp=yor+xs+yc
!
! To calculate the absolute coordinates of the model boundary points
!
! Attention:
!
! Here we assume that in x direction the boundaries are at i=1 and nx1.
! and at j=1 and ny1 in y direction.
! If in a particular solver (model), the boundaries for certain
! variables are defined at points other the above, the use needs to
! change the following loop accordingly (i.e. change ist, jst, iend
! and jend). As a result, the part in routine UPDSCL that calculates
! the points for interior points updating also need change.
!
ist=1
jst=1
iend=nx1
jend=ny1
DO j=jst,jend,jend-jst
jj=(j-jst)/(jend-jst)*(iend-ist+1)
DO i=ist,iend
ii=jj+i-ist+1
x=(i-1)*dxr
y=(j-1)*dyr
ptsrec(1,ii) = xorp+x*cosr-y*sinr
ptsrec(2,ii) = yorp+x*sinr+y*cosr
irec (1,ii) = i
irec (2,ii) = j
END DO
END DO
numpts=MAX(0, 2*(iend-ist+1))
DO i=ist,iend,iend-ist
ii=(i-ist)/(iend-ist)*(jend-jst-1)
DO j=jst+1,jend-1
jj=ii+j-jst+numpts
x=(i-1)*dxr
y=(j-1)*dyr
ptsrec(1,jj) = xorp+x*cosr-y*sinr
ptsrec(2,jj) = yorp+x*sinr+y*cosr
irec (1,jj) = i
irec (2,jj) = j
END DO
END DO
numpts=numpts+MAX(0,2*(jend-jst-1))
npts0=numpts
!
! Find the source points for receive grid points ptsrec
!
IF(samlvl) THEN
ierfl=0
iptcmp=1
ELSE
ierfl=1
iptcmp=1
END IF
CALL getsrc
(mrec,msrcs,numscs,xshift,yshift,ptsrec,numpts, &
irec,igsrc,ierfl,iptcmp)
!
! if no point can draw boudary value from listed source grids, skip through
!
IF(.NOT.samlvl) THEN
IF(numpts /= npts0) THEN
WRITE(6,'(i4,'' BOUNDARY POINTS DID NOT FIND SOURCE FOR B.C.S'' &
& )') npts0-numpts
WRITE(6,'(''IT WAS FOR '',I1,''-D SCALAR NO. '',I4)')ndim,ivar
WRITE(6,'(''JOB STOPPED IN ITPBSL.'')')
STOP
END IF
END IF
!
! Loop through all source grids
!
DO ks=1,numscs
msrc=msrcs(ks)
!
! sort rec. grid points according to their source grid
!
npts=0
DO np=1,numpts
IF(igsrc(np) == msrc) THEN
npts=npts+1
pts(1,npts)=ptsrec(1,np)
pts(2,npts)=ptsrec(2,np)
irecp(1,npts)=irec(1,np)
irecp(2,npts)=irec(2,np)
END IF
END DO
IF(npts == 0) CYCLE
nxr=node(5,mrec)
nyr=node(6,mrec)
nzr=node(14,mrec)
nxs=node(5,msrc)
nys=node(6,msrc)
nzs=node(14,msrc)
!
! Calculate weights of interpolation for points pts.
!
CALL calwht
(mrec,msrc,pts,isrc,whts,npts,nwhts,ivar,ndim)
!
! Interpolate for the boundary values for variables in list ibcvar
!
DO i=1,nbcvar
ivar=ibcvar(i)
IF(ndim == 2) THEN
n3rd=1
irptr=igtnxy(mrec,ivar,1)
istmdt=igtnxy(mrec,ibdxy(3,ivar),1)
isptr=igtnxy(msrc,ivar,1)
ELSE
n3rd=nzs
irptr=igtxyz(mrec,ivar,1)
istmdt=igtxyz(mrec,ibdxyz(3,ivar),1)
isptr=igtxyz(msrc,ivar,1)
END IF
CALL intpst1
( nxr,nyr,nxs,nys,n3rd, &
irecp,isrc,whts,npts,nw1,nw2)
IF(ndim == 2) THEN
CALL retnxy
(mrec,ivar,1,irptr,.true.)
CALL retnxy
(mrec,ibdxy(3,ivar),1,istmdt,.false.)
CALL retnxy
(msrc,ivar,1,isptr,.false.)
ELSE
CALL retxyz
(mrec,ivar,1,irptr,.false.)
CALL retxyz
(mrec,ibdxyz(3,ivar),1,istmdt,.true.)
! call retxyz(mrec,ivar,1,irptr,.true.)
! call retxyz(mrec,ibdxyz(3,ivar),1,istmdt,.false.)
CALL retxyz
(msrc,ivar,1,isptr,.false.)
END IF
END DO
!
END DO
500 CONTINUE
RETURN
END SUBROUTINE itpbsl
!
!--------------------------------------------------------------------------
!
SUBROUTINE itpbvt(mrec,msrcs,numscs,ptsrec, & 2,15
irec,isrc,igsrc,whts,pts,irecp,mxpts,nwhts, &
ibcvar,nbcvar,ndim,samlvl,tmwght)
!
DIMENSION msrcs(numscs),ptsrec(2,mxpts,2),irec(2,mxpts,2), &
isrc(2,mxpts,2),igsrc(mxpts,2),whts(nwhts,mxpts,2), &
pts(2,mxpts,2),irecp(2,mxpts,2),ibcvar(nbcvar), &
numpts(2),npts(2)
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'agrialloc.inc'
INCLUDE 'grddsc.inc'
LOGICAL :: samstg, smstg2, samlvl
nw1=nint( SQRT(FLOAT(nwhts)) )
nw2=nw1
!
IF(nbcvar == 0) RETURN
CALL chkdim
(ndim,'ITPBVT')
nxr=node(5,mrec)
nyr=node(6,mrec)
nzr=node(14,mrec)
cosr=rnode(21,mrec)
sinr=rnode(22,mrec)
dxr=rnode(9 ,mrec)
dyr=rnode(10,mrec)
xor=rnode(1,mrec)
yor=rnode(2,mrec)
!
! To loop through all vectors, which has two components stored consecutively
! in nbcvar. The 1st compt. has to be the one in x direction.
!
DO ivtr=1,nbcvar,2
ivar1=ibcvar(ivtr)
ivar2=ibcvar(ivtr+1)
samstg=smstg2(ivar1,ivar2,ndim)
IF(ndim == 2) THEN
xshft = MAX( stgxy(1,ivar1),stgxy(1,ivar2) )
yshft = MAX( stgxy(2,ivar1),stgxy(2,ivar2) )
ELSE IF(ndim == 3)THEN
xshft = MAX( stgxyz(1,ivar1),stgxyz(1,ivar2) )
yshft = MAX( stgxyz(2,ivar1),stgxyz(2,ivar2) )
END IF
!
! Find the points on rec. grid for both compnts. and the source grid for
! them.If the two compnts have the same staggering, they have the
! same coordinate points and the source grid, therefore those for the
! 2nd compnt. do not need to be recalculated.
!
DO icmpr=1,2
!
! icmpr=1 for the 1st component of vectors
! =2 for the 2nd component of vectors
!
IF(icmpr == 2.AND.samstg)THEN
numpts(2)=numpts(1)
DO ip=1,numpts(2)
ptsrec(1,ip,2)=ptsrec(1,ip,1)
ptsrec(2,ip,2)=ptsrec(2,ip,1)
irec(1,ip,2)=irec(1,ip,1)
irec(2,ip,2)=irec(2,ip,1)
igsrc(ip,2)=igsrc(ip,1)
END DO
CYCLE
END IF
ivar=ibcvar(ivtr+icmpr-1)
IF(ndim == 2) THEN
xshift=stgxy(1,ivar)
yshift=stgxy(2,ivar)
nx1=nxr-idmxy(1,ivar)
ny1=nyr-idmxy(2,ivar)
ELSE IF(ndim == 3)THEN
xshift=stgxyz(1,ivar)
yshift=stgxyz(2,ivar)
nx1=nxr-idmxyz(1,ivar)
ny1=nyr-idmxyz(2,ivar)
END IF
xc=xshift*dxr*cosr
ys=yshift*dyr*sinr
xs=xshift*dxr*sinr
yc=yshift*dyr*cosr
xorp=xor+xc-ys
yorp=yor+xs+yc
!
! To calculate the absolute coordinates of the model boundary points
!
! Attention:
!
! Here we assume that in x direction the boundaries are at i=1 and nx1.
! and at j=1 and ny1 in y direction.
! If in a particular solver (model), the boundaries for certain
! variables are defined at points other the above, the use needs to
! change the following loop accordingly (i.e. change ist, jst, iend
! and jend). As a result, the part in routine UPDVTR that calculates
! the points for interior points updating also need change.
!
ist=1
jst=1
iend=nx1
jend=ny1
DO j=jst,jend,jend-jst
jj=(j-jst)/(jend-jst)*(iend-ist+1)
DO i=ist,iend
ii=jj+i-ist+1
x=(i-1)*dxr
y=(j-1)*dyr
ptsrec(1,ii,icmpr) = xorp+x*cosr-y*sinr
ptsrec(2,ii,icmpr) = yorp+x*sinr+y*cosr
irec (1,ii,icmpr) = i
irec (2,ii,icmpr) = j
END DO
END DO
numpts(icmpr)=MAX(0,2*(iend-ist+1))
DO i=ist,iend,iend-ist
ii=(i-ist)/(iend-ist)*(jend-jst-1)
DO j=jst+1,jend-1
jj=ii+j-jst+numpts(icmpr)
x=(i-1)*dxr
y=(j-1)*dyr
ptsrec(1,jj,icmpr) = xorp+x*cosr-y*sinr
ptsrec(2,jj,icmpr) = yorp+x*sinr+y*cosr
irec (1,jj,icmpr) = i
irec (2,jj,icmpr) = j
END DO
END DO
numpts(icmpr)=numpts(icmpr)+MAX(0,2*(jend-jst-1))
npts0=numpts(icmpr)
!
! Find the source grid for receive grid points ptsrec
!
IF(samlvl) THEN
ierfl=0
iptcmp=1
ELSE
ierfl=1
iptcmp=1
END IF
CALL getsrc
(mrec,msrcs,numscs,xshft,yshft,ptsrec(1,1,icmpr), &
numpts(icmpr),irec(1,1,icmpr),igsrc(1,icmpr),ierfl,iptcmp)
!
IF(.NOT.samlvl) THEN
IF(numpts(icmpr) /= npts0) THEN
WRITE(6,'(i4,'' BOUNDARY POINTS DID NOT FIND SOURCE FOR B.C.S'' &
& )') npts0-numpts(icmpr)
WRITE(6,'(''IT WAS FOR '',I1,''-D VECTOR COMPNT. '',2I4)') &
ndim,ivar
WRITE(6,'(''JOB STOPPED IN ITPBVT.'')')
STOP
END IF
END IF
END DO
!
! find storage pointers for the receive grid variables (2 compnts of vector).
! uppack them if neccesary.
!
IF(ndim == 2) THEN
irptr1=igtnxy(mrec,ivar1,1)
irptr2=igtnxy(mrec,ivar2,1)
mrptr1=igtnxy(mrec,ibdxy(3,ivar1),1)
mrptr2=igtnxy(mrec,ibdxy(3,ivar2),1)
ELSE
irptr1=igtxyz(mrec,ivar1,1)
irptr2=igtxyz(mrec,ivar2,1)
mrptr1=igtxyz(mrec,ibdxyz(3,ivar1),1)
mrptr2=igtxyz(mrec,ibdxyz(3,ivar2),1)
END IF
!
! Loop through all source grids
!
DO ks=1,numscs
msrc=msrcs(ks)
nxs=node(5,msrc)
nys=node(6,msrc)
nzs=node(14,msrc)
coss=rnode(21,msrc)
sins=rnode(22,msrc)
cossr= coss*cosr+sins*sinr
sinsr= coss*sinr-sins*cosr
DO icmps=1,2
DO icmpr=1,2
IF(icmpr == 2.AND.samstg)THEN
npts(2)=npts(1)
DO ip=1,npts(2)
pts (1,ip,2)=pts (1,ip,1)
pts (2,ip,2)=pts (2,ip,1)
isrc (1,ip,2)=isrc (1,ip,1)
isrc (2,ip,2)=isrc (2,ip,1)
irecp(1,ip,2)=irecp(1,ip,1)
irecp(2,ip,2)=irecp(2,ip,1)
END DO
DO iw=1,nwhts
DO ip=1,npts(2)
whts(iw,ip,2)=whts(iw,ip,1)
END DO
END DO
CYCLE
END IF
!
! sort rec. grid points according to their source grid
!
npt=0
DO np=1,numpts(icmpr)
IF(igsrc(np,icmpr) == msrc) THEN
npt=npt+1
pts(1,npt,icmpr)=ptsrec(1,np,icmpr)
pts(2,npt,icmpr)=ptsrec(2,np,icmpr)
irecp(1,npt,icmpr)=irec(1,np,icmpr)
irecp(2,npt,icmpr)=irec(2,np,icmpr)
END IF
END DO
npts(icmpr)=npt
IF(npts(icmpr) == 0) CYCLE
!
! Calculate weights of interpolation for points pts.
!
ivsrc=ibcvar(ivtr+icmps-1)
CALL calwht
(mrec,msrc,pts(1,1,icmpr),isrc(1,1,icmpr), &
whts(1,1,icmpr),npts(icmpr),nwhts,ivsrc,ndim)
END DO
IF(npts(1) == 0.AND.npts(2) == 0) CYCLE
!
! find pointer for one of the components (icmps) of the source grid
! vector. If it is packed, do uppacking for this variable.
!
ivsrc=ibcvar(ivtr+icmps-1)
IF(ndim == 2) THEN
n3rd=1
isptr=igtnxy(msrc,ivsrc,1)
ELSE
n3rd=nzs
isptr=igtxyz(msrc,ivsrc,1)
END IF
IF(icmps == 1)THEN
proj1=cossr
proj2=-sinsr
ELSE
proj1=sinsr
proj2=cossr
END IF
icmpr=1
IF(npts(icmpr) /= 0) THEN
CALL intpvt1
( nxr,nyr,nxs,nys,n3rd, &
irecp(1,1,icmpr),isrc(1,1,icmpr),whts(1,1,icmpr), &
npts(icmpr),nw1,nw2,icmps,proj1)
END IF
icmpr=2
IF(npts(icmpr) /= 0) THEN
CALL intpvt1
( nxr,nyr,nxs,nys,n3rd, &
irecp(1,1,icmpr),isrc(1,1,icmpr),whts(1,1,icmpr), &
npts(icmpr),nw1,nw2,icmps,proj2)
END IF
IF(ndim == 2) THEN
CALL retnxy
(msrc,ivsrc,1,isptr,.false.)
ELSE
CALL retxyz
(msrc,ivsrc,1,isptr,.false.)
END IF
END DO
!
END DO
IF(ndim == 2) THEN
CALL retnxy
(mrec,ivar1,1,irptr1,.false.)
CALL retnxy
(mrec,ivar2,1,irptr2,.false.)
CALL retnxy
(mrec,ibdxy(3,ivar1),1,mrptr1,.true.)
CALL retnxy
(mrec,ibdxy(3,ivar2),1,mrptr2,.true.)
ELSE
CALL retxyz
(mrec,ivar1,1,irptr1,.false.)
CALL retxyz
(mrec,ivar2,1,irptr2,.false.)
CALL retxyz
(mrec,ibdxyz(3,ivar1),1,mrptr1,.true.)
CALL retxyz
(mrec,ibdxyz(3,ivar2),1,mrptr2,.true.)
END IF
END DO
RETURN
END SUBROUTINE itpbvt
!
!--------------------------------------------------------------------------
!
SUBROUTINE intpst(varrec,varsrc,varmdt,tmwght, &
nxr,nyr,nxs,nys,nz, &
irec,isrc,whts,numpts,nw1,nw2)
DIMENSION varrec(nxr,nyr,nz),varsrc(nxs,nys,nz), &
varmdt(nxr,nyr,nz), &
irec(2,numpts),isrc(2,numpts),whts(nw1,nw2,numpts)
!
itest = 0
IF(itest == 1) THEN
WRITE(6,'('' DATA DUMP IN INTRPS '',8I6)') &
nxr,nyr,nxs,nys,nz,numpts,nw1,nw2
DO ip=1,20
WRITE(7,1114) irec(1,ip),irec(2,ip),isrc(1,ip),isrc(2,ip), &
whts(1,1,ip),whts(2,1,ip),whts(1,2,ip),whts(2,2,ip)
END DO
1114 FORMAT(1X,4I7,4E10.3)
END IF
!
! set weights for temporal linear interp
!
w1 = tmwght
w2 = 1.-w1
!
IF(nw1 == 2.AND.nw2 == 2) THEN
!
! The interpolation calculations are all independent of each other,
! they can be vecotrized in inner loop and parellelized in outer loop.
!
nzswap=2
!
! if nz if small, swap the loop order
!
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varrec(ir,jr,ks)=w1*(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np)) &
+w2*varmdt(ir,jr,ks)
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varrec(ir,jr,ks)=w1*(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np)) &
+w2*varmdt(ir,jr,ks)
END DO
END DO
END IF
ELSE
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varrec(ir,jr,ks)= w1*( &
varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np) ) &
+w2*varmdt(ir,jr,ks)
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varrec(ir,jr,ks)= w1*( &
varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np) ) &
+w2*varmdt(ir,jr,ks)
END DO
END DO
END IF
END IF
RETURN
END SUBROUTINE intpst
!
!--------------------------------------------------------------------------
!
SUBROUTINE intpvt(varrec,varsrc,varmdt,tmwght, &
nxr,nyr,nxs,nys,nz, &
irec,isrc,whts,numpts,nw1,nw2,mode,proj)
DIMENSION varrec(nxr,nyr,nz),varsrc(nxs,nys,nz), &
varmdt(nxr,nyr,nz), &
irec(2,numpts),isrc(2,numpts),whts(nw1,nw2,numpts)
!
itest = 0
IF(itest == 1) THEN
WRITE(6,'('' DATA DUMP IN INTRPV '',8I6)') &
nxr,nyr,nxs,nys,nz,numpts,nw1,nw2
DO ip=1,20
WRITE(7,1114) irec(1,ip),irec(2,ip),isrc(1,ip),isrc(2,ip), &
whts(1,1,ip),whts(2,1,ip),whts(1,2,ip),whts(2,2,ip)
END DO
1114 FORMAT(1X,4I7,4E10.3)
END IF
iaccum=0
IF(mode == 2) iaccum=1
w1 = tmwght
w2 = 1.-w1
IF(nw1 == 2.AND.nw2 == 2) THEN
!
! The interpolation calculations are all independent of each other,
! they can be vecotrized in inner loop and parellelized in outer loop.
!
nzswap=2
!
! if nz if small, swap the loop order
!
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varrec(ir,jr,ks)=varrec(ir,jr,ks)*iaccum &
+(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np))*proj*w1 &
+w2*iaccum*varmdt(ir,jr,ks)
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varrec(ir,jr,ks)=varrec(ir,jr,ks)*iaccum &
+(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np))*proj*w1 &
+w2*iaccum*varmdt(ir,jr,ks)
END DO
END DO
END IF
ELSE
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varrec(ir,jr,ks)=varrec(ir,jr,ks)*iaccum &
+(varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np))*proj*w1 &
+w2*iaccum*varmdt(ir,jr,ks)
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varrec(ir,jr,ks)=varrec(ir,jr,ks)*iaccum &
+(varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np))*proj*w1 &
+w2*iaccum*varmdt(ir,jr,ks)
END DO
END DO
END IF
END IF
RETURN
END SUBROUTINE intpvt
!
!--------------------------------------------------------------------------
!
SUBROUTINE intpst1(varrec,varsrc,varmdt,tmwght, & 1
nxr,nyr,nxs,nys,nz, &
irec,isrc,whts,numpts,nw1,nw2)
DIMENSION varrec(nxr,nyr,nz),varsrc(nxs,nys,nz), &
varmdt(nxr,nyr,nz), &
irec(2,numpts),isrc(2,numpts),whts(nw1,nw2,numpts)
!
itest = 0
IF(itest == 1) THEN
WRITE(6,'('' DATA DUMP IN INTRPS '',8I6)') &
nxr,nyr,nxs,nys,nz,numpts,nw1,nw2
DO ip=1,20
WRITE(7,1114) irec(1,ip),irec(2,ip),isrc(1,ip),isrc(2,ip), &
whts(1,1,ip),whts(2,1,ip),whts(1,2,ip),whts(2,2,ip)
END DO
1114 FORMAT(1X,4I7,4E10.3)
END IF
!
! set weights for temporal linear interp
!
w1 = tmwght
w2 = 1.-w1
!
IF(nw1 == 2.AND.nw2 == 2) THEN
!
! The interpolation calculations are all independent of each other,
! they can be vecotrized in inner loop and parellelized in outer loop.
!
nzswap=2
!
! if nz if small, swap the loop order
!
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varmdt(ir,jr,ks)=(varrec(ir,jr,ks)-w1* &
(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np)))/w2
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varmdt(ir,jr,ks)=(varrec(ir,jr,ks)-w1* &
(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np)))/w2
END DO
END DO
END IF
ELSE
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varmdt(ir,jr,ks)=(varrec(ir,jr,ks)-w1*( &
varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np)))/w2
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varmdt(ir,jr,ks)=(varrec(ir,jr,ks)-w1*( &
varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np)))/w2
END DO
END DO
END IF
END IF
RETURN
END SUBROUTINE intpst1
!
!--------------------------------------------------------------------------
!
SUBROUTINE intpvt1(varrec,varsrc,varmdt,tmwght, & 2
nxr,nyr,nxs,nys,nz, &
irec,isrc,whts,numpts,nw1,nw2,mode,proj)
DIMENSION varrec(nxr,nyr,nz),varsrc(nxs,nys,nz), &
varmdt(nxr,nyr,nz), &
irec(2,numpts),isrc(2,numpts),whts(nw1,nw2,numpts)
!
itest = 0
IF(itest == 1) THEN
WRITE(6,'('' DATA DUMP IN INTRPV '',8I6)') &
nxr,nyr,nxs,nys,nz,numpts,nw1,nw2
DO ip=1,20
WRITE(7,1114) irec(1,ip),irec(2,ip),isrc(1,ip),isrc(2,ip), &
whts(1,1,ip),whts(2,1,ip),whts(1,2,ip),whts(2,2,ip)
END DO
1114 FORMAT(1X,4I7,4E10.3)
END IF
iaccum=0
IF(mode == 2) iaccum=1
w1 = tmwght
w2 = 1.-w1
IF(nw1 == 2.AND.nw2 == 2) THEN
!
! The interpolation calculations are all independent of each other,
! they can be vecotrized in inner loop and parellelized in outer loop.
!
nzswap=2
!
! if nz if small, swap the loop order
!
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varmdt(ir,jr,ks)=varmdt(ir,jr,ks)*iaccum &
+(varrec(ir,jr,ks)*iaccum &
-(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np))*proj*w1)/w2
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varmdt(ir,jr,ks)=varmdt(ir,jr,ks)*iaccum &
+(varrec(ir,jr,ks)*iaccum &
-(varsrc(is,js,ks)*whts(1,1,np) &
+varsrc(is+1,js,ks)*whts(2,1,np) &
+varsrc(is,js+1,ks)*whts(1,2,np) &
+varsrc(is+1,js+1,ks)*whts(2,2,np))*proj*w1)/w2
END DO
END DO
END IF
ELSE
IF(nz <= nzswap)THEN
!fpp$ nodepchk
DO ks=1,nz
!DIR$ IVDEP
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
varmdt(ir,jr,ks)=varmdt(ir,jr,ks)*iaccum &
+(varrec(ir,jr,ks)*iaccum &
-(varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np))*proj*w1)/w2
END DO
END DO
ELSE
!fpp$ nodepchk
DO np=1,numpts
ir = irec(1,np)
jr = irec(2,np)
is = isrc(1,np)
js = isrc(2,np)
!DIR$ IVDEP
DO ks=1,nz
varmdt(ir,jr,ks)=varmdt(ir,jr,ks)*iaccum &
+(varrec(ir,jr,ks)*iaccum &
-(varsrc(is-1,js-1,ks)*whts(1,1,np) &
+varsrc(is ,js-1,ks)*whts(2,1,np) &
+varsrc(is+1,js-1,ks)*whts(3,1,np) &
+varsrc(is-1,js ,ks)*whts(1,2,np) &
+varsrc(is ,js ,ks)*whts(2,2,np) &
+varsrc(is+1,js ,ks)*whts(3,2,np) &
+varsrc(is-1,js+1,ks)*whts(1,3,np) &
+varsrc(is ,js+1,ks)*whts(2,3,np) &
+varsrc(is+1,js+1,ks)*whts(3,3,np))*proj*w1)/w2
END DO
END DO
END IF
END IF
RETURN
END SUBROUTINE intpvt1
!
!--------------------------------------------------------------------------
!
SUBROUTINE exchbc( mptr1,mptr2 ) 1,8
!
! Exchange boundary values (conditons) by interpolation bwtween
! two overlapping grids at the same level. mptr1 and mptr2 are the
! grid number.
!
! Calling this once instead of calling updbc(mptr) twice so that each
! grid do not need to be unpacked and packed twice.
!
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'agrialloc.inc'
INCLUDE 'grddsc.inc'
INCLUDE 'agricst.inc'
DIMENSION ibcvar(nxy2d+nxyz3d),ibcsta(nxy2d+nxyz3d)
LOGICAL :: smstg2
!
IF( intrpodr == 2 ) THEN
nwhts=9
ELSE
nwhts=4
END IF
level1 = node(4,mptr1)
IF(level1 > lfine.OR.level1 == 0) GO TO 999
level2 = node(4,mptr2)
IF(level2 > lfine.OR.level2 == 0) GO TO 998
IF(level1 /= level2) GO TO 997
nx1=node(5,mptr1)
ny1=node(6,mptr1)
nx2=node(5,mptr2)
ny2=node(6,mptr2)
!
! --------------------------------------------------------------------
! calculate b.c.'s for 2-d scalars
! --------------------------------------------------------------------
!
ndim=2
!
! initialize b.c. calculating status with value 1 for each scalar
!
DO i=1,nxy2d
ibcsta(i)=1
END DO
100 CONTINUE
!
! search for variables that need b.c.'s and group them according to
! their staggering
!
nbcvar=0
DO ibc=1,nxy2d
!
! look for 1st 2-d scalar that needs b.c. calculation
!
IF(ibdxy(1,ibc) == 1.AND.ibcsta(ibc) == 1.AND. ibdxy(2,ibc) == 0)THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ibc
ibcsta(ibc)=0
DO ivar=ibc+1,nxy2d
!
! list additonal variables with matching staggering to that of ibc.
!
IF(ibdxy(1,ivar) == 1.AND.ibcsta(ivar) == 1.AND. &
ibdxy(2,ivar) == 0.AND.smstg2(ivar,ibc,ndim))THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ivar
ibcsta(ivar)=0
END IF
END DO
END IF
END DO
!
! if the b.c.'s for all 2-d scalars are calculated goto next step
!
IF(nbcvar == 0)GO TO 200
!
! perform interpolations for 2-d scalar variables listed in ibcvar(nbcvar).
!
mxpts=2*MAX( nx1+ny1,nx2+ny2 )
iptr1=igetsp(2*mxpts)
iptr2=igetsp(2*mxpts*2)
iptr3=igetsp(2*mxpts*2)
iptr4=igetsp( mxpts)
iptr5=igetsp(nwhts*mxpts*2)
iptr6=igetsp(2*mxpts)
!
CALL itpbs2
(mptr1,mptr2,a(iptr1),a(iptr2),a(iptr3),a(iptr4), &
a(iptr5),a(iptr6),mxpts,nwhts,ibcvar,nbcvar,ndim)
!
! restore tem space
CALL resett
!
! go back to check if any more variable needs b.c. calculation
!
GO TO 100
200 CONTINUE
!
! --------------------------------------------------------------------
! calculate b.c.'s for 3-d scalars
! --------------------------------------------------------------------
!
ndim=3
!
! initialize b.c. calculating status with value 1 for each vari.
!
DO i=1,nxyz3d
ibcsta(i)=1
END DO
!
! search for variables that need b.c.'s and group them according to
! their staggering
!
300 CONTINUE
nbcvar=0
DO ibc=1,nxyz3d
!
! look for 1st 3-d scalar that needs b.c. calculation
!
IF(ibdxyz(1,ibc) == 1.AND.ibcsta(ibc) == 1.AND. ibdxyz(2,ibc) == 0)THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ibc
ibcsta(ibc)=0
DO ivar=ibc+1,nxyz3d
!
! list additonal variables with matching staggering to that of ibc.
!
IF(ibdxyz(1,ivar) == 1.AND.ibcsta(ivar) == 1.AND. &
ibdxyz(2,ivar) == 0.AND.smstg2(ivar,ibc,ndim))THEN
nbcvar=nbcvar+1
ibcvar(nbcvar)=ivar
ibcsta(ivar)=0
END IF
END DO
END IF
END DO
!
! if the b.c.'s for all 3-d scalars are calculated goto next step
!
IF(nbcvar == 0) GO TO 400
!
! perform interpolations for 3-d scalars listed in ibcvar(nbcvar).
!
mxpts=2*MAX( nx1+ny1,nx2+ny2 )
iptr1=igetsp(2*mxpts)
iptr2=igetsp(2*mxpts*2)
iptr3=igetsp(2*mxpts*2)
iptr4=igetsp( mxpts)
iptr5=igetsp(nwhts*mxpts*2)
iptr6=igetsp(2*mxpts)
!
CALL itpbs2
(mptr1,mptr2,a(iptr1),a(iptr2),a(iptr3),a(iptr4), &
a(iptr5),a(iptr6),mxpts,nwhts,ibcvar,nbcvar,ndim)
!
! restore tem space
!
CALL resett
!
! go back to check if any more 3-d scalars needs b.c. calculation
!
GO TO 300
400 CONTINUE
!
!-------------------------------------------------------------------
! calculate b.c.'s for 2-d vectors
! --------------------------------------------------------------------
!
ndim=2
!
! initialize b.c. calculating status with value 1 for each vari.
!
DO i=1,nxy2d
ibcsta(i)=1
END DO
500 CONTINUE
!
! search for 2-d vectors that need b.c.'s
!
nbcvar=0
DO ibc=1,nxy2d
!
! look for 1st 2-d vector that needs b.c. calculation
!
IF(ibdxy(1,ibc) == 1.AND.ibdxy(2,ibc) > 0.AND. ibcsta(ibc) == 1) THEN
ibc2=ibdxy(2,ibc)
ibcvar(nbcvar+1)=ibc
ibcvar(nbcvar+2)=ibc2
nbcvar=nbcvar+2
ibcsta(ibc )=0
ibcsta(ibc2)=0
EXIT
END IF
END DO
67 CONTINUE
!
! if the b.c.'s for all 2-d vectors are calculated goto next step
!
IF(nbcvar == 0)GO TO 600
!
! perform interpolations for 2-d vectors listed in ibcvar(nbcvar).
!
mxpts=2*MAX( nx1+ny1,nx2+ny2 )
iptr1=igetsp(4*mxpts)
iptr2=igetsp(4*mxpts)
iptr3=igetsp(4*mxpts)
iptr4=igetsp(2*mxpts)
iptr5=igetsp(nwhts*mxpts*2)
iptr6=igetsp(4*mxpts)
!
CALL itpbv2
(mptr1,mptr2,a(iptr1),a(iptr2),a(iptr3),a(iptr4), &
a(iptr5),a(iptr6),mxpts,nwhts,ibcvar,nbcvar,ndim)
!
! restore tem space
!
CALL resett
!
! go back to check if any more 2-d vector needs b.c. calculation
!
GO TO 500
600 CONTINUE
!
!------------------------------------------------------------------
! calculate b.c.'s for 3-d vectors
!------------------------------------------------------------------
!
ndim=3
!
! initialize b.c. calculating status with value 1 for each vari.
!
DO i=1,nxyz3d
ibcsta(i)=1
END DO
!
700 CONTINUE
!
! search for variables that need b.c.'s and group them according to
! their staggering
!
nbcvar=0
DO ibc=1,nxyz3d
!
! look for 1st 3-d vector that needs b.c. calculation
!
IF(ibdxyz(1,ibc) == 1.AND.ibdxyz(2,ibc) > 0.AND. ibcsta(ibc) == 1) THEN
ibc2=ibdxyz(2,ibc)
ibcvar(nbcvar+1)=ibc
ibcvar(nbcvar+2)=ibc2
nbcvar=nbcvar+2
ibcsta(ibc )=0
ibcsta(ibc2)=0
EXIT
END IF
END DO
92 CONTINUE
!
! if the b.c.'s for all 3-d vectors are calculated goto next step
!
IF(nbcvar == 0)GO TO 800
!
! perform interpolations for 3-d vectors listed in ibcvar(nbcvar).
!
mxpts=2*MAX( nx1+ny1,nx2+ny2 )
iptr1=igetsp(4*mxpts)
iptr2=igetsp(4*mxpts)
iptr3=igetsp(4*mxpts)
iptr4=igetsp(2*mxpts)
iptr5=igetsp(nwhts*mxpts*2)
iptr6=igetsp(4*mxpts)
!
CALL itpbv2
(mptr1,mptr2,a(iptr1),a(iptr2),a(iptr3),a(iptr4), &
a(iptr5),a(iptr6),mxpts,nwhts,ibcvar,nbcvar,ndim)
!
! restore tem space
!
CALL resett
!
! go back to check if any more 3-d vector needs b.c. calculation
!
GO TO 700
800 CONTINUE
RETURN
999 WRITE(6,'(''ERROR IN EXCHBC: THE LEVEL OF INPUT GRID'', &
& '' <1 OR >lfine, the grid was #'',i3)') mptr1
RETURN
998 WRITE(6,'(''ERROR IN EXCHBC: THE LEVEL OF INPUT GRID'', &
& '' <1 OR >lfine, the grid was #'',i3)') mptr2
RETURN
997 WRITE(6,'(''ERROR IN EXCHBC: THE INPUT GRIDS MUST ON SAME LEVEL'' &
& ,'' the levels were '',2I4)') level1, level2
RETURN
END SUBROUTINE exchbc
!
!--------------------------------------------------------------------------
!
SUBROUTINE itpbs2(mptr1,mptr2,pts,irec,isrc,igsrc,whts, & 2,9
pts1,mxpts,nwhts,ibcvar,nbcvar,ndim)
!
! exchange boundary values for scalar varibales in IBCVAR(NBCVAR)
! between two overlapping grids at the same level by interpolation.
!
DIMENSION pts(2,mxpts),irec(2,mxpts,2),isrc(2,mxpts,2), &
igsrc(mxpts),whts(nwhts,mxpts,2),pts1(2,mxpts),ibcvar(nbcvar)
INTEGER :: npts(2)
LOGICAL :: change(2)
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'agrialloc.inc'
INCLUDE 'grddsc.inc'
nw1=nint( SQRT(FLOAT(nwhts)) )
nw2=nw1
!
! Calculate and store the coor's of p-points of the receive
! grid in array pts with their indecies stored in irec.
!
IF(nbcvar == 0) RETURN
CALL chkdim
(ndim,'ITPBS2')
ivar=ibcvar(1)
DO igrid=1,2
ig=igrid
IF( igrid == 1) THEN
mrec=mptr1
msrc=mptr2
ELSE
mrec=mptr2
msrc=mptr1
END IF
nx=node(5,mrec)
ny=node(6,mrec)
cosr=rnode(21,mrec)
sinr=rnode(22,mrec)
dxr=rnode(9 ,mrec)
dyr=rnode(10,mrec)
xor=rnode(1,mrec)
yor=rnode(2,mrec)
IF(ndim == 2) THEN
xshift=stgxy(1,ivar)
yshift=stgxy(2,ivar)
nx1=nx-idmxy(1,ivar)
ny1=ny-idmxy(2,ivar)
ELSE
xshift=stgxyz(1,ivar)
yshift=stgxyz(2,ivar)
nx1=nx-idmxyz(1,ivar)
ny1=ny-idmxyz(2,ivar)
END IF
xc=xshift*dxr*cosr
ys=yshift*dyr*sinr
xs=xshift*dxr*sinr
yc=yshift*dyr*cosr
xorp=xor+xc-ys
yorp=yor+xs+yc
!
! To calculate the absolute coordinates of the model boundary points
!
! Attention:
!
! Here we assume that in x direction the boundaries are at i=1 and nx1.
! and at j=1 and ny1 in y direction.
! If in a particular solver (model), the boundaries for certain
! variables are defined at points other the above, the use needs to
! change the following loop accordingly (i.e. change ist, jst, iend
! and jend). As a result, the part in routine UPDSCL that calculates
! the points for interior points updating also need change.
!
ist=1
jst=1
iend=nx1
jend=ny1
DO j=jst,jend,jend-jst
jj=(j-jst)/(jend-jst)*(iend-ist+1)
DO i=ist,iend
ii=jj+i-ist+1
x=(i-1)*dxr
y=(j-1)*dyr
pts (1,ii) = xorp+x*cosr-y*sinr
pts (2,ii) = yorp+x*sinr+y*cosr
irec(1,ii,ig) = i
irec(2,ii,ig) = j
END DO
END DO
numpts=MAX(0, 2*(iend-ist+1))
DO i=ist,iend,iend-ist
ii=(i-ist)/(iend-ist)*(jend-jst-1)
DO j=jst+1,jend-1
jj=ii+j-jst+numpts
x=(i-1)*dxr
y=(j-1)*dyr
pts (1,jj) = xorp+x*cosr-y*sinr
pts (2,jj) = yorp+x*sinr+y*cosr
irec(1,jj,ig) = i
irec(2,jj,ig) = j
END DO
END DO
numpts=numpts+MAX(0,2*(jend-jst-1))
!
! Find the source points for receive grid points pts
!
ierfl=0
iptcmp=1
CALL getsrc
(mrec,msrc,1,xshift,yshift,pts,numpts, &
irec(1,1,ig),igsrc,ierfl,iptcmp)
!
npts(ig)=numpts
change(ig) = .true.
IF(npts(ig) == 0) change(ig)=.false.
IF(npts(ig) == 0) CYCLE
DO ip=1,npts(ig)
pts1(1,ip)=pts(1,ip)
pts1(2,ip)=pts(2,ip)
END DO
!
! Calculate weights of interpolation for points pts. The points are
! stored in pts1 so that pts remained unchanged after calwht.
!
CALL calwht
(mrec,msrc,pts1,isrc(1,1,ig),whts(1,1,ig),npts(ig), &
nwhts,ivar,ndim)
END DO
!
! now we have receive points in irec, source points in isrc and interpolation
! weights wghts, for both interp from grid 2 to 1 and from grid 1 to 2.
! The actual interpolations are done in INTPSL for scalars listed in
! ibcvar. Interp bwteen the two grids are done at the same time
! so that packing unpacking only has to be done once.
!
nx1=node(5 ,mptr1)
ny1=node(6 ,mptr1)
nz1=node(14,mptr1)
nx2=node(5 ,mptr2)
ny2=node(6 ,mptr2)
nz2=node(14,mptr2)
IF(nz1 /= nz2)GO TO 999
!
! Interpolate for the boundary values for variables in list ibcvar
!
IF(npts(1) == 0.AND.npts(2) == 0) GO TO 500
DO i=1,nbcvar
ivar=ibcvar(i)
IF(ndim == 2) THEN
n3rd=1
iptr1=igtnxy(mptr1,ivar,1)
iptr2=igtnxy(mptr2,ivar,1)
ELSE
n3rd=nz1
iptr1=igtxyz(mptr1,ivar,1)
iptr2=igtxyz(mptr2,ivar,1)
END IF
IF(change(1)) CALL intrps
( irec(1,1,1),isrc(1,1,1),whts(1,1,1),npts(1),nw1,nw2)
IF(change(2)) CALL intrps
( irec(1,1,2),isrc(1,1,2),whts(1,1,2),npts(2),nw1,nw2)
!
! if change is false, no need to pack the data back.
!
IF(ndim == 2) THEN
CALL retnxy
(mptr1,ivar,1,iptr1,change(1))
CALL retnxy
(mptr2,ivar,1,iptr2,change(2))
ELSE
CALL retxyz
(mptr1,ivar,1,iptr1,change(1))
CALL retxyz
(mptr2,ivar,1,iptr2,change(2))
END IF
END DO
500 CONTINUE
RETURN
999 WRITE(6,'(''ERROR IN INTPS2: NZ FOR TWO GRIDS NOT EQUAL,'')')
WRITE(6,'(''THEY WERE '',2I5,'', JOB ABORTED.'')')nz1,nz2
STOP
END SUBROUTINE itpbs2
!
!--------------------------------------------------------------------------
!
SUBROUTINE itpbv2(mptr1,mptr2,pts,irec,isrc,igsrc,whts, & 2,13
pts1,mxpts,nwhts,ibcvar,nbcvar,ndim)
!
! exchange boundary values for vectors in IBCVAR(NBCVAR)
! between two overlapping grids at the same level by interpolation.
!
! input list: mptr1, mptr2,mxpts,nwhts,ibcvar,nbcvar,ndim
! working arrays: pts,irec,isrc,igsrc,whts
! index implication for e.g. irec:
!
! irec(2, mxpts, 2 )
! | | |
! x or y points components
!
! similarly for isrc, igsrc,whts.
!
DIMENSION pts(2,mxpts,2),irec(2,mxpts,2),isrc(2,mxpts,2), &
igsrc(mxpts,2),whts(nwhts,mxpts,2),pts1(2,mxpts,2), &
ibcvar(nbcvar),numpts(2),npts(2)
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'agrialloc.inc'
INCLUDE 'grddsc.inc'
LOGICAL :: samstg, smstg2, change(2)
nw1=nint( SQRT(FLOAT(nwhts)) )
nw2=nw1
!
IF(nbcvar == 0) RETURN
CALL chkdim
(ndim,'ITPBV2')
!
! To loop through all vectors, which has two components stored consecutively
! in nbcvar. The 1st compt. has to be the one in x direction.
! Here, interp weights are calculated for each vector.
!
DO ivtr=1,nbcvar,2
ivar1=ibcvar(ivtr)
ivar2=ibcvar(ivtr+1)
samstg=smstg2(ivar1,ivar2,ndim)
!
! find storage pointers for the receive grid variables (2 compnts of vector).
! uppack them if neccesary.
!
IF(ndim == 2) THEN
xshft = MAX( stgxy(1,ivar1),stgxy(1,ivar2) )
yshft = MAX( stgxy(2,ivar1),stgxy(2,ivar2) )
iptr11=igtnxy(mptr1,ivar1,1)
iptr12=igtnxy(mptr1,ivar2,1)
iptr21=igtnxy(mptr2,ivar1,1)
iptr22=igtnxy(mptr2,ivar2,1)
ELSE
xshft = MAX( stgxyz(1,ivar1),stgxyz(1,ivar2) )
yshft = MAX( stgxyz(2,ivar1),stgxyz(2,ivar2) )
iptr11=igtxyz(mptr1,ivar1,1)
iptr12=igtxyz(mptr1,ivar2,1)
iptr21=igtxyz(mptr2,ivar1,1)
iptr22=igtxyz(mptr2,ivar2,1)
END IF
DO igrid=1,2
ig=igrid
change(ig)=.true.
IF(ig == 1)THEN
mrec=mptr1
msrc=mptr2
ELSE
mrec=mptr2
msrc=mptr1
END IF
nxr=node(5,mrec)
nyr=node(6,mrec)
nzr=node(14,mrec)
cosr=rnode(21,mrec)
sinr=rnode(22,mrec)
dxr=rnode(9 ,mrec)
dyr=rnode(10,mrec)
xor=rnode(1,mrec)
yor=rnode(2,mrec)
nxs=node(5,msrc)
nys=node(6,msrc)
nzs=node(14,msrc)
coss=rnode(21,msrc)
sins=rnode(22,msrc)
cossr= coss*cosr+sins*sinr
sinsr= coss*sinr-sins*cosr
IF(nzr /= nzs) GO TO 999
n3rd=nzs
IF(ndim == 2) n3rd=1
!
! Find the points on rec. grid for both compnts. and the source grid for
! them.If the two compnts have the same staggering, they have the
! same coordinate points and the source grid, therefore those for the
! 2nd compnt. do not need to be recalculated.
!
DO icmpr=1,2
!
! icmpr=1 for the 1st component of vectors
! =2 for the 2nd component of vectors
!
IF(icmpr == 2.AND.samstg)THEN
numpts(2)=numpts(1)
DO ip=1,numpts(2)
pts (1,ip,2)=pts (1,ip,1)
pts (2,ip,2)=pts (2,ip,1)
irec(1,ip,2)=irec(1,ip,1)
irec(2,ip,2)=irec(2,ip,1)
igsrc(ip,2)=igsrc(ip,1)
END DO
CYCLE
END IF
ivar=ibcvar(ivtr+icmpr-1)
IF(ndim == 2) THEN
xshift=stgxy(1,ivar)
yshift=stgxy(2,ivar)
nx1=nxr-idmxy(1,ivar)
ny1=nyr-idmxy(2,ivar)
ELSE IF(ndim == 3)THEN
xshift=stgxyz(1,ivar)
yshift=stgxyz(2,ivar)
nx1=nxr-idmxyz(1,ivar)
ny1=nyr-idmxyz(2,ivar)
END IF
xc=xshift*dxr*cosr
ys=yshift*dyr*sinr
xs=xshift*dxr*sinr
yc=yshift*dyr*cosr
xorp=xor+xc-ys
yorp=yor+xs+yc
!
! To calculate the absolute coordinates of the model boundary points
!
! Attention:
!
! Here we assume that in x direction the boundaries are at i=1 and nx1.
! and at j=1 and ny1 in y direction.
! If in a particular solver (model), the boundaries for certain
! variables are defined at points other the above, the use needs to
! change the following loop accordingly (i.e. change ist, jst, iend
! and jend). As a result, the part in routine UPDVTR that calculates
! the points for interior points updating also need change.
!
ist=1
jst=1
iend=nx1
jend=ny1
DO j=jst,jend,jend-jst
jj=(j-jst)/(jend-jst)*(iend-ist+1)
DO i=ist,iend
ii=jj+i-ist+1
x=(i-1)*dxr
y=(j-1)*dyr
pts (1,ii,icmpr) = xorp+x*cosr-y*sinr
pts (2,ii,icmpr) = yorp+x*sinr+y*cosr
irec(1,ii,icmpr) = i
irec(2,ii,icmpr) = j
END DO
END DO
numpts(icmpr)=MAX(0,2*(iend-ist+1))
DO i=ist,iend,iend-ist
ii=(i-ist)/(iend-ist)*(jend-jst-1)
DO j=jst+1,jend-1
jj=ii+j-jst+numpts(icmpr)
x=(i-1)*dxr
y=(j-1)*dyr
pts (1,jj,icmpr) = xorp+x*cosr-y*sinr
pts (2,jj,icmpr) = yorp+x*sinr+y*cosr
irec(1,jj,icmpr) = i
irec(2,jj,icmpr) = j
END DO
END DO
numpts(icmpr)=numpts(icmpr)+MAX(0,2*(jend-jst-1))
!
! Find the source grid for receive grid points pts
!
ierfl=0
iptcmp=1
CALL getsrc
(mrec,msrc,1,xshft,yshft,pts(1,1,icmpr), &
numpts(icmpr),irec(1,1,icmpr),igsrc(1,icmpr),ierfl,iptcmp)
npts(icmpr)=numpts(icmpr)
IF(npts(icmpr) == 0) change(ig)=.false.
!
END DO
DO icmps=1,2
DO icmpr=1,2
IF(icmpr == 2.AND.samstg)THEN
! if the two components have the same staggering, simply make a copy for the
! 2nd component.
npts(2)=npts(1)
IF(npts(2) == 0) CYCLE
DO ip=1,npts(2)
pts (1,ip,2)=pts (1,ip,1)
pts (2,ip,2)=pts (2,ip,1)
isrc(1,ip,2)=isrc(1,ip,1)
isrc(2,ip,2)=isrc(2,ip,1)
irec(1,ip,2)=irec(1,ip,1)
irec(2,ip,2)=irec(2,ip,1)
END DO
DO iw=1,nwhts
DO ip=1,npts(2)
whts(iw,ip,2)=whts(iw,ip,1)
END DO
END DO
CYCLE
END IF
IF(npts(icmpr) == 0) CYCLE
DO ip=1,npts(icmpr)
pts1(1,ip,icmpr)=pts(1,ip,icmpr)
pts1(2,ip,icmpr)=pts(2,ip,icmpr)
END DO
!
! Calculate weights of interpolation for points pts. The points are
! stored in pts1 so that pts remained unchanged after calwht.
!
ivsrc=ibcvar(ivtr+icmps-1)
CALL calwht
(mrec,msrc,pts1(1,1,icmpr),isrc(1,1,icmpr), &
whts(1,1,icmpr),npts(icmpr),nwhts,ivsrc,ndim)
END DO
IF(npts(1) == 0.AND.npts(2) == 0) CYCLE
!
! find pointer for one of the components (icmps) of the source grid
! vector. If it is packed, do uppacking for this variable.
!
ivsrc=ibcvar(ivtr+icmps-1)
IF(icmps == 1)THEN
proj1=cossr
proj2=-sinsr
ELSE
proj1=sinsr
proj2=cossr
END IF
IF(ig == 1)THEN
irptr1=iptr11
irptr2=iptr12
isptr1=iptr21
isptr2=iptr22
ELSE
irptr1=iptr21
irptr2=iptr22
isptr1=iptr11
isptr2=iptr12
END IF
IF(icmps == 1) THEN
isptr=isptr1
ELSE
isptr=isptr2
END IF
icmpr=1
IF(npts(icmpr) /= 0) THEN
CALL intrpv
( irec(1,1,icmpr),isrc(1,1,icmpr),whts(1,1,icmpr), &
npts(icmpr),nw1,nw2,icmps,proj1)
END IF
icmpr=2
IF(npts(icmpr) /= 0) THEN
CALL intrpv
( irec(1,1,icmpr),isrc(1,1,icmpr),whts(1,1,icmpr), &
npts(icmpr),nw1,nw2,icmps,proj2)
END IF
END DO
END DO
IF(ndim == 2) THEN
CALL retnxy
(mptr1,ivar1,1,iptr11,change(1))
CALL retnxy
(mptr1,ivar2,1,iptr12,change(1))
CALL retnxy
(mptr2,ivar1,1,iptr21,change(2))
CALL retnxy
(mptr2,ivar2,1,iptr22,change(2))
ELSE
CALL retxyz
(mptr1,ivar1,1,iptr11,change(1))
CALL retxyz
(mptr1,ivar2,1,iptr12,change(1))
CALL retxyz
(mptr2,ivar1,1,iptr21,change(2))
CALL retxyz
(mptr2,ivar2,1,iptr22,change(2))
END IF
END DO
RETURN
999 WRITE(6,'(''ERROR IN INTPV2: NZ FOR TWO GRIDS NOT EQUAL,'')')
WRITE(6,'(''THEY WERE '',2I5,'', JOB ABORTED.'')')nzr,nzs
STOP
END SUBROUTINE itpbv2