copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
***s/r nesajr -- Horizontal blending between internal model solution(fn)
*                and external boundary conditions(bcs,bcn,bcw.bce).
*

      subroutine nesajr  (fn,bcs,bcn,bcw,bce,lminx,lmaxx,lminy,lmaxy, 14
     $       lminxs,lmaxxs,lminys,lmaxys,lminxw,lmaxxw,lminyw,lmaxyw,
     $                                                     lnk,is,js)
      implicit none
*
      integer lminx,lmaxx,lminy,lmaxy,lminxs,lmaxxs,lminys,lmaxys,
     $        lminxw,lmaxxw,lminyw,lmaxyw,lnk,is,js
      real fn(lminx:lmaxx,lminy:lmaxy,lnk),
     $     bcs(lminxs:lmaxxs,lminys:lmaxys,*),
     $     bcn(lminxs:lmaxxs,lminys:lmaxys,*),
     $     bcw(lminxw:lmaxxw,lminyw:lmaxyw,*),
     $     bce(lminxw:lmaxxw,lminyw:lmaxyw,*)
*
*OBJECT
*     This subroutine blends the field "fn", updated by a dynamic
*     timestep, with the nesting field "bc*". Nesting is performed
*     according to equation (4.1.1) on hblen_x points along the x-axis 
*     and on hblen_y points along the y-axis. The attenuation function
*     "p" is as described by equations 4.1.2, 4.1.3 abd 4.1.4.
*
*ARGUMENTS
*    NAMES     I/O  TYPE  A/S        DESCRIPTION
*
*    fn         O     R    A    field to be nested
*    bc*        I     R    A    nesting fields
*    is         I     I    S    staggering parameter along x
*    js         I     I    S    staggering parameter along y
*
*IMPLICIT
#include "lcldim.cdk"
#include "bcsdim.cdk"
#include "nestpnt.cdk"
*
**
      integer i,j,k,nit,njt,il,ih,jl,jh
      real*8 p,one
      parameter (one = 1.0d0)
*----------------------------------------------------------------------
*
      if (south+north+west+east.eq.0) return
*
      nit  = ldni-is* east
      njt  = ldnj-js*north
*
      il   = 1   + hblen_x* west
      ih   = nit - hblen_x* east
      jl   = 1   + hblen_y*south
      jh   = njt - hblen_y*north
*
      if (west.eq.1) then
         do k=1,lnk
         do j=jl,jh
         do i=1,il-1
	    p = wh_w(i)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcw(i,j,k)
         end do
         end do
         end do
      endif
      if (east.eq.1) then
         do k=1,lnk
         do j=jl,jh
         do i=ih+1,nit
            p = wh_e(i+is)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bce(i-bcs_ofi,j,k)
         end do
         end do
         do j=minye,maxye
         do i=nit+1,ldni
            fn(i,j,k) = bce(i-bcs_ofi,j,k)
         end do
         end do
         end do
      endif
      if (south.eq.1) then
         do k=1,lnk
         do j=1,jl-1
         do i=il,ih
            p = wh_s(j)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcs(i,j,k)
         end do
         end do
         end do
      endif
      if (north.eq.1) then
         do k=1,lnk
         do j=jh+1,njt
         do i=il,ih
            p = wh_n(j+js)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcn(i,j-bcs_ofj,k)
         end do
         end do
         do j=njt+1,ldnj
         do i=1,ldni
            fn(i,j,k) = bcn(i,j-bcs_ofj,k)
         end do
         end do
         end do
      endif
*
      if (south*west.eq.1) then
         do k=1,lnk
         do j=1,jl-1
         do i=1,il-1
            p = wh_sw(i,j)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcs(i,j,k)
         end do
         end do
         end do
      endif
      if (south*east.eq.1) then
         do k=1,lnk
         do j=1,jl-1
         do i=ih+1,nit
            p = wh_se(i+is,j)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcs(i,j,k)
         end do
         end do
         do i=nit+1,ldni
         do j=minys,maxys
            fn(i,j,k) = bcs(i,j,k)
         end do
         end do
         end do
      endif
      if (north*west.eq.1) then
         do k=1,lnk
         do i=1,il-1
         do j=jh+1,njt
            p = wh_nw(i,j+js)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcn(i,j-bcs_ofj,k)
         end do
         end do
         end do
      endif
      if (north*east.eq.1) then
         do k=1,lnk
         do i=ih+1,nit
         do j=jh+1,njt
             p = wh_ne(i+is,j+js)
            fn(i,j,k)= (one-p)*fn(i,j,k)+p*bcn(i,j-bcs_ofj,k)
         end do
         end do
         do i=nit+1,ldni
         do j=minyn,maxyn
            fn(i,j,k) = bcn(i,j-bcs_ofj,k)
         end do
         end do
         end do
      endif
*
*----------------------------------------------------------------------
      return
      end