!
!--------------------------------------------------------------------------
!
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