!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE UPDEXBC                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE updexbc(mptr,iexbcbuf) 1,12
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate the external boundary data from base grid to grid mptr
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  04/15/1997
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!  mptr       Grid ID number
!
!  iexbcbuf   Pointer that points to the location of external
!             boundary data for base grid.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'nodal.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'agrialloc.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'exbc.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: lisvar(22)
  INTEGER :: mptr, iexbcbuf

  INTEGER :: nxs,nys,nzs, nxys, nxyzs
  INTEGER :: nx,ny,nz, nxy,nxyz
  INTEGER :: i,j,k,ijk

  INTEGER :: iu0exbc,iv0exbc,iw0exbc,ipt0exbc
  INTEGER :: iudtexbc,ivdtexbc,iwdtexbc,iptdtexbc
  INTEGER :: iuexbc,ivexbc,iwexbc,iptexbc

  INTEGER :: iptr1,iptr2,iptr3,iptr4,iptr5,iptr6

  INTEGER :: item3d,nwhts

  REAL :: delta_t
!
!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------
!
  INTEGER :: igetsp,igtexbc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF( intrpodr == 2 ) THEN
    nwhts=9
  ELSE
    nwhts=4
  END IF

  nxs=node(5,1)
  nys=node(6,1)
  nzs=node(14,1)

  nx=node(5,mptr)
  ny=node(6,mptr)
  nz=node(14,mptr)

  nxy =nx*ny
  nxyz=nxy*nz

  nxys =nxs*nys
  nxyzs=nxys*nzs

  iptr1=igetsp(2*nxy)
  iptr2=igetsp(2*nxy)
  iptr3=igetsp(2*nxy)
  iptr4=igetsp(nxy)
  iptr5=igetsp(nwhts*nxy)
  iptr6=igetsp(2*nxy)

  item3d = igetsp( nxyzs )

  delta_t = curtim - ( abstfcst0 - abstinit )
!
!-----------------------------------------------------------------------
!
!  Interpolate uexbc
!
!-----------------------------------------------------------------------
!
  iu0exbc  = igtexbc(1,1,1)
  iudtexbc = igtexbc(1,12,1)

  DO k=1,nzs-1
    DO j=1,nys-1
      DO i=1,nxs
        ijk = (k-1)*nxs*nys + (j-1)*nxs + i
        a(item3d+ijk-1) = a(iu0exbc+ijk-1) + a(iudtexbc+ijk-1)*delta_t
      END DO
    END DO
  END DO

  iuexbc = iexbcbuf
  lisvar(1) = 1
  CALL itpnexbc(mptr,lisvar,1,item3d,iuexbc,                            &
       a(iptr1),a(iptr2),a(iptr3),                                      &
       a(iptr4),a(iptr5),a(iptr6),nxy,nwhts)
!
!-----------------------------------------------------------------------
!
!  Interpolate vexbc
!
!-----------------------------------------------------------------------
!
  iv0exbc  = igtexbc(1,2,1)
  ivdtexbc = igtexbc(1,13,1)

  DO k=1,nzs-1
    DO j=1,nys
      DO i=1,nxs-1
        ijk = (k-1)*nxs*nys + (j-1)*nxs + i
        a(item3d+ijk-1) = a(iv0exbc+ijk-1) + a(ivdtexbc+ijk-1)*delta_t
      END DO
    END DO
  END DO

  ivexbc = iexbcbuf + 1*nxyz
  lisvar(1) = 2
  CALL itpnexbc(mptr,lisvar,1,item3d,ivexbc,                            &
       a(iptr1),a(iptr2),a(iptr3),                                      &
       a(iptr4),a(iptr5),a(iptr6),nxy,nwhts)
!
!-----------------------------------------------------------------------
!
!  Interpolate wexbc
!
!-----------------------------------------------------------------------
!
  iw0exbc = igtexbc(1,3,1)
  iwdtexbc = igtexbc(1,14,1)

  DO k=1,nzs
    DO j=1,nys-1
      DO i=1,nxs-1
        ijk = (k-1)*nxs*nys + (j-1)*nxs + i
        a(item3d+ijk-1) = a(iw0exbc+ijk-1) + a(iwdtexbc+ijk-1)*delta_t
      END DO
    END DO
  END DO

  iwexbc = iexbcbuf + 2*nxyz
  lisvar(1) = 3
  CALL itpnexbc(mptr,lisvar,1,item3d,iwexbc,                            &
       a(iptr1),a(iptr2),a(iptr3),                                      &
       a(iptr4),a(iptr5),a(iptr6),nxy,nwhts)
!
!-----------------------------------------------------------------------
!
!  Interpolate scalar ptexbc.
!
!-----------------------------------------------------------------------
!
  ipt0exbc  = igtexbc(1,4,1)
  iptdtexbc = igtexbc(1,15,1)

  DO k=1,nzs-1
    DO j=1,nys-1
      DO i=1,nxs-1
        ijk = (k-1)*nxs*nys + (j-1)*nxs + i
        a(item3d+ijk-1) = a(ipt0exbc+ijk-1)+a(iptdtexbc+ijk-1)*delta_t
      END DO
    END DO
  END DO

  iptexbc = iexbcbuf + 3*nxyz
  lisvar(1) = 4

  CALL itpnexbc(mptr,lisvar,1,item3d,iptexbc,                           &
       a(iptr1),a(iptr2),a(iptr3),                                      &
       a(iptr4),a(iptr5),a(iptr6),nxy,nwhts)

  CALL retexbc(1,1,1,iu0exbc,   .true.)
  CALL retexbc(1,2,1,iv0exbc,   .true.)
  CALL retexbc(1,3,1,iw0exbc,   .true.)
  CALL retexbc(1,4,1,ipt0exbc,  .true.)
  CALL retexbc(1,12,1,iudtexbc, .true.)
  CALL retexbc(1,13,1,ivdtexbc, .true.)
  CALL retexbc(1,14,1,iwdtexbc, .true.)
  CALL retexbc(1,15,1,iptdtexbc,.true.)

  RETURN
END SUBROUTINE updexbc
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE ITPNEXBC                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE itpnexbc(mrec,lisvar,nvar,isptr,irptr,                       & 4,3
           pts,irec,isrc,igsrc,whts,pts1,nptrec,nwhts)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate the external boundary data from base grid to grid mrec
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  04/15/1997
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: mrec,nptrec,nwhts
  INTEGER :: isptr,irptr, nvar
  INTEGER :: irec(2,nptrec),isrc(2,nptrec),igsrc(nptrec)
  INTEGER :: lisvar(nvar)

  REAL :: pts(2,nptrec),whts(nwhts,nptrec),pts1(2,nptrec)

  INTEGER :: nxr,nyr,nzr
  INTEGER :: nxs,nys,nzs
  INTEGER :: nx1,ny1
  INTEGER :: ist,jst,iend,jend

  INTEGER :: ivar
  INTEGER :: nw1,nw2
  INTEGER :: numpts,npts0
  INTEGER :: n3rd
  INTEGER :: ierfl,iptcmp

  INTEGER :: i,ii,j,jj, ip

  REAL :: cosr,sinr
  REAL :: x,y,xc,yc,xs,ys
  REAL :: xor,yor,xorp,yorp
  REAL :: dxr,dyr
  REAL :: xshift,yshift
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'nodal.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'agrialloc.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nw1=nint( SQRT(FLOAT(nwhts)) )
  nw2=nw1

  IF(nvar == 0) RETURN

  nxr=node(5,mrec)
  nyr=node(6,mrec)
  nzr=node(14,mrec)

  nxs=node(5,1)
  nys=node(6,1)
  nzs=node(14,1)

  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=lisvar(1)
  xshift=stgexbc(1,ivar)
  yshift=stgexbc(2,ivar)
  nx1=nxr-idmexbc(1,ivar)
  ny1=nyr-idmexbc(2,ivar)

  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 all grid points including
! boundary and interior.
!
  ist=1
  jst=1
  iend=nx1
  jend=ny1
  DO j=jst,jend
    jj=(j-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) = i
      irec(2,ii) = j
    END DO
  END DO
  numpts=(iend-ist+1)*(jend-jst+1)
  npts0=numpts
!
! Find the source points for receive grid points pts
!
  ierfl=0
  iptcmp=1

  CALL getsrc(mrec,1, 1 ,xshift,yshift,pts,numpts,                      &
              irec,igsrc,ierfl,iptcmp)

  IF( (numpts /= npts0) .AND. verbose6 ) THEN
    WRITE(6,'('' WARNING: IN INTNSL,'',I4,'' POINTS ON GRID '',I4,      &
    &   '' got no source from source grid '',i4)')                      &
        npts0-numpts, mrec, 1
    WRITE(6,'('' IT WAS FOR EXBC VARIABLE NO. '',I4)') ivar
  END IF
!
! Calculate weights of interpolation for points pts.
!
  DO ip=1,numpts
    pts1(1,ip)=pts(1,ip)
    pts1(2,ip)=pts(2,ip)
  END DO

  CALL calwht_exbc(mrec,pts1,isrc,whts,numpts,nwhts,ivar)
!
! Interpolate for variables in list lisvar from source grid MSRC
! to receive grid MREC.
!
  DO i=1,nvar
    ivar=lisvar(i)
    n3rd=nzs

    CALL intrps(        irec,isrc,whts,numpts,nw1,nw2)

  END DO

  500   CONTINUE

  RETURN
END SUBROUTINE itpnexbc
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE CALWHT_EXBC              ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE calwht_exbc(mrec,points,isrc,whts,numpts,nwhts,ivar) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  this routine calculates the weights for the interpolation
!  of mass points or vertical velocity (stuff that doesn't need
!  any rotation from base grid to grid mrec
!
!  if nw = 4 bilinear interp., = 9 bi-quadratic
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  04/15/1997
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!  mrec       Grid ID number
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: mrec
  INTEGER :: numpts,nwhts

  INTEGER :: isrc(2,numpts)
  REAL :: points(2,numpts),whts(nwhts,numpts)

  REAL :: xshift,yshift, yp
  REAL :: alpha,alpha2

  REAL :: tempx1,tempx2,tempx3
  REAL :: tempy1,tempy2,tempy3

  INTEGER :: i,ivar
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'nodal.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'agrialloc.inc'
  INCLUDE 'grddsc.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  REAL :: cosc,sinc
  REAL :: dxc,dyc
  REAL :: odxc,odyc
  REAL :: xor,yor,xorp,yorp
  REAL :: xc,xs,yc,ys,x,y
!
!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------
!
  REAL :: fq, fl
  REAL :: xp,x1,x2,x3,x4

  fq(xp,x1,x2,x3,x4)=(xp-x2)*(xp-x3)/((x1-x2)*(x1-x3))+x4
  fl(xp,x1,x2)=(xp-x2)/(x1-x2)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  cosc = rnode( 21,1)
  sinc = rnode( 22,1)
  dxc  = rnode(9,1)
  dyc  = rnode(10,1)
  odxc = 1./dxc
  odyc = 1./dyc
  xor  = rnode(1,1)
  yor  = rnode(2,1)

  xshift=stgexbc(1,ivar)
  yshift=stgexbc(2,ivar)

  xc=xshift*dxc*cosc
  ys=yshift*dyc*sinc
  xs=xshift*dxc*sinc
  yc=yshift*dyc*cosc

  xorp=xor+xc-ys
  yorp=yor+xs+yc

  DO i=1,numpts
    x=points(1,i)
    y=points(2,i)
    points(1,i)=1+( cosc*(x-xorp)+sinc*(y-yorp))*odxc
    points(2,i)=1+(-sinc*(x-xorp)+cosc*(y-yorp))*odyc
  END DO

!
!   next, calculate weights for the interpolation
!
  IF ( nwhts == 9 ) THEN
    alpha=0.
    alpha2=0.
    DO i=1,numpts
      xp=points(1,i)-FLOAT( nint(points(1,i)) )
      yp=points(2,i)-FLOAT( nint(points(2,i)) )
      tempy1=fq(yp,-1.,0.,1.,alpha)
      tempy2=fq(yp,0.,-1.,1.,alpha2)
      tempy3=fq(yp,1.,-1.,0.,alpha)
      tempx1=fq(xp,-1.,0.,1.,alpha)
      tempx2=fq(xp,0.,-1.,1.,alpha2)
      tempx3=fq(xp,1.,-1.,0.,alpha)
      whts(1,i)=tempy1*tempx1
      whts(2,i)=tempy1*tempx2
      whts(3,i)=tempy1*tempx3
      whts(4,i)=tempy2*tempx1
      whts(5,i)=tempy2*tempx2
      whts(6,i)=tempy2*tempx3
      whts(7,i)=tempy3*tempx1
      whts(8,i)=tempy3*tempx2
      whts(9,i)=tempy3*tempx3
      isrc(1,i)=nint(points(1,i))
      isrc(2,i)=nint(points(2,i))
    END DO
  ELSE
    DO i=1,numpts
      xp=points(1,i)-FLOAT( INT(points(1,i)) )
      yp=points(2,i)-FLOAT( INT(points(2,i)) )
      tempy1=fl(yp,0.,1.)
      tempy2=fl(yp,1.,0.)
      tempx1=fl(xp,0.,1.)
      tempx2=fl(xp,1.,0.)
      whts(1,i)=tempy1*tempx1
      whts(2,i)=tempy1*tempx2
      whts(3,i)=tempy2*tempx1
      whts(4,i)=tempy2*tempx2
      isrc(1,i)=INT(points(1,i))
      isrc(2,i)=INT(points(2,i))
    END DO
  END IF

  RETURN
END SUBROUTINE calwht_exbc