---

Absorbets

      subroutine absorber(z,x,y,tau)

      dimension tau(l,1-ih:np+ih,1-ih:mp+ih,2)

      common/davies/ relx(np,mp),rely(np,mp),zab,towxL,towxR,towy,towz,
     1               dxabL,dxabR,dyab,iab,iabth,iabqw

compute absorbers at x boundaries
      if ( irlx.eq.1 ) then
      toliL=1./towxL
      toliR=1./towxR
      do j=1,mp
       xab1=xr1(j)*xc1(j)+dxabL  !VEC xab1=xcr(1,j)*cosa(1,j)+dxabL
       xab2=xr2(j)*xc2(j)-dxabR  !VEC xab2=xcr(n,j)*cosa(n,j)-dxabR
       do i=1,np
        ia=(npos-1)*np + i
        r0=xcr(i,j)*cosa(i,j)
        r1=amax1(0., xab1-r0)/dxabL
        r2=amax1(0., r0-xab2)/dxabR
        relb=toliL*r1+toliR*r2-2.*toliL*toliR*r1*r2/(toliL+toliR) 
         tmp(i,j)=irlx*relb
         relx(i,j)=tmp(i,j)

compute absorbers at y boundaries
      if (j3.eq.1) then
      if (irly.eq.1) then
      toli=1./towy
      do j=1,mp
        do i=1,np
         yab1=yr1(i)+dyab         !VEC yab1=ycr(i,1 )+dyab
         yab2=yr2(i)-dyab         !VEC yab2=ycr(i,mp)-dyab
         r0=ycr(i,j)
         r1=amax1(0.,-r0+yab1)/dyab
         r2=amax1(0., r0-yab2)/dyab
         relb=toli*(r1+r2-r1*r2)
         tmp(i,j)=irly*relb
         rely(i,j)=tmp(i,j)

c     --- compute absorbers at upper boundary
      towi=1./towz                             
      abscl=5.7e3        ! for negative zab
      do k=1,l
        do j=1,mp
          do i=1,np
            zl=zstr(k)/gi(i,j)+zs(i,j)
            if(zab.ge.0.) then
             t1=iab*amax1(0.,zl-zab)
             tau(k,i,j,1)=towi*t1/(zb-zab)
             t2=iabth*amax1(0.,zl-zab)
             tau(k,i,j,2)=towi*t2/(zb-zab)
            else
             tau(k,i,j,1)=iab*towi*exp((zl-zb)/abscl)
             tau(k,i,j,2)=iabth*towi*exp((zl-zb)/abscl)
            endif

BC for advection

      subroutine mp3bc(x,iflg,bcx,bcy,n1,n2,n3)

      dimension x(1-ih:np+ih,1-ih:mp+ih,l),
     .        bcx(1-ih:mp+ih,n3,2),
     .        bcy(1-ih:np+ih,n3,2),

      common/profl/ th0(1-ih:np+ih,1-ih:mp+ih,l),
     .              rho(1-ih:np+ih,1-ih:mp+ih,l),
     .              prf(1-ih:np+ih,1-ih:mp+ih,l,3),
     .              zcr(l)

      if(iflg.eq.1) then                   ! ---> th
      if(implgw.eq.1) then
         if (leftedge.eq.1) then
            do k=1,n3
               do j=1,n2
                bcx(j,k,1)=0.
                bcx(j,k,2)=0.
               enddo
            enddo
         else if (rightedge.eq.1) then

 ....................

      if(iflg.eq.2.or.iflg.eq.3) then      ! ---> u0.v0
         if (leftedge.eq.1) then
            do k=1,n3
               do j=1,n2
                bcx(j,k,1)=prf( 1,j,k,iflg)
                bcx(j,k,2)=prf( 0,j,k,iflg)

......................

#if (MOISTMOD > 0)
      if(iflg.eq.6) then                   ! ---> qv
         if (leftedge.eq.1) then
            do k=1,n3
               do j=1,n2
                bcx(j,k,1)=qve( 1,j,k)
                bcx(j,k,2)=qve( 0,j,k)

BC for pressure solver

      subroutine velbc(ue,ve,rho)

C----------------------------------------------------------------------------
C   This routine sets the boundary conditions for normal solenoidal
C   velocities at the model boundaries. There is no universal way to 
C   do this. There are two important formulae to consider:  
C      a) V_solenoidal=V_contravariant- partial X_bar/partial t
C      b) V_solenoidal= [Metric-Coefficient Matrix]*V_physical. 
C   Both are equivalent but their relative utility depends on the 
C   problem at hand.
C
C   In the code below vertical-boundary solenoidal velocity is set assuming
C   that the curvilinear boundary is impermeable, whereupon contravariant 
C   velocity vanishes there (fluid adheres to the boundary), so relation (a)
C   is employed. Horizontal components exploit (b) formulae while assuming
C   known physical velocities at the boundary (here profiles ue and ve,
C   but anything desired can be passed by the subroutine call). Alternate
C   formulations are included (but commented out) for the use with special
C   applications (e.g., flow past chanels with vibrating lateral walls)
C----------------------------------------------------------------------------

      do 3200 k=1,l,l-1
      kk=1+k/l
      do 320 j=1,mp
      do 320 i=1,np
         g33=gi(i,j)*gmus(k)
      ob(i,j,kk)=gmul(k)*zsd(i,j)/zb*g33-(gmul(k)/zb-1.)*zhd(i,j)*g33  !(a)
C        g110=1./((1-icylind)*gmm(i,j,k)*cosa(i,j)+icylind*1.)
C        g220=1./gmm(i,j,k)
C        g13=(s13(i,j)*gmul(k)-h13(i,j))*gmus(k)*g110
C        g23=(s23(i,j)*gmul(k)-h23(i,j))*gmus(k)*g220
C     ob(i,j,kk)=g13*ue(i,j,k)+g23*ve(i,j,k)+g33*we        !alternative (b)
  320 continue
 3200 continue

      if (leftedge.eq.1) then
      do 321 k=1,l
      do 321 j=1,mp
         ja=(mpos-1)*mp + j
          g110w=1./((1-icylind)*gmm(1,j,k)*cosa(1,j)+icylind*1.)
          g220w=1./gmm(1,j,k)
          g11w=strxx(1,j)*g110w
          g21w=strxy(1,j)*g220w
      ub(j,k,1)=g11w*ue(1,j,k)+g21w*ve(1,j,k)
C     ub(j,k,1)=-strxd(1,j)                              !alternative (a)
#if (PARALLEL > 0)
          g110w=1./((1-icylind)*gmm(0,j,k)*cosa(0,j)+icylind*1.)
          g220w=1./gmm(0,j,k)
          g11w=strxx(0,j)*g110w
          g21w=strxy(0,j)*g220w
      ub(j,k,2)=g11w*ue(0,j,k)+g21w*ve(0,j,k)
C     ub(j,k,2)=-strxd(0,j)                              !alternative (a)
#endif

      .........

       vb(i,k,2)=g12n*ue(i,mp,k)+g22n*ve(i,mp,k)
C      vb(i,k,2)=-stryd(i,mp)                            !alternative (a)

      ........

      call vbcad(rho)

      return
      end

Corrections for outflow velocities in open BC mode

      subroutine vbcad(d)

      if((ibcx.eq.1).and.(ibcy.eq.1.or.j3.eq.0)) return

character of the adjustement of velocities at the boundaries:
c iflg=0 multiplicative adjustment of outflow velocities only;
c iflg.ne.0 additive adjustement of both inflow and outflow velocities
c z boundary has been coded in for special purpose exepriments but
commented out as omega=0 is required at z=0 and h
      iflg=0

customarily there are no fluxes through z boundaries (except for that
caused by time dependend orography); typically there is no need for a
code computing adjustemnt at z boundaries; however, in special cases such
code may be helpful. izflg=0 baypasses the special code, izlg=1 activates it.
      izflg=1
      izflg=0   ! default

compute adjustement of outflow velocities at the boundaries:

cofluxes through x boundaries

         do k=1,l
            do j=1,mp
             tmp1(1,j,k)=wgz(k)*wgy(j)*
     .         (d(np,j,k)*pp(ub(j,k,2))-d(np+1,j,k)*pn(ub(j,k,1)))
             tmp2(2,j,k)=wgz(k)*wgy(j)*
     .         (d(np,j,k)*pn(ub(j,k,2))-d(np+1,j,k)*pp(ub(j,k,1)))
             tmp3(3,j,k)=wgz(k)*wgy(j)*(d(np,j,k)+d(np+1,j,k))
            end do
         end do
      end if

      udb = udb*dy*dz
      uout=uout*dy*dz
      uinf=uinf*dy*dz

cofluxes through y boundaries
         do k=1,l
            do i=1,np
               tmp1(i,1,k)=wgz(k)*wgx(i)*
     .           (d(i,mp,k)*pp(vb(i,k,2))-d(i,mp+1,k)*pn(vb(i,k,1)))
               tmp2(i,2,k)=wgz(k)*wgx(i)*
     .           (d(i,mp,k)*pn(vb(i,k,2))-d(i,mp+1,k)*pp(vb(i,k,1)))
               tmp3(i,3,k)=wgz(k)*wgx(i)*(d(i,mp,k) + d(i,mp+1,k))
            end do
         end do

      vdb = vdb*dx*dz
      vout=vout*dx*dz
      vinf=vinf*dx*dz
   20 continue

cofluxes through z boundaries
      if(izflg.eq.1) then

      do j=1,mp
         do i=1,np
            tmp1(i,j,1)=wgx(i)*wgy(j)*
     .           (d(i,j,l)*pp(ob(i,j,2))-d(i,j,1)*pn(ob(i,j,1)))
            tmp2(i,j,2)=wgx(i)*wgy(j)*
     .           (d(i,j,l)*pn(ob(i,j,2))-d(i,j,1)*pp(ob(i,j,1)))
            tmp3(i,j,3) = wgx(i)*wgy(j)*(d(i,j,l)+d(i,j,1))
         end do
      end do

      do j=1,mp
         do i=1,np
           zbrdt= gmul(1)/zb*gi(i,j)*gmus(1)*zsd(i,j)
     .      -(gmul(1)/zb-1.)*zhd(i,j)*gi(i,j)*gmus(1)
           zhrdt=-(gmul(l)/zb-1.)*zhd(i,j)*gi(i,j)*gmus(l)
     .      +gmul(l)/zb*gi(i,j)*gmus(l)*zsd(i,j)
           tmp1(i,j,1) = wgx(i)*wgy(j)*(d(i,j,l)*zhrdt-d(i,j,1)*zbrdt)
         end do
      end do

      odb = odb*dx*dy
      oout=oout*dx*dy
      oinf=oinf*dx*dy
      tflx=tflx*dx*dy

constant of mass adjustement
      if(uout+vout+oout.eq.0.) iflg=1
      epsia=0.
      if(iflg.eq.0) then
      epsim=-(uinf+vinf+oinf+tflx)/(uout+vout+oout)
      else
      epsim=1.
      epsia=-(uinf+vinf+oinf+tflx+uout+vout+oout)/(udb+vdb+odb)
      endif

corrected outflow velocities
      if(ibcx.eq.1) goto 111
      if (leftedge.eq.1) then
      do 11 k=1,l
      do 11 j=1,mp
   11 ub(j,k,1)=pp(ub(j,k,1))+epsim*pn(ub(j,k,1))-epsia
      endif
      if(rightedge.eq.1) then
      do 12 k=1,l
      do 12 j=1,mp
   12 ub(j,k,2)=pn(ub(j,k,2))+epsim*pp(ub(j,k,2))+epsia
      endif
  111 continue

      if(ibcy.eq.1.or.j3.eq.0) goto 222
      if (botedge.eq.1) then
      do 22 k=1,l
      do 22 i=1,np
   22 vb(i,k,1)=pp(vb(i,k,1))+epsim*pn(vb(i,k,1))-epsia
      endif
      if (topedge.eq.1) then
      do 23 k=1,l
      do 23 i=1,np
   23 vb(i,k,2)=pn(vb(i,k,2))+epsim*pp(vb(i,k,2))+epsia
      endif
  222 continue

      do 33 j=1,mp
      do 33 i=1,np
      ob(i,j,2)=pn(ob(i,j,2))+epsim*pp(ob(i,j,2))+epsia
   33 ob(i,j,1)=pp(ob(i,j,1))+epsim*pn(ob(i,j,1))-epsia