c
c     @(#)proect.f	1.1
c
      SUBROUTINE PROECT(RS,K,IVT)
c
C----------------------------------------------------------------------
C
C  THIS SUBROUTINE COMPUTES BOUNDARY CONDITIONS FOR VORTICITY.
C
C
C  ARGUMENTS:
C            L=140km -  corr radius
c            a=2.1038*dhor/140.0          
c            r=x-ct  - coordinate
C        k - level   
c        ivt = 1/0  -switcher between vorticity/topdensity
c
c     INCLUDE 'qg_par.inc'
      INCLUDE 'mladv.inc'
C

      real chh,cxx,cyy,czz,a,ar,ar3
      dimension rs(mn0)
      COMMON/PARM/KZ,M,MM1,MM2,MP1,MSQ,TIM,TMAX,TSTART,HY,REA,DT,
     +   XBASIN,YY(MX0),ICALC,KZP1,KZM1,IFTOP,IFBOT,IFTVV,IFREL,THETA,
     +   ifsbl,ifeva
      COMMON/THREED/HZ(K0),SCSIG,SIGZ(K01),EIGVAL(K0),AINV(K0,K0),
     +      AMAT(K0,K0)
      COMMON/EXTPAR/RLAT0,RLNG0,T0,V0,DHOR,HT,TIME0,
     +                R0,F0STAR
      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)
      dimension hz1(k0)
C
C --- FOR BOUNDARY COMPUTATION
c
c --- compute boundary condition for vorticity as sum of second 
C --- derivative of streamfunction along boundary, extention 
C --- streamfunction across boundary using corr function and thermal vorticity
c
c ---  computute corr function for vorticity
c
      DATA inach/0/
        if(inach.ne.0)go to 60

         inach=1

         a=2.1038*dhor/140.0
         ar=a*hy
         ar3=ar*ar*ar
         chh=(1+ar-ar3/3)*exp(-ar)
c ---         chh=-a*a
           do 50 i=2,kz
             hz1(i)=0.5*(hz(i-1)+hz(i))
   50   continue

   60   continue

           if(ivt.eq.0)go to 199
c
c --- southern boundary
c
            IP = 1
          DO 160 I = 2, MM1
            IP = IP+1
            cxx=(stream(ip+1,k)-2*stream(ip,k)+stream(ip-1,k))/hy/hy
            cyy=(stream(ip+m,k)+(chh-2)*stream(ip,k))/hy/hy

             if(k.eq.1)then
          czz=scsig/hz(1)*(topden(ip)-sigz(2)*(stream(ip,1)-
     +    stream(ip,2))/hz1(2))
              go to 100
             endif
              
             if(k.eq.kz)then
          czz=scsig/hz(kz)*(sigz(kz)*(stream(ip,kzm1)-
     +    stream(ip,kz))/hz1(kz))
            go to 100
          endif

          czz=scsig/hz(k)*(sigz(k)*(stream(ip,k-1)-stream(ip,k))
     +     /hz1(k)-sigz(k+1)*(stream(ip,k)-stream(ip,k+1))/hz1(k+1))
  100   continue
           rs(ip)=cxx+czz

  160     CONTINUE
c
c --- northern boundary
c 
           IP = MSQ - M + 1
          do 170 i=2,mm1
            IP=IP+1
             cxx=(stream(ip+1,k)-2*stream(ip,k)+
     +       stream(ip-1,k))/hy/hy
             cyy=((chh-2)*stream(ip,k)+stream(ip+m,k))/hy/hy

             if(k.eq.1)then
          czz=scsig/hz(1)*(topden(ip)-sigz(2)*(stream(ip,1)-
     +    stream(ip,2))/hz1(2))
              go to 165
             endif

             if(k.eq.kz)then
          czz=scsig/hz(kz)*(sigz(kz)*(stream(ip,kzm1)-
     +    stream(ip,kz))/hz1(kz))
            go to 165
          endif

          czz=scsig/hz(k)*(sigz(k)*(stream(ip,k-1)-stream(ip,k))
     +    /hz1(k)-sigz(k+1)*(stream(ip,k)-stream(ip,k+1))/hz1(k+1))
  165   continue
           rs(ip)=cxx+czz

  170     CONTINUE
C
c --- western boundary
c
            IP = 1
          DO 180 J = 2, MM1
            IP = IP + M  
         cyy=(stream(ip+m,k)-2*stream(ip,k)+stream(ip-m,k))/hy/hy
         cxx=(stream(ip+1,k)+(chh-2)*stream(ip,k))/hy/hy

             if(k.eq.1)then
          czz=scsig/hz(1)*(topden(ip)-sigz(2)*(stream(ip,1)-
     +    stream(ip,2))/hz1(2))
              go to 175
             endif

             if(k.eq.kz)then
          czz=scsig/hz(kz)*(sigz(kz)*(stream(ip,kzm1)-
     +    stream(ip,kz))/hz1(kz))
            go to 175
          endif

        czz=scsig/hz(k)*(sigz(k)*(stream(ip,k-1)-stream(ip,k))
     +  /hz1(k)-sigz(k+1)*(stream(ip,k)-stream(ip,k+1))/hz1(k+1))
  175   continue
           rs(ip)=cyy+czz

  180     CONTINUE
c
c --- eastern boundary
c
            IP = M
          do 190 j=2,mm1  
            IP=IP+M
         cyy=(stream(ip+m,k)-2*stream(ip,k)+stream(ip-m,k))/hy/hy
         cxx=((chh-2)*stream(ip,k)+stream(ip-1,k))/hy/hy

             if(k.eq.1)then
          czz=scsig/hz(1)*(topden(ip)-sigz(2)*(stream(ip,1)-
     +    stream(ip,2))/hz1(2))
              go to 185
             endif

             if(k.eq.kz)then
          czz=scsig/hz(kz)*(sigz(kz)*(stream(ip,kzm1)-
     +    stream(ip,kz))/hz1(kz))
            go to 185
          endif

        czz=scsig/hz(k)*(sigz(k)*(stream(ip,k-1)-stream(ip,k))
     +  /hz1(k)-sigz(k+1)*(stream(ip,k)-stream(ip,k+1))/hz1(k+1))
  185   continue
           rs(ip)=cyy+czz

  190     CONTINUE
            ip=1
            rs(ip)=0.5*(rs(ip+1)+rs(ip+m))
            ip=m
            rs(ip)=0.5*(rs(ip-1)+rs(ip+m))
            ip=msq-mm1
            rs(ip)=0.5*(rs(ip-m)+rs(ip+1))
            ip=msq
            rs(ip)=0.5*(rs(ip-1)+rs(ip-m))
          CONTINUE
         RETURN
c
c ---  topdensity computation
c
  199   continue
c
c --- southern boundary
c
            IP = 1
          DO 200 I = 2, MM1
            IP = IP+1
          rs(ip)=scsig*sigz(1)*(stream(ip,1)-stream(ip,2))/hz1(2)
  200  continue
c
c --- northern boundary
c 
           IP = MSQ - M + 1
          do 220 i=2,mm1
            IP=IP+1
          rs(ip)=scsig*sigz(1)*(stream(ip,1)-stream(ip,2))/hz1(2)
  220   continue
C
c --- western boundary
c
            IP = 1
          DO 240 J = 2, MM1
            IP = IP + M
          rs(ip)=scsig*sigz(1)*(stream(ip,1)-stream(ip,2))/hz1(2)
  240   continue
c
c --- eastern boundary
c
            IP = M
          do 260 j=2,mm1  
            IP=IP+M
          rs(ip)=scsig*sigz(1)*(stream(ip,1)-stream(ip,2))/hz1(2)
  260   continue
            ip=1
            rs(ip)=0.5*(rs(ip+1)+rs(ip+m))
            ip=m
            rs(ip)=0.5*(rs(ip-1)+rs(ip+m))
            ip=msq-mm1
            rs(ip)=0.5*(rs(ip-m)+rs(ip+1))
            ip=msq
            rs(ip)=0.5*(rs(ip-1)+rs(ip-m))
  279   return

         end
