      subroutine robc
c
c=======================================================================
c                                                                    ===
c  This routine sets open boundary condition for tracers, velocity,  ===
c  tranport streamfuntion, and barotropic relative vorticity.        ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     IBDY    boundary to radiate (integer):                         ===
c             IBDY=1  western  boundary                              ===
c             IBDY=2  southern boundary                              ===
c             IBDY=3  eastern  boundary                              ===
c             IBDY=4  northern boundary                              ===
c     J       current slab number (integer).                         ===
c                                                                    ===
c  Entries:                                                          ===
c                                                                    ===
c  TROBC_EXT   One-time level, spatial extrapolation on tracers.     ===
c  TROBC_MOI   Modified Orlanski radiation (implicit) on tracers.    ===
c  TROBC_ORI   Orlanski radiation (implicit) on tracers.             ===
c  VROBC_EXT   One-time level, spatial extrapolation on velocity.    ===
c  VROBC_MOI   Modified Orlanski radiation (implicit) on velocity.   ===
c  VROBC_ORI   Orlanski radiation (implicit) on velocity.            ===
c  VROBC_SR    Spall and Robinson boundary conditions on velocity.   ===
c  PROBC_MOI   Modified Orlanski radiation (implicit) on transport.  ===
c  PROBC_ORI   Orlanski radiation (implicit) on transport.           ===
c  ZROBC_ORI   Orlanski radiation (implicit) on vorticity tendency.  ===
c  ZROBC_EXT   Temporal extrapolation on vorticity.                  ===
c  ZROBC_RPH   Reduced physics on vorticity.                         ===
c  ZROBC_SR    Spall and Robinson boundary conditions on vorticity.  ===
#ifdef coast
c  Call        IND_BDY                                               ===
#endif
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <scalar.h>
#include <onedim.h>
#include <fields.h>
#include <workspa.h>
#include <vertslabs.h>
#include <bndata.h>
#include <fullwd.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,ib,ibdy,j,jb,k,n,ic
      FLOAT
     *      centz,delbar,delbar1,delbar2,dzdx,dzdy,eastz,fxa,fxb,fxc,mu,
     *      small,soutz,uvel,utot,vtot,vvel,westz,zfe,zfn,zfs,zfw,znrth,
     *      zpc,zpe,zpn,zps,zpw,uub_ibic,uub_ib,uu_ibic,
     *      vvb_ibic,vvb_ib,vv_ibic,uub_i,uubm_i,uu_i,vvb_i,vvbm_i,vv_i,
     *      uubp_i,vvbp_i
      FLOAT
     *      deltax(km),deltay(km),zetanew(km),zetaold(km)
      parameter (small=c1em35)

#ifdef coast
      integer l, ind_b(4,2,mcseg)
#endif

c
c=======================================================================
c  Set Orlanski Radiation Implicit (ORI) open boundary conditions
c  on tracers. 
c=======================================================================
c
      entry trobc_ori (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        ic=1
        do 100 n=1,nt
        do 100 k=1,km
          fxa=(tb(ib+ic,k,n)-ta(ib+ic,k,n))
          fxb=ta(ib+ic,k,n)+tb(ib+ic,k,n)-c2*t(ib+2*ic,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                ta(ib,k,n)=tb(ib,k,n)
              else
                mu=fxa/fxb
                ta(ib,k,n)=((c1-mu)*tb(ib,k,n)+c2*mu*t(ib+ic,k,n))
     *                     /(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ta(ib,k,n)=fxc*tb(ib,k,n)+(c1-fxc)*t(ib+ic,k,n)
#endif
            else
              ta(ib,k,n)=t(ib+ic,k,n)
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ta(ib,k,n)=t(ib+ic,k,n)
#endif
          else
            ta(ib,k,n)=tb(ib,k,n)
          endif
 100    continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt
        ic=-1
        do 110 n=1,nt
        do 110 k=1,km
          fxa=(tb(ib+ic,k,n)-ta(ib+ic,k,n))
          fxb=ta(ib+ic,k,n)+tb(ib+ic,k,n)-c2*t(ib+2*ic,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                ta(ib,k,n)=tb(ib,k,n)
              else
                mu=fxa/fxb
                ta(ib,k,n)=((c1-mu)*tb(ib,k,n)+c2*mu*t(ib+ic,k,n))
     *                     /(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ta(ib,k,n)=fxc*tb(ib,k,n)+(c1-fxc)*t(ib+ic,k,n)
#endif
            else
              ta(ib,k,n)=t(ib+ic,k,n)
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ta(ib,k,n)=t(ib+ic,k,n)
#endif
          else
            ta(ib,k,n)=tb(ib,k,n)
          endif
 110    continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        do 130 n=1,nt
        do 130 k=1,km
#ifdef pjhcrnr
          fxa=(tb(2,k,n)-ta(2,k,n))
          fxb=ta(2,k,n)+tb(2,k,n)-c2*tp(3,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ta(1,k,n)=tbm(1,k,n)
              else
                mu=fxa/fxb
                ta(1,k,n)=((c1-mu)*tbm(1,k,n)+c2*mu*t(2,k,n))/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ta(1,k,n)=fxc*tbm(1,k,n)+(c1-fxc)*t(2,k,n)
# endif
            else
              ta(1,k,n)=t(2,k,n)
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ta(1,k,n)=t(2,k,n)
# endif
          else
            ta(1,k,n)=tbm(1,k,n)
          endif
c
          fxa=(tb(imtm1,k,n)-ta(imtm1,k,n))
          fxb=ta(imtm1,k,n)+tb(imtm1,k,n)-c2*tp(imtm2,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ta(imt,k,n)=tbm(imt,k,n)
              else
                mu=fxa/fxb
                ta(imt,k,n)=((c1-mu)*tbm(imt,k,n)+c2*mu*t(imtm1,k,n))
     &                      /(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ta(imt,k,n)=fxc*tbm(imt,k,n)+(c1-fxc)*t(imtm1,k,n)
# endif
            else
              ta(imt,k,n)=t(imtm1,k,n)
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ta(imt,k,n)=t(imtm1,k,n)
# endif
          else
            ta(imt,k,n)=tbm(imt,k,n)
          endif
#endif
          do 120 i=2,imtm1
            fxa=(tb(i,k,n)-ta(i,k,n))
            fxb=ta(i,k,n)+tb(i,k,n)-c2*tp(i,k,n)
            if(fxa*fxb.gt.c0) then
              if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
                if(abs(fxa).lt.small) then
                  ta(i,k,n)=tbm(i,k,n)
                else
                  mu=fxa/fxb
                  ta(i,k,n)=((c1-mu)*tbm(i,k,n)+c2*mu*t(i,k,n))/(c1+mu)
                endif
#else
                mu=fxa/fxb
                fxc=(c1-mu)/(c1+mu)
                ta(i,k,n)=fxc*tbm(i,k,n)+(c1-fxc)*t(i,k,n)
#endif
              else
                ta(i,k,n)=t(i,k,n)
              endif 
#ifdef pjhori
            elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
              ta(i,k,n)=t(i,k,n)
#endif
            else
              ta(i,k,n)=tbm(i,k,n)
            endif
 120      continue
#ifndef pjhcrnr
          ta(1,k,n)=p5*(ta(1,k,n)+ta(2,k,n))
          ta(imt,k,n)=p5*(ta(imt,k,n)+ta(imtm1,k,n))
#endif
 130    continue    
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        do 150 n=1,nt
        do 150 k=1,km
#ifdef pjhcrnr
          fxa=(tb(2,k,n)-ta(2,k,n))
          fxb=ta(2,k,n)+tb(2,k,n)-c2*tm(3,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ta(1,k,n)=tbp(1,k,n)
              else
                mu=fxa/fxb
                ta(1,k,n)=((c1-mu)*tbp(1,k,n)+c2*mu*t(2,k,n))/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ta(1,k,n)=fxc*tbp(1,k,n)+(c1-fxc)*t(2,k,n)
# endif
            else
              ta(1,k,n)=t(2,k,n)
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ta(1,k,n)=t(2,k,n)
# endif
          else
            ta(1,k,n)=tbp(1,k,n)
          endif
c
          fxa=(tb(imtm1,k,n)-ta(imtm1,k,n))
          fxb=ta(imtm1,k,n)+tb(imtm1,k,n)-c2*tm(imtm2,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ta(imt,k,n)=tbp(imt,k,n)
              else
                mu=fxa/fxb
                ta(imt,k,n)=((c1-mu)*tbp(imt,k,n)+c2*mu*t(imtm1,k,n))
     &                      /(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ta(imt,k,n)=fxc*tbp(imt,k,n)+(c1-fxc)*t(imtm1,k,n)
# endif
            else
              ta(imt,k,n)=t(imtm1,k,n)
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ta(imt,k,n)=t(imtm1,k,n)
# endif
          else
            ta(imt,k,n)=tbp(imt,k,n)
          endif
#endif
          do 140 i=2,imtm1
            fxa=(tb(i,k,n)-ta(i,k,n))
            fxb=ta(i,k,n)+tb(i,k,n)-c2*tm(i,k,n)
            if(fxa*fxb.gt.c0) then
              if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
                if(abs(fxa).lt.small) then
                  ta(i,k,n)=tbp(i,k,n)
                else
                  mu=fxa/fxb
                  ta(i,k,n)=((c1-mu)*tbp(i,k,n)+c2*mu*t(i,k,n))/(c1+mu)
                endif
#else
                mu=fxa/fxb
                fxc=(c1-mu)/(c1+mu)
                ta(i,k,n)=fxc*tbp(i,k,n)+(c1-fxc)*t(i,k,n)
#endif
              else
                ta(i,k,n)=t(i,k,n)
              endif 
#ifdef pjhori
            elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
              ta(i,k,n)=t(i,k,n)
#endif
            else
              ta(i,k,n)=tbp(i,k,n)
            endif
 140      continue
#ifndef pjhcrnr
          ta(1,k,n)=p5*(ta(1,k,n)+ta(2,k,n))
          ta(imt,k,n)=p5*(ta(imt,k,n)+ta(imtm1,k,n))
#endif
 150    continue    
      endif
      return
c
c=======================================================================
c  Set Modified Orlanski Implicit (MOI) radiation open boundary
c  conditions on tracers. 
c=======================================================================
c
      entry trobc_moi (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        ic=1
        do 200 n=1,nt
        do 200 k=1,km
          fxa=tb(ib+ic,k,n)-ta(ib+ic,k,n)
          fxb=ta(ib+ic,k,n)+tb(ib+ic,k,n)-c2*t(ib+2*ic,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              ta(ib,k,n)=tb(ib,k,n)
            else
              ta(ib,k,n)=t(ib+ic,k,n)
            endif 
          else
            ta(ib,k,n)=tb(ib,k,n)
          endif
 200    continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt
        ic=-1
        do 210 n=1,nt
        do 210 k=1,km
          fxa=tb(ib+ic,k,n)-ta(ib+ic,k,n)
          fxb=ta(ib+ic,k,n)+tb(ib+ic,k,n)-c2*t(ib+2*ic,k,n)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              ta(ib,k,n)=tb(ib,k,n)
            else
              ta(ib,k,n)=t(ib+ic,k,n)
            endif 
          else
            ta(ib,k,n)=tb(ib,k,n)
          endif
 210    continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        do 230 n=1,nt
        do 230 k=1,km
          do 220 i=2,imtm1
            fxa=tb(i,k,n)-ta(i,k,n)
            fxb=ta(i,k,n)+tb(i,k,n)-c2*tp(i,k,n)
            if(fxa*fxb.gt.c0) then
              if(abs(fxb).gt.abs(fxa)) then
                ta(i,k,n)=tbm(i,k,n)
              else
                ta(i,k,n)=t(i,k,n)
              endif 
            else
              ta(i,k,n)=tbm(i,k,n)
            endif
 220      continue
          ta(1,k,n)=p5*(ta(1,k,n)+ta(2,k,n))
          ta(imt,k,n)=p5*(ta(imt,k,n)+ta(imtm1,k,n))
 230    continue    
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        do 250 n=1,nt
        do 250 k=1,km
          do 240 i=2,imtm1
            fxa=tb(i,k,n)-ta(i,k,n)
            fxb=ta(i,k,n)+tb(i,k,n)-c2*tm(i,k,n)
            if(fxa*fxb.gt.c0) then
              if(abs(fxb).gt.abs(fxa)) then
                ta(i,k,n)=tbp(i,k,n)
              else
                ta(i,k,n)=t(i,k,n)
              endif 
            else
              ta(i,k,n)=tbp(i,k,n)
            endif
 240      continue
          ta(1,k,n)=p5*(ta(1,k,n)+ta(2,k,n))
          ta(imt,k,n)=p5*(ta(imt,k,n)+ta(imtm1,k,n))
 250    continue    
      endif
      return
c
c=======================================================================
c  Set one time level, spatial extrapolation boundary conditions on
c  tracers.
c=======================================================================
c
      entry trobc_ext (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        do 300 n=1,nt
        do 300 k=1,km
          utot=p5*(u(ib,k)+um(ib,k))
          if(utot.lt.c0) then    
            ta(ib,k,n)=ta(ib+1,k,n)
          else
            ta(ib,k,n)=tb(ib,k,n)
          endif
 300    continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt
        do 310 n=1,nt
        do 310 k=1,km
          utot=p5*(u(ib-1,k)+um(ib-1,k))
          if(utot.gt.c0) then
            ta(ib,k,n)=ta(ib-1,k,n)
          else
            ta(ib,k,n)=tb(ib,k,n)
          endif
 310    continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        do 330 n=1,nt
        do 330 k=1,km
          do 320 i=2,imtm1
            vtot = p5*(vm(i-1,k)+vm(i,k))
            if(vtot.lt.c0) then
              ta(i,k,n)=ta(i,k,n)
            else
              ta(i,k,n)=tbm(i,k,n)
            endif
 320     continue    
         ta(1,k,n)=p5*(ta(1,k,n)+ta(2,k,n))
         ta(imt,k,n)=p5*(ta(imt,k,n)+ta(imtm1,k,n))
 330   continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        do 350 n=1,nt
        do 350 k=1,km
          do 340 i=2,imtm1
            vtot=p5*(v(i-1,k)+v(i,k))
            if(vtot.gt.c0) then
              ta(i,k,n)=ta(i,k,n)
            else
              ta(i,k,n)=tbp(i,k,n)
            endif
 340      continue    
          ta(1,k,n)=p5*(ta(1,k,n)+ta(2,k,n))
          ta(imt,k,n)=p5*(ta(imt,k,n)+ta(imtm1,k,n))
 350    continue
      endif
      return
c
c=======================================================================
c  Set one time level, spatial extrapolation boundary conditions on
c  velocity.
c=======================================================================
c
      entry vrobc_ext (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        do 400 k=1,km
          utot=u(ib+1,k)
          if(utot.lt.c0) then    
            ua(ib,k)=c2*(u (ib+1,k)-ubar (ib+1))-
     *                  (ub(ib+2,k)-ubarb(ib+2))
            va(ib,k)=va(ib+1,k)
          else
            ua(ib,k)=c2*(ub(ib+1,k)-ubarb(ib+1))-
     *                  (ub(ib+2,k)-ubarb(ib+2))
            va(ib,k)=c2*(vb(ib+1,k)-vbarb(ib+1))-
     *                  (vb(ib+2,k)-vbarb(ib+2))
          endif
 400    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 402 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 402   continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 403 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 403   continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imtm1
        do 410 k=1,km
          utot=u(ib-1,k)
          if(utot.gt.c0) then
            ua(ib,k)=c2*(u (ib-1,k)-ubar (ib-1))-
     *                  (ub(ib-2,k)-ubarb(ib-2))
            va(ib,k)=va(ib-1,k)
          else
            ua(ib,k)=c2*(ub(ib-1,k)-ubarb(ib-1))-
     *                  (ub(ib-2,k)-ubarb(ib-2))
            va(ib,k)=c2*(vb(ib-1,k)-vbarb(ib-1))-
     *                  (vb(ib-2,k)-vbarb(ib-2))
          endif
 410    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 412 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 412    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 414 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 414   continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        do 430 k=1,km
          do 420 i=2,imtm2
            vtot=v(i,k)
            if(vtot.lt.c0) then
              ua(i,k)=ua(i,k)
              va(i,k)=c2*(v (i,k)-vbar (i))-(vb (i,k)-vbarb(i))
            else
              ua(i,k)=c2*(ub(i,k)-ubarb(i))-(ubp(i,k)-ubarbp(i))
              va(i,k)=c2*(vb(i,k)-vbarb(i))-(vbp(i,k)-vbarbp(i))
            endif
 420      continue
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 430    continue    
c
c  Determine and subtract incorrect vertical mean.
c
        do 436 i=1,imtm1
        fxa=c0
        fxb=c0
        do 432 k=1,km
           fxa=fxa+ua(i,k)*dzvqz(i,k,0)
           fxb=fxb+va(i,k)*dzvqz(i,k,0)
 432    continue
        fxa=fxa*hvav(i,j)
        fxb=fxb*hvav(i,j)
        do 434 k=1,km
           ua(i,k)=ua(i,k)-fxa
           va(i,k)=va(i,k)-fxb
 434    continue
 436    continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        do 450 k=1,km
          do 440 i=2,imtm2
            vtot=v(i,k)
            if(vtot.gt.c0) then
              ua(i,k)=ua(i,k)
              va(i,k)=c2*(v (i,k)-vbar (i))-(vb (i,k)-vbarb(i))
            else
              ua(i,k)=c2*(ub(i,k)-ubarb(i))-(ubm(i,k)-ubarbm(i))
              va(i,k)=c2*(vb(i,k)-vbarb(i))-(vbm(i,k)-vbarbm(i))
            endif
 440      continue
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 450    continue    
c
c  Determine and subtract incorrect vertical mean.
c
        do 456 i=1,imtm1
        fxa=c0
        fxb=c0
        do 452 k=1,km
           fxa=fxa+ua(i,k)*dzvqz(i,k,0)
           fxb=fxb+va(i,k)*dzvqz(i,k,0)
 452    continue
        fxa=fxa*hvav(i,j)
        fxb=fxb*hvav(i,j)
        do 454 k=1,km
           ua(i,k)=ua(i,k)-fxa
           va(i,k)=va(i,k)-fxb
 454    continue
 456    continue
c
      endif
      return
c
c=======================================================================
c  Set Orlanski Radiation Implicit (ORI) open boundary conditions
c  on velocities. 
c=======================================================================
c
      entry vrobc_ori (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        ic=1
        do 500 k=1,km
          uub_ibic=ub(ib+ic,k)-ubarb(ib+ic)
          uub_ib=ub(ib,k)-ubarb(ib)
          uu_ibic=u(ib+ic,k)-ubar(ib+ic)
c
          fxa=(uub_ibic-ua(ib+ic,k))
          fxb=ua(ib+ic,k) + uub_ibic -
     &         c2*(u(ib+2*ic,k)-ubar(ib+2*ic))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(ib,k)=uub_ib
              else
                mu=fxa/fxb
                ua(ib,k)=((c1-mu)*uub_ib+c2*mu*uu_ibic)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(ib,k)=fxc*uub_ib+(c1-fxc)*uu_ibic
#endif
            else
              ua(ib,k)=uu_ibic
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(ib,k)=uu_ibic
#endif
          else
            ua(ib,k)=uub_ib
          endif
 500    continue
        do 501 k=1,km
          vvb_ibic=vb(ib+ic,k)-vbarb(ib+ic)
          vvb_ib=vb(ib,k)-vbarb(ib)
          vv_ibic=v(ib+ic,k)-vbar(ib+ic)
c
          fxa=(vvb_ibic-va(ib+ic,k))
          fxb=va(ib+ic,k) + vvb_ibic -
     &         c2*(v(ib+2*ic,k)-vbar(ib+2*ic))

          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                va(ib,k)=vvb_ib
              else
                mu=fxa/fxb
                va(ib,k)=((c1-mu)*vvb_ib+c2*mu*vv_ibic)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(ib,k)=fxc*vvb_ib+(c1-fxc)*vv_ibic
#endif
            else
              va(ib,k)=vv_ibic
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(ib,k)=vv_ibic
#endif
          else
            va(ib,k)=vvb_ib
          endif
 501    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 510 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 510    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 520 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 520    continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imu
        ic=-1
        do 530 k=1,km
          uub_ibic=ub(ib+ic,k)-ubarb(ib+ic)
          uub_ib=ub(ib,k)-ubarb(ib)
          uu_ibic=u(ib+ic,k)-ubar(ib+ic)
c
          fxa=(uub_ibic-ua(ib+ic,k))
          fxb=ua(ib+ic,k) + uub_ibic -
     &         c2*(u(ib+2*ic,k)-ubar(ib+2*ic))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(ib,k)=uub_ib
              else
                mu=fxa/fxb
                ua(ib,k)=((c1-mu)*uub_ib+c2*mu*uu_ibic)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(ib,k)=fxc*uub_ib+(c1-fxc)*uu_ibic
#endif
            else
              ua(ib,k)=uu_ibic
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(ib,k)=uu_ibic
#endif
          else
            ua(ib,k)=uub_ib
          endif
 530    continue
        do 531 k=1,km
          vvb_ibic=vb(ib+ic,k)-vbarb(ib+ic)
          vvb_ib=vb(ib,k)-vbarb(ib)
          vv_ibic=v(ib+ic,k)-vbar(ib+ic)
c
          fxa=(vvb_ibic-va(ib+ic,k))
          fxb=va(ib+ic,k) + vvb_ibic -
     &         c2*(v(ib+2*ic,k)-vbar(ib+2*ic))

          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                va(ib,k)=vvb_ib
              else
                mu=fxa/fxb
                va(ib,k)=((c1-mu)*vvb_ib+c2*mu*vv_ibic)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(ib,k)=fxc*vvb_ib+(c1-fxc)*vv_ibic
#endif
            else
              va(ib,k)=vv_ibic
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(ib,k)=vv_ibic
#endif
          else
            va(ib,k)=vvb_ib
          endif
 531    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 540 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 540    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 550 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 550    continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
c  Zonal Component.
c
        do 560 k=1,km
#ifdef pjhcrnr
c
c    SW corner
c
          uub_i=ub(2,k)-ubarb(2)
          uubm_i=ubm(1,k)-ubarbm(1)
          uu_i=u(2,k)-ubar(2)
c
          fxa=(uub_i-ua(2,k))
          fxb=ua(2,k) + uub_i -c2*(up(3,k)-ubarp(3))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(1,k)=uubm_i
              else
                mu=fxa/fxb
                ua(1,k)=((c1-mu)*uubm_i+c2*mu*uu_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(1,k)=fxc*uubm_i+(c1-fxc)*uu_i
# endif
            else
              ua(1,k)=uu_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(1,k)=uu_i
# endif
          else
            ua(1,k)=uubm_i
          endif
c
c    SE corner
c
          uub_i=ub(imum1,k)-ubarb(imum1)
          uubm_i=ubm(imu,k)-ubarbm(imu)
          uu_i=u(imum1,k)-ubar(imum1)
c
          fxa=(uub_i-ua(imum1,k))
          fxb=ua(imum1,k) + uub_i -c2*(up(imum2,k)-ubarp(imum2))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(imu,k)=uubm_i
              else
                mu=fxa/fxb
                ua(imu,k)=((c1-mu)*uubm_i+c2*mu*uu_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(imu,k)=fxc*uubm_i+(c1-fxc)*uu_i
# endif
            else
              ua(imu,k)=uu_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(imu,k)=uu_i
# endif
          else
            ua(imu,k)=uubm_i
          endif
c
c    Interior southern boundary
c
#endif
        do 560 i=2,imum1
          uub_i=ub(i,k)-ubarb(i)
          uubm_i=ubm(i,k)-ubarbm(i)
          uu_i=u(i,k)-ubar(i)
c
          fxa=(uub_i-ua(i,k))
          fxb=ua(i,k) + uub_i -c2*(up(i,k)-ubarp(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(i,k)=uubm_i
              else
                mu=fxa/fxb
                ua(i,k)=((c1-mu)*uubm_i+c2*mu*uu_i)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(i,k)=fxc*uubm_i+(c1-fxc)*uu_i
#endif
            else
              ua(i,k)=uu_i
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(i,k)=uu_i
#endif
          else
            ua(i,k)=uubm_i
          endif
 560    continue
c
c  Meridional Component.
c
        do 561 k=1,km
#ifdef pjhcrnr
c
c    SW corner
c
          vvb_i=vb(2,k)-vbarb(2)
          vvbm_i=vbm(1,k)-vbarbm(1)
          vv_i=v(2,k)-vbar(2)
c
          fxa=(vvb_i-va(2,k))
          fxb=va(2,k)+vvb_i-c2*(vp(3,k)-vbarp(3))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                va(1,k)=vvbm_i
              else
                mu=fxa/fxb
                va(1,k)=((c1-mu)*vvbm_i+c2*mu*vv_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(1,k)=fxc*vvbm_i+(c1-fxc)*vv_i
# endif
            else
              va(1,k)=vv_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(1,k)=vv_i
# endif
          else
            va(1,k)=vvbm_i
          endif
c
c    SE corner
c
          vvb_i=vb(imum1,k)-vbarb(imum1)
          vvbm_i=vbm(imu,k)-vbarbm(imu)
          vv_i=v(imum1,k)-vbar(imum1)
c
          fxa=(vvb_i-va(imum1,k))
          fxb=va(imum1,k)+vvb_i-c2*(vp(imum2,k)-vbarp(imum2))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                va(imu,k)=vvbm_i
              else
                mu=fxa/fxb
                va(imu,k)=((c1-mu)*vvbm_i+c2*mu*vv_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(imu,k)=fxc*vvbm_i+(c1-fxc)*vv_i
# endif
            else
              va(imu,k)=vv_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(imu,k)=vv_i
# endif
          else
            va(imu,k)=vvbm_i
          endif
c
c    Interior southern boundary
c
#endif
        do 561 i=2,imum1
          vvb_i=vb(i,k)-vbarb(i)
          vvbm_i=vbm(i,k)-vbarbm(i)
          vv_i=v(i,k)-vbar(i)
c
          fxa=(vvb_i-va(i,k))
          fxb=va(i,k)+vvb_i-c2*(vp(i,k)-vbarp(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                va(i,k)=vvbm_i
              else
                mu=fxa/fxb
                va(i,k)=((c1-mu)*vvbm_i+c2*mu*vv_i)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(i,k)=fxc*vvbm_i+(c1-fxc)*vv_i
#endif
            else
              va(i,k)=vv_i
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(i,k)=vv_i
#endif
          else
            va(i,k)=vvbm_i
          endif
 561    continue
#ifndef pjhcrnr
c
c  Set end points.
c
        do 570 k=1,km
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 570    continue
#endif
c
c  Determine and subtract incorrect vertical mean.
c
        do 600 i=1,imtm1
           fxa=c0
           fxb=c0
           do 580 k=1,km
              fxa=fxa+ua(i,k)*dzvqz(i,k,0)
              fxb=fxb+va(i,k)*dzvqz(i,k,0)
 580       continue
           fxa=fxa*hvav(i,j)
           fxb=fxb*hvav(i,j)
           do 590 k=1,km
              ua(i,k)=ua(i,k)-fxa
              va(i,k)=va(i,k)-fxb
 590       continue
 600    continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
c  Zonal Component.
c
        do 610 k=1,km
#ifdef pjhcrnr
c
c    NW corner
c
          uub_i=ub(2,k)-ubarb(2)
          uubp_i=ubp(1,k)-ubarbp(1)
          uu_i=u(2,k)-ubar(2)
c
          fxa=(uub_i-ua(2,k))
          fxb=ua(2,k)+uub_i-c2*(um(3,k)-ubarm(3))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(1,k)=uubp_i
              else
                mu=fxa/fxb
                ua(1,k)=((c1-mu)*uubp_i+c2*mu*uu_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(1,k)=fxc*uubp_i+(c1-fxc)*uu_i
# endif
            else
              ua(1,k)=uu_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(1,k)=uu_i
# endif
          else
            ua(1,k)=uubp_i
          endif
c
c    NE corner
c
          uub_i=ub(imum1,k)-ubarb(imum1)
          uubp_i=ubp(imu,k)-ubarbp(imu)
          uu_i=u(imum1,k)-ubar(imum1)
c
          fxa=(uub_i-ua(imum1,k))
          fxb=ua(imum1,k)+uub_i-c2*(um(imum2,k)-ubarm(imum2))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(imu,k)=uubp_i
              else
                mu=fxa/fxb
                ua(imu,k)=((c1-mu)*uubp_i+c2*mu*uu_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(imu,k)=fxc*uubp_i+(c1-fxc)*uu_i
# endif
            else
              ua(imu,k)=uu_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(imu,k)=uu_i
# endif
          else
            ua(imu,k)=uubp_i
          endif
c
c    Interior northern boundary
c
#endif
        do 610 i=2,imum1
          uub_i=ub(i,k)-ubarb(i)
          uubp_i=ubp(i,k)-ubarbp(i)
          uu_i=u(i,k)-ubar(i)
c
          fxa=(uub_i-ua(i,k))
          fxb=ua(i,k)+uub_i-c2*(um(i,k)-ubarm(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                ua(i,k)=uubp_i
              else
                mu=fxa/fxb
                ua(i,k)=((c1-mu)*uubp_i+c2*mu*uu_i)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              ua(i,k)=fxc*uubp_i+(c1-fxc)*uu_i
#endif
            else
              ua(i,k)=uu_i
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            ua(i,k)=uu_i
#endif
          else
            ua(i,k)=uubp_i
          endif
 610    continue
c
c  Meridional Component.
c
        do 611 k=1,km
#ifdef pjhcrnr
c
c    NW corner
c
          vvb_i=vb(2,k)-vbarb(2)
          vvbp_i=vbp(1,k)-vbarbp(1)
          vv_i=v(2,k)-vbar(2)
c 
          fxa=(vvb_i-va(2,k))
          fxb=va(2,k)+vvb_i-c2*(vm(3,k)-vbarm(3))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                va(1,k)=vvbp_i
              else
                mu=fxa/fxb
                va(1,k)=((c1-mu)*vvbp_i+c2*mu*vv_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(1,k)=fxc*vvbp_i+(c1-fxc)*vv_i
# endif
            else
              va(1,k)=vv_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(1,k)=vv_i
# endif
          else
            va(1,k)=vvbp_i
          endif
c
c    NE corner
c
          vvb_i=vb(imum1,k)-vbarb(imum1)
          vvbp_i=vbp(imu,k)-vbarbp(imu)
          vv_i=v(imum1,k)-vbar(imum1)
c 
          fxa=(vvb_i-va(imum1,k))
          fxb=va(imum1,k)+vvb_i-c2*(vm(imum2,k)-vbarm(imum2))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small) then
                va(imu,k)=vvbp_i
              else
                mu=fxa/fxb
                va(imu,k)=((c1-mu)*vvbp_i+c2*mu*vv_i)/(c1+mu)
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(imu,k)=fxc*vvbp_i+(c1-fxc)*vv_i
# endif
            else
              va(imu,k)=vv_i
            endif 
# ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(imu,k)=vv_i
# endif
          else
            va(imu,k)=vvbp_i
          endif
c
c    Interior northern boundary
c
#endif
        do 611 i=2,imum1
          vvb_i=vb(i,k)-vbarb(i)
          vvbp_i=vbp(i,k)-vbarbp(i)
          vv_i=v(i,k)-vbar(i)
c 
          fxa=(vvb_i-va(i,k))
          fxb=va(i,k)+vvb_i-c2*(vm(i,k)-vbarm(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small) then
                va(i,k)=vvbp_i
              else
                mu=fxa/fxb
                va(i,k)=((c1-mu)*vvbp_i+c2*mu*vv_i)/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              va(i,k)=fxc*vvbp_i+(c1-fxc)*vv_i
#endif
            else
              va(i,k)=vv_i
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            va(i,k)=vv_i
#endif
          else
            va(i,k)=vvbp_i
          endif
 611    continue
#ifndef pjhcrnr
c
c  Set end points.
c
        do 620 k=1,km
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 620   continue
#endif
c
c  Determine and subtract incorrect vertical mean.
c
        do 650 i=1,imtm1
          fxa=c0
          fxb=c0
          do 630 k=1,km
            fxa=fxa+ua(i,k)*dzvqz(i,k,0)
            fxb=fxb+va(i,k)*dzvqz(i,k,0)
 630      continue
          fxa=fxa*hvav(i,j)
          fxb=fxb*hvav(i,j)
          do 640 k=1,km
            ua(i,k)=ua(i,k)-fxa
            va(i,k)=va(i,k)-fxb
 640      continue
 650    continue
      endif
      return
c
c=======================================================================
c  Set Modified Orlanski Implicit (MOI) radiation open boundary
c  conditions on internal velocities. 
c=======================================================================
c
      entry vrobc_moi (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        ic=1
        do 700 k=1,km
          uub_ibic=ub(ib+ic,k)-ubarb(ib+ic)
          uub_ib=ub(ib,k)-ubarb(ib)
          uu_ibic=u(ib+ic,k)-ubar(ib+ic)
c
          fxa=c2*(uub_ibic-ua(ib+ic,k))
          fxb=ua(ib+ic,k) + uub_ibic -
     &         c2*(u(ib+2*ic,k)-ubar(ib+2*ic))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              ua(ib,k)=uub_ib
            else
              ua(ib,k)=uu_ibic
            endif 
          else
            ua(ib,k)=uub_ibic
          endif
 700    continue
c
        do 701 k=1,km
          vvb_ibic=vb(ib+ic,k)-vbarb(ib+ic)
          vvb_ib=vb(ib,k)-vbarb(ib)
          vv_ibic=v(ib+ic,k)-vbar(ib+ic)
c
          fxa=c2*(vvb_ibic-va(ib+ic,k))
          fxb=va(ib+ic,k) + vvb_ibic -
     &         c2*(v(ib+2*ic,k)-vbar(ib+2*ic))

          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              va(ib,k)=vvb_ib
            else
              va(ib,k)=vv_ibic
            endif 
          else
            va(ib,k)=vvb_ibic
          endif
 701    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 710 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 710    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 720 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 720    continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imu
        ic=-1
        do 730 k=1,km
          uub_ibic=ub(ib+ic,k)-ubarb(ib+ic)
          uub_ib=ub(ib,k)-ubarb(ib)
          uu_ibic=u(ib+ic,k)-ubar(ib+ic)
c
          fxa=c2*(uub_ibic-ua(ib+ic,k))
          fxb=ua(ib+ic,k) + uub_ibic -
     &         c2*(u(ib+2*ic,k)-ubar(ib+2*ic))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              ua(ib,k)=uub_ib
            else
              ua(ib,k)=uu_ibic
            endif 
          else
            ua(ib,k)=uub_ibic
          endif
 730    continue
c
        do 731 k=1,km
          vvb_ibic=vb(ib+ic,k)-vbarb(ib+ic)
          vvb_ib=vb(ib,k)-vbarb(ib)
          vv_ibic=v(ib+ic,k)-vbar(ib+ic)
c
          fxa=c2*(vvb_ibic-va(ib+ic,k))
          fxb=va(ib+ic,k) + vvb_ibic -
     &         c2*(v(ib+2*ic,k)-vbar(ib+2*ic))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              va(ib,k)=vvb_ib
            else
              va(ib,k)=vv_ibic
            endif 
          else
            va(ib,k)=vvb_ibic
          endif
 731    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 740 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 740    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 750 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 750    continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        do 760 k=1,km
        do 760 i=2,imtm2
          uub_i=ub(i,k)-ubarb(i)
          uubm_i=ubm(i,k)-ubarbm(i)
          uu_i=u(i,k)-ubar(i)
c
          fxa=c2*(uub_i-ua(i,k))
          fxb=ua(i,k) + uub_i -c2*(up(i,k)-ubarp(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              ua(i,k)=uub_i
            else
              ua(i,k)=uu_i
            endif 
          else
            ua(i,k)=uubm_i
          endif
 760    continue
c
        do 761 k=1,km
        do 761 i=2,imtm2
          vvb_i=vb(i,k)-vbarb(i)
          vvbm_i=vbm(i,k)-vbarbm(i)
          vv_i=v(i,k)-vbar(i)
c
          fxa=c2*(vvb_i-va(i,k))
          fxb=va(i,k)+vvb_i-c2*(vp(i,k)-vbarp(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              va(i,k)=vvbm_i
            else
              va(i,k)=vv_i
            endif 
          else
            va(i,k)=vvbm_i
          endif
 761    continue
c  
c  Set end points.
c
        do 770 k=1,km
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 770    continue    
c
c  Determine and subtract incorrect vertical mean.
c
        do 800 i=1,imtm1
          fxa=c0
          fxb=c0
          do 780 k=1,km
            fxa=fxa+ua(i,k)*dzvqz(i,k,0)
            fxb=fxb+va(i,k)*dzvqz(i,k,0)
 780      continue
          fxa=fxa*hvav(i,j)
          fxb=fxb*hvav(i,j)
          do 790 k=1,km
            ua(i,k)=ua(i,k)-fxa
            va(i,k)=va(i,k)-fxb
 790      continue
 800    continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        do 810 k=1,km
        do 810 i=2,imtm2
          uub_i=ub(i,k)-ubarb(i)
          uubp_i=ubp(i,k)-ubarbp(i)
          uu_i=u(i,k)-ubar(i)
c
          fxa=c2*(uub_i-ua(i,k))
          fxb=ua(i,k)+uub_i-c2*(um(i,k)-ubarm(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              ua(i,k)=uubp_i
            else
              ua(i,k)=uu_i
            endif 
          else
            ua(i,k)=uubp_i
          endif
 810    continue
c
        do 811 k=1,km
        do 811 i=2,imtm2
          vvb_i=vb(i,k)-vbarb(i)
          vvbp_i=vbp(i,k)-vbarbp(i)
          vv_i=v(i,k)-vbar(i)
c 
          fxa=c2*(vvb_i-va(i,k))
          fxb=va(i,k)+vvb_i-c2*(vm(i,k)-vbarm(i))
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
              va(i,k)=vvb_i
            else
              va(i,k)=vv_i
            endif 
          else
            va(i,k)=vvbp_i
          endif
 811    continue
c
c  Set end points.
c
        do 820 k=1,km
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 820    continue    
c
c  Determine and subtract incorrect vertical mean.
c
        do 850 i=1,imtm1
          fxa=c0
          fxb=c0
          do 830 k=1,km
            fxa=fxa+ua(i,k)*dzvqz(i,k,0)
            fxb=fxb+va(i,k)*dzvqz(i,k,0)
 830      continue
          fxa=fxa*hvav(i,j)
          fxb=fxb*hvav(i,j)
          do 840 k=1,km
            ua(i,k)=ua(i,k)-fxa
            va(i,k)=va(i,k)-fxb
 840      continue
 850    continue
      endif
      return
c
c=======================================================================
c  Set Spall and Robinson open boundary conditions on internal mode
c  velocity (CFvN type).
c=======================================================================
c
      entry vrobc_sr (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        do 900 k=1,km
          deltax(k)=p5*(u(ib,k)+u(ib+1,k))*dtts*dxur(ib+1)*csr(j)
          if(v(ib,k).le.c0) then
            deltay(k)=dtts*dyu2r(j)*(v(ib,k)+vp(ib,k))
          else
            deltay(k)=dtts*dyu2r(j)*(v(ib,k)+vm(ib,k))
          endif
 900    continue
c
        do 910 k=1,km
          if(v(ib,k).le.c0) then
            zetaold(k)=(v (ib+1,k)-v (ib,k))*dxur(ib)*csr(j)-
     *                 (up(ib  ,k)-u (ib,k))*dyur(j)
          else
            zetaold(k)=(v (ib+1,k)-v (ib,k))*dxur(ib)*csr(j)-
     *                 (u (ib  ,k)-um(ib,k))*dyur(j)
          endif
          zetanew(k)=(c1-deltax(k))*zetaold(k)+
     *               ((v (ib+2,k)-v (ib  ,k))*dxu2r(ib+1)*csr(j)-
     *                (up(ib+1,k)-um(ib+1,k))*dyu2r(j))*deltax(k)
     *               -deltay(k)*p5*((vp(ib+1,k)+vm(ib,k)-vp(ib,k)
     *                              -vm(ib+1,k))*dxu2r(ib+1)*csr(j)
     *               -((up(ib,k)+um(ib,k)-c2*u(ib,k))*dyu2r(j)))
 910    continue
c
        do 920 k=1,km
          if(v(ib,k).le.c0) then
            va(ib,k)=va(ib+1,k)-(up(ib,k)-u(ib,k))
     *                        -dxu(ib)*cs(j)*zetanew(k)
          else
            va(ib,k)=va(ib+1,k)-(u(ib,k)-um(ib,k))
     *                         -dxu(ib)*cs(j)*zetanew(k)
          endif
 920    continue
        do 930 k=1,km
          ua(ib,k)=uo(j,k,west)
          if(deltax(k).ge.c0) then
            va(ib,k)=vo(j,k,west)
          endif
 930    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 940 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
 940    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 950 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
 950    continue
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imtm1
        do 960 k=1,km
          deltax(k)=p5*(u(ib,k)+u(ib-1,k))*dtts*dxur(ib)*csr(j)
          if(v(ib,k).le.c0) then
            deltay(k)=dtts*dyu2r(j)*(v(ib,k)+vp(ib,k))
          else
            deltay(k)=dtts*dyu2r(j)*(v(ib,k)+vm(ib,k))
          endif
 960    continue
c
        do 970 k=1,km
          if(v(ib,k).le.c0) then
            zetaold(k)=(v (ib,k)-v (ib-1,k))*dxur(ib)*csr(j)-
     *                 (up(ib,k)-u (ib  ,k))*dyur(j)
          else
            zetaold(k)=(v (ib,k)-v (ib-1,k))*dxur(ib)*csr(j)-
     *                 (u (ib,k)-um(ib  ,k))*dyur(j)
          endif
          zetanew(k)=(c1-deltax(k))*zetaold(k)+
     *               ((v (ib  ,k)-v (ib-2,k))*dxu2r(ib)*csr(j)-
     *                (up(ib-1,k)-um(ib-1,k))*dyu2r(j))*deltax(k)
     *               -deltay(k)*p5*((vp(ib,k)+vm(ib-1,k)-vp(ib-1,k)
     *                              -vm(ib,k))*dxu2r(ib)*csr(j)
     *               -((up(ib,k)+um(ib,k)-c2*u(ib,k))*dyu2r(j)))
 970    continue
c
        do 980 k=1,km
          if(v(ib,k).le.c0) then
            va(ib,k)=va(ib-1,k)+(up(ib,k)-u(ib,k))
     *                         +dxu(ib)*cs(j)*zetanew(k)
          else
            va(ib,k)=va(ib-1,k)+(u(ib,k)-um(ib,k))
     *                        +dxu(ib)*cs(j)*zetanew(k)
          endif
 980    continue
        do 990 k=1,km
          ua(ib,k)=uo(j,k,east)
          if(deltax(k).le.c0) then
            va(ib,k)=vo(j,k,east)
          endif
 990    continue
c
c  Determine and subtract incorrect vertical mean.
c
        fxa=c0
        fxb=c0
        do 1000 k=1,km
          fxa=fxa+ua(ib,k)*dzvqz(ib,k,0)
          fxb=fxb+va(ib,k)*dzvqz(ib,k,0)
1000    continue
        fxa=fxa*hvav(ib,j)
        fxb=fxb*hvav(ib,j)
        do 1010 k=1,km
          ua(ib,k)=ua(ib,k)-fxa
          va(ib,k)=va(ib,k)-fxb
1010    continue
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        jb=1
        do 1060 i=2,imtm2
          do 1020 k=1,km
            deltay(k)=p5*(vm(i,k)+v(i,k))*dtts*dyur(jb)
            if(um(i,k).le.c0) then
              deltax(k)=dtts*dxu2r(i)*(um(i,k)+um(i+1,k))*csr(jb)
            else
              deltax(k)=dtts*dxu2r(i)*(um(i,k)+um(i-1,k))*csr(jb)
            endif
1020      continue
c
          do 1030 k=1,km
            if(um(i,k).le.c0) then
              zetaold(k)=(vm(i+1,k)-vm(i,k))*dxur(i)*csr(jb)-
     *                   (u (i  ,k)-um(i,k))*dyur(jb)
            else
              zetaold(k)=(vm(i,k)-vm(i-1,k))*dxur(i)*csr(jb)-
     *                   (u (i,k)-um(i  ,k))*dyur(jb)
            endif
            zetanew(k)=(c1-deltax(k))*zetaold(k)+
     *                 ((vm(i+1,k)-vm(i-1,k))*dxur(i)*csr(jb)-
     *                  (u (i  ,k)-um(i  ,k))*dyu2r(jb))*deltax(k)
     *                 -deltay(k)*p5*((v (i+1,k)+vm(i-1,k)-v (i-1,k)-
     *                                 vm(i+1,k))*dxu2r(i)*csr(jb)
     *                 -((up(i,k)+um(i,k)-c2*u(i,k))*dyu2r(jb)))
1030      continue
c
          do 1040 k=1,km
            if(u(i,k).le.c0) then
              ua(i,k)=ua(i,k)-(vm(i+1,k)-vm(i  ,k))
     *                       +dxu(i)*cs(jb)*zetanew(k)
            else
              ua(i,k)=ua(i,k)-(vm(i  ,k)-vm(i-1,k))
     *                       +dxu(i)*cs(jb)*zetanew(k)
            endif
1040      continue
          do 1050 k=1,km
            if(deltay(k).ge.c0) then
              ua(i,k)=uo(i,k,south)
            endif
            va(i,k)=vo(i,k,south)
1050      continue
1060    continue
c
c  Set end points.
c
        do 1070 k=1,km
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
 1070  continue
c
c  Determine and subtract incorrect vertical mean.
c
          do 1095 i=1,imtm1
             fxa=c0
             fxb=c0
             do 1080 k=1,km
                fxa=fxa+ua(i,k)*dzvqz(i,k,0)
                fxb=fxb+va(i,k)*dzvqz(i,k,0)
 1080        continue
             fxa=fxa*hvav(i,jb)
             fxb=fxb*hvav(i,jb)
             do 1090 k=1,km
                ua(i,k)=ua(i,k)-fxa
                va(i,k)=va(i,k)-fxb
 1090        continue
 1095     continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        jb=jmtm1
        do 1140 i=2,imtm2
          do 1100 k=1,km
            deltay(k)=p5*(v(i,k)+vp(i,k))*dtts*dyur(jb)
            if(up(i,k).le.c0) then
              deltax(k)=dtts*dxu2r(i)*(up(i,k)+up(i+1,k))*csr(jb)
            else
              deltax(k)=dtts*dxu2r(i)*(up(i,k)+up(i-1,k))*csr(jb)
            endif
1100      continue
c
          do 1110 k=1,km
            if(up(i,k).le.c0) then
              zetaold(k)=(vp(i+1,k)-vp(i,k))*dxur(i)*csr(jb)-
     *                   (up(i  ,k)-u (i,k))*dyur(jb)
            else
              zetaold(k)=(vp(i,k)-vp(i-1,k))*dxur(i)*csr(jb)-
     *                   (up(i,k)-u (i  ,k))*dyur(jb)
            endif
            zetanew(k)=(c1-deltax(k))*zetaold(k)+
     *                 ((vp(i+1,k)-vp(i-1,k))*dxu2r(i)*csr(jb)-
     *                  (up(i  ,k)-u (i  ,k))*dyur(jb))*deltax(k)
     *                 -deltay(k)*p5*((vp(i+1,k)+v(i-1,k)-vp(i-1,k)-
     *                                 v (i+1,k))*dxu2r(i)*csr(jb)
     *                 -((up(i,k)+um(i,k)-c2*u(i,k))*dyu2r(jb)))
1110      continue
c
          do 1120 k=1,km
            if(up(i,k).le.c0) then
              ua(i,k)=u(i,k)+(vp(i+1,k)-vp(i  ,k))
     *                      -dxu(i)*cs(jb)*zetanew(k)
            else
              ua(i,k)=u(i,k)+(vp(i  ,k)-vp(i-1,k))
     *                      -dxu(i)*cs(jb)*zetanew(k)
            endif
1120      continue
          do 1130 k=1,km
            if(deltay(k).le.c0) then
              ua(i,k)=uo(i,k,north)
            endif
            va(i,k)=vo(i,k,north)
1130      continue
1140    continue
c
c  Set end points.
c
        do 1150 k=1,km
          ua(1,k)=p5*(ua(1,k)+ua(2,k))
          ua(imtm1,k)=p5*(ua(imtm1,k)+ua(imtm2,k))
          va(1,k)=p5*(va(1,k)+va(2,k))
          va(imtm1,k)=p5*(va(imtm1,k)+va(imtm2,k))
1150    continue
c
c  Determine and subtract incorrect vertical mean.
c
        do 1180 i=1,imtm1
          fxa=c0
          fxb=c0
          do 1160 k=1,km
            fxa=fxa+ua(i,k)*dzvqz(i,k,0)
            fxb=fxb+va(i,k)*dzvqz(i,k,0)
1160      continue
          fxa=fxa*hvav(i,jb)
          fxb=fxb*hvav(i,jb)
          do 1170 k=1,km
            ua(i,k)=ua(i,k)-fxa
            va(i,k)=va(i,k)-fxb
1170      continue
1180    continue
      endif
      return
c
c=======================================================================
c  Set Orlanski Radiation Implicit (ORI) open boundary conditions
c  on transport streamfunction. 
c=======================================================================
c
      entry probc_ori (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        ic=1
        fxa=(pb(ib+ic,j)-p(ib+ic,j))
        fxb=p(ib+ic,j)+pb(ib+ic,j)-c2*p(ib+2*ic,j)
        if(fxa*fxb.gt.c0) then
           if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small*abs(fxb)) then
                 po(j,1,west)=pb(ib,j)
              else
                 mu=fxa/fxb
                 po(j,1,west)=((c1-mu)*pb(ib,j)+c2*mu*p(ib+ic,j))/
     *                        (c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              po(j,1,west)=fxc*pb(ib,j)+(c1-fxc)*p(ib+ic,j)
#endif
           else
              po(j,1,west)=p(ib+ic,j)
           endif 
#ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
          po(j,1,west)=p(ib+ic,j)
#endif
        else
           po(j,1,west)=pb(ib,j)
        endif
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt
        ic=-1
        fxa=(pb(ib+ic,j)-p(ib+ic,j))
        fxb=p(ib+ic,j)+pb(ib+ic,j)-c2*p(ib+2*ic,j)
        if(fxa*fxb.gt.c0) then
           if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small*abs(fxb)) then
                 po(j,1,east)=pb(ib,j)
              else
                 mu=fxa/fxb
                 po(j,1,east)=((c1-mu)*pb(ib,j)+c2*mu*p(ib+ic,j))/
     *                        (c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              po(j,1,east)=fxc*pb(ib,j)+(c1-fxc)*p(ib+ic,j)
#endif
           else
              po(j,1,east)=p(ib+ic,j)
           endif 
#ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
          po(j,1,east)=p(ib+ic,j)
#endif
        else
           po(j,1,east)=pb(ib,j)
        endif
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
#ifndef pjhcrnr
        do 1200 i=1,imt
#else
        do 1200 i=2,imtm1
#endif
          fxa=(pb(i,2)-p(i,2))
          fxb=p(i,2)+pb(i,2)-c2*p(i,3)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small*abs(fxb)) then
                po(i,1,south)=pb(i,1)
              else
                mu=fxa/fxb
                po(i,1,south)=((c1-mu)*pb(i,1)+c2*mu*p(i,2))/(c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              po(i,1,south)=fxc*pb(i,1)+(c1-fxc)*p(i,2)
#endif
            else
              po(i,1,south)=p(i,2)
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
             po(i,1,south)=p(i,2)
#endif
          else
             po(i,1,south)=pb(i,1)
          endif
1200    continue    
c
#ifdef pjhcrnr
c  SW corner
c
        fxa=(pb(2,2)-p(2,2))
        fxb=p(2,2)+pb(2,2)-c2*p(3,3)
        if(fxa*fxb.gt.c0) then
          if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
            if(abs(fxa).lt.small*abs(fxb)) then
              po(1,1,south)=pb(1,1)
            else
              mu=fxa/fxb
              po(1,1,south)=((c1-mu)*pb(1,1)+c2*mu*p(2,2))/(c1+mu)
            endif
# else
            mu=fxa/fxb
            fxc=(c1-mu)/(c1+mu)
            po(1,1,south)=fxc*pb(1,1)+(c1-fxc)*p(2,2)
# endif
          else
            po(1,1,south)=p(2,2)
          endif 
# ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
           po(1,1,south)=p(2,2)
# endif
        else
           po(1,1,south)=pb(1,1)
        endif
        po(1,1,west)=po(1,1,south)
c
c  SE corner
c
        fxa=(pb(imtm1,2)-p(imtm1,2))
        fxb=p(imtm1,2)+pb(imtm1,2)-c2*p(imtm2,3)
        if(fxa*fxb.gt.c0) then
          if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
            if(abs(fxa).lt.small*abs(fxb)) then
              po(imt,1,south)=pb(imt,1)
            else
              mu=fxa/fxb
              po(imt,1,south)=((c1-mu)*pb(imt,1)+c2*mu*p(imtm1,2))
     &                        /(c1+mu)
            endif
# else
            mu=fxa/fxb
            fxc=(c1-mu)/(c1+mu)
            po(imt,1,south)=fxc*pb(imt,1)+(c1-fxc)*p(imtm1,2)
# endif
          else
            po(imt,1,south)=p(imtm1,2)
          endif 
# ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
           po(imt,1,south)=p(imtm1,2)
# endif
        else
           po(imt,1,south)=pb(imt,1)
        endif
        po(1,1,east)=po(imt,1,south)
c
#endif
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
#ifndef pjhcrnr
        do 1210 i=1,imt
#else
        do 1210 i=2,imtm1
#endif
          fxa=(pb(i,jmtm1)-p(i,jmtm1))
          fxb=p(i,jmtm1)+pb(i,jmtm1)-c2*p(i,jmtm2)
          if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
              if(abs(fxa).lt.small*abs(fxb)) then
                po(i,1,north)=pb(i,jmt)
              else
                mu=fxa/fxb
                po(i,1,north)=((c1-mu)*pb(i,jmt)+c2*mu*p(i,jmtm1))/
     *                        (c1+mu)
              endif
#else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              po(i,1,north)=fxc*pb(i,jmt)+(c1-fxc)*p(i,jmtm1)
#endif
            else
              po(i,1,north)=p(i,jmtm1)
            endif 
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
             po(i,1,north)=p(i,jmtm1)
#endif
          else
             po(i,1,north)=pb(i,jmt)
          endif
1210    continue
c
#ifndef pjhcrnr
        po(1,1,west)=p5*( po(2,1,south)+ po(2,1,west))
        po(1,1,south)=po(1,1,west)
c
        po(1,1,east)=p5*( po(imtm1,1,south)+ po(2,1,east) )
        po(imt,1,south)=po(1,1,east)
c
        po(jmt,1,west)=p5*( po(2,1,north)+ po(jmtm1,1,west) )
        po(1,1,north)=po(jmt,1,west)
c
        po(jmt,1,east)=p5*( po(imtm1,1,north)+ po(jmtm1,1,east) )
        po(imt,1,north)= po(jmt,1,east)
#else
c
c  NW corner
c
        fxa=(pb(2,jmtm1)-p(2,jmtm1))
        fxb=p(2,jmtm1)+pb(2,jmtm1)-c2*p(3,jmtm2)
        if(fxa*fxb.gt.c0) then
          if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
            if(abs(fxa).lt.small*abs(fxb)) then
              po(1,1,north)=pb(1,jmt)
            else
              mu=fxa/fxb
              po(1,1,north)=((c1-mu)*pb(1,jmt)+c2*mu*p(2,jmtm1))/
     *                      (c1+mu)
            endif
# else
            mu=fxa/fxb
            fxc=(c1-mu)/(c1+mu)
            po(1,1,north)=fxc*pb(1,jmt)+(c1-fxc)*p(2,jmtm1)
# endif
          else
            po(1,1,north)=p(2,jmtm1)
          endif 
# ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
           po(1,1,north)=p(2,jmtm1)
# endif
        else
           po(1,1,north)=pb(1,jmt)
        endif
        po(jmt,1,west)=po(1,1,north)
c
c  NE corner
c
        fxa=(pb(imtm1,jmtm1)-p(imtm1,jmtm1))
        fxb=p(imtm1,jmtm1)+pb(imtm1,jmtm1)-c2*p(imtm2,jmtm2)
        if(fxa*fxb.gt.c0) then
          if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
            if(abs(fxa).lt.small*abs(fxb)) then
              po(imt,1,north)=pb(imt,jmt)
            else
              mu=fxa/fxb
              po(imt,1,north)=((c1-mu)*pb(imt,jmt)+c2*mu*p(imtm1,jmtm1))
     *                        /(c1+mu)
            endif
# else
            mu=fxa/fxb
            fxc=(c1-mu)/(c1+mu)
            po(imt,1,north)=fxc*pb(imt,jmt)+(c1-fxc)*p(imtm1,jmtm1)
# endif
          else
            po(imt,1,north)=p(imtm1,jmtm1)
          endif 
# ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
           po(imt,1,north)=p(imtm1,jmtm1)
# endif
        else
           po(imt,1,north)=pb(imt,jmt)
        endif
        po(jmt,1,east)= po(imt,1,north)
#endif
# ifdef coast
c
        call ind_bdy(ind_b)
#   if defined islands
        do 1211 l=1,ncseg-nisle
#   else
        do 1211 l=1,ncseg
#   endif
           po(ind_b(2,1,l),1,ind_b(2,2,l)) = (c1/c3)*
     &                  ( po(ind_b(1,1,l),1,ind_b(1,2,l)) +
     &                    po(ind_b(2,1,l),1,ind_b(2,2,l)) +
     &                    po(ind_b(4,1,l),1,ind_b(4,2,l)) )
           po(ind_b(3,1,l),1,ind_b(3,2,l)) =
     &                    po(ind_b(2,1,l),1,ind_b(2,2,l))
 1211   continue
# endif
      endif
      return
c
c=======================================================================
c  Set Modified Orlanski Implicit (MOI) radiation open boundary
c  conditions on transport streamfunction.
c=======================================================================
c
      entry probc_moi (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        ic=1
        fxa=pb(ib+ic,j)-p (ib+ic,j)
        fxb=p (ib+ic,j)+pb(ib+ic,j)-c2*p(ib+2*ic,j)
        if(fxa*fxb.gt.c0) then
           po(j,1,west)=p(ib+ic,j)
        else
           po(j,1,west)=pb(ib,j)
        endif
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt
        ic=-1
        fxa=pb(ib+ic,j)-p (ib+ic,j)
        fxb=p (ib+ic,j)+pb(ib+ic,j)-c2*p(ib+2*ic,j)
        if(fxa*fxb.gt.c0) then
           po(j,1,east)=p(ib+ic,j)
        else
           po(j,1,east)=pb(ib,j)
        endif
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        do 1300 i=1,imt
          fxa=pb(i,2)-p (i,2)
          fxb=p (i,2)+pb(i,2)-c2*p(i,3)
          if(fxa*fxb.gt.c0) then
             po(i,1,south)=p(i,2)
          else
             po(i,1,south)=pb(i,1)
          endif
1300    continue    
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        do 1310 i=1,imt
          fxa=pb(i,jmtm1)-p (i,jmtm1)
          fxb=p (i,jmtm1)+pb(i,jmtm1)-c2*p(i,jmtm2)
          if(fxa*fxb.gt.c0) then
             po(i,1,north)=p(i,jmtm1)
          else
             po(i,1,north)=pb(i,jmt)
          endif
1310    continue    
c
        po(1,1,west)=p5*( po(2,1,south)+ po(2,1,west))
        po(1,1,south)=po(1,1,west)
c
        po(1,1,east)=p5*( po(imtm1,1,south)+ po(2,1,east) )
        po(imt,1,south)=po(1,1,east)
c
        po(jmt,1,west)=p5*( po(2,1,north)+ po(jmtm1,1,west) )
        po(1,1,north)=po(jmt,1,west)
c
        po(jmt,1,east)=p5*( po(imtm1,1,north)+ po(jmtm1,1,east) )
        po(imt,1,north)= po(jmt,1,east)
# ifdef coast
        call ind_bdy(ind_b)
#   if defined islands
        do 1311 l=1,ncseg-nisle
#   else
        do 1311 l=1,ncseg
#   endif
           po(ind_b(2,1,l),1,ind_b(2,2,l)) = (c1/c3)*
     &                  ( po(ind_b(1,1,l),1,ind_b(1,2,l)) +
     &                    po(ind_b(2,1,l),1,ind_b(2,2,l)) +
     &                    po(ind_b(4,1,l),1,ind_b(4,2,l)) )
           po(ind_b(3,1,l),1,ind_b(3,2,l)) =
     &                    po(ind_b(2,1,l),1,ind_b(2,2,l))
 1311   continue
# endif
      endif
      return
c
c=======================================================================
c  Set Orlanski Radiation Implicit (ORI) open boundary conditions
c  on vorticity tendency.
c=======================================================================
c
      entry zrobc_ori (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=2
        ic=1

          fxa=(ztdb(ib+ic,j)-ztd(ib+ic,j))
          fxb=ztd(ib+ic,j)+ztdb(ib+ic,j)-c2*ztd(ib+2*ic,j)
          if(fxa*fxb.gt.c0) then
             if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
                if(abs(fxa).lt.small*abs(fxb)) then
                   po(j,2,west)=ztdb(ib,j)/dtsf
                else
                   mu=fxa/fxb
                   po(j,2,west)=
     *                  ((c1-mu)*ztdb(ib,j)+c2*mu*ztd(ib+ic,j))/
     *                  (c1+mu)/dtsf
                endif
#else
                mu=fxa/fxb
                fxc=(c1-mu)/(c1+mu)
                po(j,2,west)=(fxc*ztdb(ib,j)+(c1-fxc)*ztd(ib+ic,j))
     *                       /dtsf
#endif
             else
                po(j,2,west)=ztd(ib+ic,j)/dtsf
             endif
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
             po(j,2,west)=ztd(ib+ic,j)/dtsf
#endif
          else
             po(j,2,west)=ztdb(ib,j)/dtsf
          endif

        ztd(ib,j)=po(j,2,west)*dtsf
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt-1
        ic=-1
c
          fxa=(ztdb(ib+ic,j)-ztd(ib+ic,j))
          fxb=ztd(ib+ic,j)+ztdb(ib+ic,j)-c2*ztd(ib+2*ic,j)
          if(fxa*fxb.gt.c0) then
             if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
                if(abs(fxa).lt.small*abs(fxb)) then
                   po(j,2,east)=ztdb(ib,j)/dtsf
                else
                   mu=fxa/fxb
                   po(j,2,east)=
     *                  ((c1-mu)*ztdb(ib,j)+c2*mu*ztd(ib+ic,j))/
     *                  (c1+mu)/dtsf
                endif
#else
                mu=fxa/fxb
                fxc=(c1-mu)/(c1+mu)
                po(j,2,east)=(fxc*ztdb(ib,j)+(c1-fxc)*ztd(ib+ic,j))
     *                       /dtsf
#endif
             else
                po(j,2,east)=ztd(ib+ic,j)/dtsf
             endif
#ifdef pjhori
          elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
             po(j,2,east)=ztd(ib+ic,j)/dtsf
#endif
          else
            po(j,2,east)=ztdb(ib,j)/dtsf
          endif
        ztd(ib,j)=po(j,2,east)*dtsf
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
#ifndef pjhcrnr
        do 1460 i=2,imtm1
#else
        do 1460 i=3,imtm2
#endif
c
            fxa=(ztdb(i,3)-ztd(i,3))
            fxb=ztd(i,3)+ztdb(i,3)-c2*ztd(i,4)
            if(fxa*fxb.gt.c0) then
               if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
                  if(abs(fxa).lt.small*abs(fxb)) then
                     po(i,2,south)=ztdb(i,2)/dtsf
                  else
                     mu=fxa/fxb
                     po(i,2,south)=((c1-mu)*ztdb(i,2)+c2*mu*ztd(i,3))/
     *                    (c1+mu)/dtsf
                  endif
#else
                  mu=fxa/fxb
                  fxc=(c1-mu)/(c1+mu)
                  po(i,2,south)=(fxc*ztdb(i,2)+(c1-fxc)*ztd(i,3))
     *                          /dtsf
#endif
               else
                  po(i,2,south)=ztd(i,3)/dtsf
               endif
#ifdef pjhori
            elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
               po(i,2,south)=ztd(i,3)/dtsf
#endif
            else
               po(i,2,south)=ztdb(i,2)/dtsf
            endif
          ztd(i,2)=po(i,2,south)*dtsf
1460    continue    
#ifdef pjhcrnr
c
        fxa=(ztdb(3,3)-ztd(3,3))
        fxb=ztd(3,3)+ztdb(3,3)-c2*ztd(4,4)
        if(fxa*fxb.gt.c0) then
           if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small*abs(fxb)) then
                 po(2,2,south)=ztdb(2,2)/dtsf
              else
                 mu=fxa/fxb
                 po(2,2,south)=((c1-mu)*ztdb(2,2)+c2*mu*ztd(3,3))/
     *                         (c1+mu)/dtsf
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              po(2,2,south)=(fxc*ztdb(2,2)+(c1-fxc)*ztd(3,3))
     *                       /dtsf
# endif
           else
              po(2,2,south)=ztd(3,3)/dtsf
           endif
# ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
           po(2,2,south)=ztd(3,3)/dtsf
# endif
        else
           po(2,2,south)=ztdb(2,2)/dtsf
        endif
        ztd(2,2)=po(2,2,south)*dtsf
        po(2,2,west)=po(2,2,south)
c
        fxa=(ztdb(imtm2,3)-ztd(imtm2,3))
        fxb=ztd(imtm2,3)+ztdb(imtm2,3)-c2*ztd(imt-3,4)
        if(fxa*fxb.gt.c0) then
           if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
              if(abs(fxa).lt.small*abs(fxb)) then
                 po(imtm1,2,south)=ztdb(imtm1,2)/dtsf
              else
                 mu=fxa/fxb
                 po(imtm1,2,south)=((c1-mu)*ztdb(2,2)+
     *                              c2*mu*ztd(imtm2,3))/(c1+mu)/dtsf
              endif
# else
              mu=fxa/fxb
              fxc=(c1-mu)/(c1+mu)
              po(imtm1,2,south)=(fxc*ztdb(2,2)+(c1-fxc)*ztd(imtm2,3))
     *                          /dtsf
# endif
           else
              po(imtm1,2,south)=ztd(imtm2,3)/dtsf
           endif
# ifdef pjhori
        elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
           po(imtm1,2,south)=ztd(imtm2,3)/dtsf
# endif
        else
           po(imtm1,2,south)=ztdb(imtm1,2)/dtsf
        endif
        ztd(imtm1,2)=po(imtm1,2,south)*dtsf
        po(2,2,east)=po(imtm1,2,south)
#endif
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
#ifndef pjhcrnr
         do 1465 i=2,imtm1
#else
         do 1465 i=3,imtm2
#endif
               fxa=(ztdb(i,jmt-2)-ztd(i,jmt-2))
               fxb=ztd(i,jmt-2)+ztdb(i,jmt-2)-c2*ztd(i,jmt-3)
               if(fxa*fxb.gt.c0) then
                  if(abs(fxb).gt.abs(fxa)) then
#ifndef pjhori
                     if(abs(fxa).lt.small*abs(fxb)) then
                        po(i,2,north)=ztdb(i,jmt-1)/dtsf
                     else
                        mu=fxa/fxb
                        po(i,2,north)=((c1-mu)*ztdb(i,jmt-1)+
     *                          c2*mu* ztd (i,jmt-2))/
     *                        (c1+mu)/dtsf
                     endif
#else
                     mu=fxa/fxb
                     fxc=(c1-mu)/(c1+mu)
                     po(i,2,north)=(fxc*ztdb(i,jmt-1)+
     *                             (c1-fxc)*ztd(i,jmt-2))/dtsf
#endif
                  else
                     po(i,2,north)=ztd(i,jmt-2)/dtsf
                  endif
#ifdef pjhori
               elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
                  po(i,2,north)=ztd(i,jmt-2)/dtsf
#endif
               else
                  po(i,2,north)=ztdb(i,jmt-1)/dtsf
               endif
            ztd(i,jmtm1)=po(i,2,north)*dtsf
 1465    continue
#ifdef pjhcrnr
c
         fxa=(ztdb(3,jmtm2)-ztd(3,jmtm2))
         fxb=ztd(3,jmtm2)+ztdb(3,jmtm2)-c2*ztd(4,jmt-3)
         if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
               if(abs(fxa).lt.small*abs(fxb)) then
                  po(2,2,north)=ztdb(2,jmtm1)/dtsf
               else
                  mu=fxa/fxb
                  po(2,2,north)=((c1-mu)*ztdb(2,jmtm1)+
     *                              c2*mu* ztd (3,jmtm2))/
     *                              (c1+mu)/dtsf
               endif
# else
               mu=fxa/fxb
               fxc=(c1-mu)/(c1+mu)
               po(2,2,north)=(fxc*ztdb(2,jmtm1)+
     *                        (c1-fxc)*ztd(3,jmtm2))/dtsf
# endif
            else
               po(2,2,north)=ztd(3,jmtm2)/dtsf
            endif
# ifdef pjhori
         elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            po(2,2,north)=ztd(3,jmtm2)/dtsf
# endif
         else
            po(2,2,north)=ztdb(2,jmtm1)/dtsf
         endif
         ztd(2,jmtm1)=po(2,2,north)*dtsf
         po(jmtm1,2,west)=po(2,2,north)
c
         fxa=(ztdb(imtm2,jmtm2)-ztd(imtm2,jmtm2))
         fxb=ztd(imtm2,jmtm2)+ztdb(imtm2,jmtm2)-c2*ztd(imt-3,jmt-3)
         if(fxa*fxb.gt.c0) then
            if(abs(fxb).gt.abs(fxa)) then
# ifndef pjhori
               if(abs(fxa).lt.small*abs(fxb)) then
                  po(imtm1,2,north)=ztdb(imtm1,jmtm1)/dtsf
               else
                  mu=fxa/fxb
                  po(imtm1,2,north)=((c1-mu)*ztdb(imtm1,jmtm1)+
     *                              c2*mu* ztd (imtm2,jmtm2))/
     *                              (c1+mu)/dtsf
               endif
# else
               mu=fxa/fxb
               fxc=(c1-mu)/(c1+mu)
               po(imtm1,2,north)=(fxc*ztdb(imtm1,jmtm1)+
     *                            (c1-fxc)*ztd(imtm2,jmtm2))/dtsf
# endif
            else
               po(imtm1,2,north)=ztd(imtm2,jmtm2)/dtsf
            endif
# ifdef pjhori
         elseif ((fxb.eq.c0).and.(fxa.gt.c0)) then
            po(imtm1,2,north)=ztd(imtm2,jmtm2)/dtsf
# endif
         else
            po(imtm1,2,north)=ztdb(imtm1,jmtm1)/dtsf
         endif
         ztd(imtm1,jmtm1)=po(imtm1,2,north)*dtsf
         po(jmtm1,2,east)=po(imtm1,2,north)
#endif
      endif
      return
c
c=======================================================================
c  Set Spall and Robinson open boundary conditions on barotropic
c  vorticity (CFvN type).
c=======================================================================
c
      entry zrobc_sr (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        delbar=(po(j-1,1,west)-po(j+1,1,west))*
     *         hr(ib,j)*dxu2r(ib)*csr(j)
        if(j.eq.2) then
          delbar1=(p(ib+2,j  )-p(ib  ,j  ))*hr(ib+1,j-1)*dyt2r(j)
          delbar2=(p(ib+1,j-1)-p(ib+1,j+1))*hr(ib  ,j  )*dyt2r(j)
          if((delbar1.gt.c0).or.(delbar2.gt.c0)) delbar=c1
        elseif(j.eq.jmtm1) then
          delbar1=(p(ib+2,j  )-p(ib  ,j  ))*hr(ib+1,j+1)*dyt2r(j)
          delbar2=(p(ib+1,j-1)-p(ib+1,j+1))*hr(ib  ,j  )*dyt2r(j)
          if((delbar1.lt.c0).or.(delbar2.gt.c0)) delbar=c1
        endif
c
        if(delbar.lt.c0) then
          uvel=(p(ib+1,j-1)-p(ib+1,j+1))*dyt2r(j)*hr(ib+1,j)
          vvel=(p(ib+2,j  )-p(ib  ,j  ))*dyt2r(j)*hr(ib+1,j)
          call zetabar(ib+2,j  ,eastz)
          call zetabar(ib+1,j  ,centz)
          call zetabar(ib+1,j-1,soutz)
          call zetabar(ib+1,j+1,znrth)
          dzdx=(eastz-centz)*csr(j)*dxtr(ib+1)
          if(vvel.gt.c0) then
            dzdy=(centz-soutz)*dytr(j)
          else
            dzdy=(znrth-centz)*dytr(j)
          endif
          po(j,2,west)=-uvel*dzdx-vvel*dzdy
          ztd(ib+1,j)=c2dtsf*po(j,2,west)
        else 
          ztd(ib+1,j)=c2dtsf*po(j,2,west)
        endif
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imt
        delbar=(po(j-1,1,east)-po(j+1,1,east))*
     *         hr(ib,j)*dxu2r(ib-1)*csr(j)
        if(j.eq.2) then
          delbar1=(p(ib  ,j  )-p(ib-2,j  ))*hr(ib-1,j-1)*dyt2r(j)
          delbar2=(p(ib-1,j-1)-p(ib-1,j+1))*hr(ib  ,j  )*dyt2r(j)
          if((delbar1.gt.c0).or.(delbar2.lt.c0)) delbar=cm1
        elseif(j.eq.jmtm1) then
          delbar1=(p(ib,j)-p(ib-2,j))*hr(ib-1,j+1)*dyt2r(j)
          delbar2=(p(ib,j)-p(ib-2,j))*hr(ib  ,j  )*dyt2r(j)
          if((delbar1.lt.c0).or.(delbar2.lt.c0)) delbar=cm1
        endif
c
        if(delbar.gt.c0) then
          uvel=(p(ib-1,j-1)-p(ib-1,j+1))*dyt2r(j)*hr(ib-1,j)
          vvel=(p(ib  ,j  )-p(ib-2,j  ))*dyt2r(j)*hr(ib-1,j)
          call zetabar(ib-1,j  ,centz)
          call zetabar(ib-2,j  ,westz)
          call zetabar(ib-1,j-1,soutz)
          call zetabar(ib-1,j+1,znrth)
          dzdx=(centz-westz)*csr(j)*dxtr(imt-1)
          if(vvel.gt.c0) then
            dzdy=(centz-soutz)*dytr(j)
          else
            dzdy=(znrth-centz)*dytr(j)
          endif
          po(j,2,east)=-uvel*dzdx-vvel*dzdy
          ztd(ib-1,j)=c2dtsf*po(j,2,east)
        else
          ztd(ib-1,j)=c2dtsf*po(j,2,east)
        endif
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        jb=1
        do 1400 i=2,imtm1
          delbar=(po(i+1,1,south)-po(i-1,1,south))*
     *           hr(i,jb)*dxu2r(i)*csr(jb)
          if(delbar.lt.c0) then
            uvel=(p(i  ,jb  )-p(i  ,jb+2))*dyt2r(jb+1)*hr(i,jb+1)
            vvel=(p(i+1,jb+1)-p(i-1,jb+1))*dyt2r(jb+1)*hr(i,jb+1)
            call zetabar(i  ,jb+1,centz)
            call zetabar(i-1,jb+1,westz)
            call zetabar(i+1,jb+1,eastz)
            call zetabar(i  ,jb+2,znrth)
            if(uvel.gt.c0) then
              dzdx=(centz-westz)*csr(jb+1)*dxtr(i)
            else
              dzdx=(eastz-centz)*csr(jb+1)*dxtr(i)
            endif
            dzdy=(znrth-centz)*dytr(jb+1)
            po(i,2,south)=-uvel*dzdx-vvel*dzdy
            ztd(i,jb+1)=c2dtsf*po(i,2,south)
          else
            ztd(i,jb+1)=c2dtsf*po(i,2,south)
          endif
1400    continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        jb=jmt
        do 1410 i=2,imtm1
          delbar=(po(i+1,1,north)-po(i-1,1,north))*
     *           hr(i,jb)*dxu2r(i)*csr(jb-1)
          vvel=(p(i+1,jb-1)-p(i-1,jb-1))*dyt2r(jb-1)*hr(i,jb-1)
          if(vvel.gt.c0) then
            uvel=(p(i  ,jb-2)-p(i  ,jb  ))*dyt2r(jb-1)*hr(i,jb-1)
            vvel=(p(i+1,jb-1)-p(i-1,jb-1))*dyt2r(jb-1)*hr(i,jb-1)
            call zetabar(i  ,jb-1,centz)
            call zetabar(i-1,jb-1,westz)
            call zetabar(i+1,jb-1,eastz)
            call zetabar(i  ,jb-2,soutz)
            if(uvel.gt.c0) then
              dzdx=(centz-westz)*csr(jb-1)*dxtr(i)
            else
              dzdx=(eastz-centz)*csr(jb-1)*dxtr(i)
            endif
            dzdy=(centz-soutz)*dytr(jb-1)
            po(i,2,north)=-uvel*dzdx-vvel*dzdy
            ztd(i,jb-1)=c2dtsf*po(i,2,north)
          else
            ztd(i,jb-1)=c2dtsf*po(i,2,north)
          endif
1410    continue
      endif
      return
c
c=======================================================================
c  Set time-extrapolated boundary conditions on the rate of change of
c  barotropic vorticity at the first interior point.
c=======================================================================
c
      entry zrobc_ext (j,ibdy)
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=2
        fxa=c2*cstr(j)*cstr(j)
        fxb=c2*cs  (j)*cstr(j)*dytr(j)*dyur(j  )
        fxc=c2*cs(j-1)*cstr(j)*dytr(j)*dyur(j-1)
        zfn=fxb/(hdv(ib-1,j  )+hdv(ib,j  ))
        zfs=fxc/(hdv(ib-1,j-1)+hdv(ib,j-1))
        zfe=fxa*dxur(ib  )*dxtr(ib)/(hdv(ib  ,j)+hdv(ib  ,j-1))
        zfw=fxa*dxur(ib-1)*dxtr(ib)/(hdv(ib-1,j)+hdv(ib-1,j-1))
        zpn=p(ib  ,j+1)-pb(ib  ,j+1)
        zps=p(ib  ,j-1)-pb(ib  ,j-1)
        zpe=p(ib+1,j  )-pb(ib+1,j  )
        zpw=p(ib-1,j  )-pb(ib-1,j  )
        zpc=p(ib  ,j  )-pb(ib  ,j  )
        ztd(ib,j)=zfn*zpn+zfs*zps+zfe*zpe+zfw*zpw-(zfn+zfs+zfe+zfw)*zpc
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imtm1
        fxa=c2*cstr(j)*cstr(j)
        fxb=c2*cs  (j)*cstr(j)*dytr(j)*dyur(j  )
        fxc=c2*cs(j-1)*cstr(j)*dytr(j)*dyur(j-1)
        zfn=fxb/(hdv(ib-1,j  )+hdv(ib,j  ))
        zfs=fxc/(hdv(ib-1,j-1)+hdv(ib,j-1))
        zfe=fxa*dxur(ib  )*dxtr(ib)/(hdv(ib  ,j)+hdv(ib  ,j-1))
        zfw=fxa*dxur(ib-1)*dxtr(ib)/(hdv(ib-1,j)+hdv(ib-1,j-1))
        zpn=p(ib  ,j+1)-pb(ib  ,j+1)
        zps=p(ib  ,j-1)-pb(ib  ,j-1)
        zpe=p(ib+1,j  )-pb(ib+1,j  )
        zpw=p(ib-1,j  )-pb(ib-1,j  )
        zpc=p(ib  ,j  )-pb(ib  ,j  )
        ztd(ib,j)=zfn*zpn+zfs*zps+zfe*zpe+zfw*zpw-(zfn+zfs+zfe+zfw)*zpc
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        jb=2
        do 1500 i=2,imtm1
          fxa=c2*cstr(jb)*cstr(jb)
          fxb=c2*cs  (jb)*cstr(jb)*dytr(jb)*dyur(jb  )
          fxc=c2*cs(jb-1)*cstr(jb)*dytr(jb)*dyur(jb-1)
          zfn=fxb/(hdv(i-1,jb  )+hdv(i,jb  ))
          zfs=fxc/(hdv(i-1,jb-1)+hdv(i,jb-1))
          zfe=fxa*dxur(i  )*dxtr(i)/(hdv(i  ,jb)+hdv(i  ,jb-1))
          zfw=fxa*dxur(i-1)*dxtr(i)/(hdv(i-1,jb)+hdv(i-1,jb-1))
          zpn=p(i  ,jb+1)-pb(i  ,jb+1)
          zps=p(i  ,jb-1)-pb(i  ,jb-1)
          zpe=p(i+1,jb  )-pb(i+1,jb  )
          zpw=p(i-1,jb  )-pb(i-1,jb  )
          zpc=p(i  ,jb  )-pb(i  ,jb  )
          ztd(i,jb)=zfn*zpn+zfs*zps+zfe*zpe+zfw*zpw-
     *              (zfn+zfs+zfe+zfw)*zpc
1500    continue
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        jb=jmtm1
        do 1510 i=2,imtm1
          fxa=c2*cstr(jb)*cstr(jb)
          fxb=c2*cs  (jb)*cstr(jb)*dytr(jb)*dyur(jb  )
          fxc=c2*cs(jb-1)*cstr(jb)*dytr(jb)*dyur(jb-1)
          zfn=fxb/(hdv(i-1,jb  )+hdv(i,jb  ))
          zfs=fxc/(hdv(i-1,jb-1)+hdv(i,jb-1))
          zfe=fxa*dxur(i  )*dxtr(i)/(hdv(i  ,jb)+hdv(i  ,jb-1))
          zfw=fxa*dxur(i-1)*dxtr(i)/(hdv(i-1,jb)+hdv(i-1,jb-1))
          zpn=p(i  ,jb+1)-pb(i  ,jb+1)
          zps=p(i  ,jb-1)-pb(i  ,jb-1)
          zpe=p(i+1,jb  )-pb(i+1,jb  )
          zpw=p(i-1,jb  )-pb(i-1,jb  )
          zpc=p(i  ,jb  )-pb(i  ,jb  )
          ztd(i,jb)=zfn*zpn+zfs*zps+zfe*zpe+zfw*zpw-
     *              (zfn+zfs+zfe+zfw)*zpc
1510    continue
      endif
      return
c
c=======================================================================
c  Set reduced physics boundary conditions on the rate of change of
c  barotropic vorticity at the first interior point.  This variant
c  uses a half timestep (so the barotropic component is available) and
c  multiplies the difference by 2.
c=======================================================================
c
      entry zrobc_rph (j,ibdy)
c
      mu = c2dtts/dtts
c
      if(ibdy.eq.west) then
c
c-----------------------------------------------------------------------
c  Western boundary.
c-----------------------------------------------------------------------
c
        ib=1
        zunz(j,1,west)=c0
        zvnz(j,1,west)=c0
        do 1600 k=1,km
          fxa=u(ib,k)-ub(ib,k)
          fxb=v(ib,k)-vb(ib,k)
          zunz(j,1,west)=zunz(j,1,west)+fxa*dzvqz(ib,k,0)
          zvnz(j,1,west)=zvnz(j,1,west)+fxb*dzvqz(ib,k,0)
1600    continue
        zunz(j,1,west)=zunz(j,1,west)*hvav(ib,j)*mu
        zvnz(j,1,west)=zvnz(j,1,west)*hvav(ib,j)*mu
c
c  If J=2, bootstrap the computation of ZUSZ and ZVSZ because the 
c  southern boundary conditions have not been computed yet.
c
        if(j.eq.2) then
          zusz(j,1,west)=zunz(j,1,west)
          zvsz(j,1,west)=zvnz(j,1,west)
        endif
c
c  Compute rate of change of barotropic vorticity at the boundary.
c
        ztd(ib+1,j)=(zunz(j,2,west)*dxu(ib+1)+
     *               zunz(j,1,west)*dxu(ib  ))*cs(j  )-
     *              (zusz(j,2,west)*dxu(ib+1)+
     *               zusz(j,1,west)*dxu(ib  ))*cs(j-1)
        ztd(ib+1,j)=(((zvnz(j,2,west)-zvnz(j,1,west))*dyu(j  )+
     *                (zvsz(j,2,west)-zvsz(j,1,west))*dyu(j-1)-
     *                ztd(ib+1,j))*dxt2r(ib+1)*dytr(j))*cstr(j)
c
c  Tranfer quantites computed to the north of the present row
c  to be defined to the south in the computation of the next row.
c
        zusz(j+1,1,west)=zunz(j,1,west)
        zvsz(j+1,1,west)=zvnz(j,1,west)
c
      elseif(ibdy.eq.east) then
c
c-----------------------------------------------------------------------
c  Eastern boundary.
c-----------------------------------------------------------------------
c
        ib=imu
        zunz(j,1,east)=c0
        zvnz(j,1,east)=c0
        do 1610 k = 1,km
          fxa=u(ib,k)-ub(ib,k)
          fxb=v(ib,k)-vb(ib,k)
          zunz(j,1,east)=zunz(j,1,east)+fxa*dzvqz(ib,k,0)
          zvnz(j,1,east)=zvnz(j,1,east)+fxb*dzvqz(ib,k,0)       
1610    continue
        zunz(j,1,east)=zunz(j,1,east)*hvav(ib,j)*mu
        zvnz(j,1,east)=zvnz(j,1,east)*hvav(ib,j)*mu
c
c  If J=2, bootstrap the computation of ZUSZ and ZVSZ because the 
c  southern boundary conditions have not been computed yet.
c
        if(j.eq.2) then
           zusz(j,1,east)=zunz(j,1,east)
           zvsz(j,1,east)=zvnz(j,1,east)
        endif
c
c  Compute rate of change of barotropic vorticity at the boundary.
c
        ztd(ib,j) = (zunz(j,1,east)*dxu(ib  )+
     *               zunz(j,2,east)*dxu(ib-1))*cs(j  )-
     *              (zusz(j,1,east)*dxu(ib  )+
     *               zusz(j,2,east)*dxu(ib-1))*cs(j-1)
        ztd(ib,j) = (((zvnz(j,1,east)-zvnz(j,2,east))*dyu(j  )+
     *                (zvsz(j,1,east)-zvsz(j,2,east))*dyu(j-1)-
     *                ztd(ib,j))*dxt2r(ib)*dytr(j))*cstr(j)
c
c  Tranfer quantities computed to the north of the present row
c  to be defined to the south in the computation of the next row.
c
        zusz(j+1,1,east)=zunz(j,1,east)
        zvsz(j+1,1,east)=zvnz(j,1,east)
c
      elseif(ibdy.eq.south) then
c
c-----------------------------------------------------------------------
c  Southern boundary.
c-----------------------------------------------------------------------
c
        jb=2
        do 1630 i=2,imum1
          zusz(i,1,south)=c0
          zvsz(i,1,south)=c0
          do 1620 k = 1,km
            fxa=um(i,k)-ubm(i,k)
            fxb=vm(i,k)-vbm(i,k)
            zusz(i,1,south)=zusz(i,1,south)+fxa*dzvqz(i,k,1)
            zvsz(i,1,south)=zvsz(i,1,south)+fxb*dzvqz(i,k,1)       
1620      continue
          zusz(i,1,south)=zusz(i,1,south)*hvav(i,jb-1)*mu
          zvsz(i,1,south)=zvsz(i,1,south)*hvav(i,jb-1)*mu
c
c  Compute rate of change of barotropic vorticity at the boundary.
c
          ztd(i,jb)=(zunz(i  ,2,south)*dxu(i  )+
     *               zunz(i-1,2,south)*dxu(i-1))*cs(jb  )-
     *              (zusz(i  ,1,south)*dxu(i  )+
     *               zusz(i-1,1,south)*dxu(i-1))*cs(jb-1)
          ztd(i,jb)=(((zvnz(i,2,south)-zvnz(i-1,2,south))*dyu(jb  )+
     *                (zvsz(i,1,south)-zvsz(i-1,1,south))*dyu(jb-1)-
     *                ztd(i,jb))*dxt2r(i)*dytr(jb))*cstr(jb)
1630    continue
#ifndef cyclic
c
c  SW corner
c
        i = 1
        zusz(i,1,south)=c0
        zvsz(i,1,south)=c0
        do 1640 k = 1,km
            fxa=um(i,k)-ubm(i,k)
            fxb=vm(i,k)-vbm(i,k)
            zusz(i,1,south)=zusz(i,1,south)+fxa*dzvqz(i,k,1)
            zvsz(i,1,south)=zvsz(i,1,south)+fxb*dzvqz(i,k,1)       
1640      continue
        zusz(i,1,south)=zusz(i,1,south)*hvav(i,jb-1)*mu
        zvsz(i,1,south)=zvsz(i,1,south)*hvav(i,jb-1)*mu
c
        i = 2
        ztd(i,jb) = (zunz(i  ,2,south)*dxu(i  )+
     *               zunz(jb,1,west)*dxu(i-1))*cs(jb  )-
     *              (zusz(i  ,1,south)*dxu(i  )+
     *               zusz(i-1,1,south)*dxu(i-1))*cs(jb-1)
        ztd(i,jb) = (((zvnz(i,2,south)-zvnz(jb,1,west))*dyu(jb  )+
     *                (zvsz(i,1,south)-zvsz(i-1,1,south))*dyu(jb-1)-
     *                ztd(i,jb))*dxt2r(i)*dytr(jb))*cstr(jb)
c
c  SE corner
c
        i = imu
        zusz(i,1,south)=c0
        zvsz(i,1,south)=c0
        do 1650 k = 1,km
            fxa=um(i,k)-ubm(i,k)
            fxb=vm(i,k)-vbm(i,k)
            zusz(i,1,south)=zusz(i,1,south)+fxa*dzvqz(i,k,1)
            zvsz(i,1,south)=zvsz(i,1,south)+fxb*dzvqz(i,k,1)       
1650      continue
        zusz(i,1,south)=zusz(i,1,south)*hvav(i,jb-1)*mu
        zvsz(i,1,south)=zvsz(i,1,south)*hvav(i,jb-1)*mu
c
        ztd(i,jb) = (zunz(jb,1,east)*dxu(i  )+
     *               zunz(i-1,2,south)*dxu(i-1))*cs(jb  )-
     *              (zusz(i  ,1,south)*dxu(i  )+
     *               zusz(i-1,1,south)*dxu(i-1))*cs(jb-1)
        ztd(i,jb) = (((zvnz(jb,1,east)-zvnz(i-1,2,south))*dyu(jb  )+
     *                (zvsz(i,1,south)-zvsz(i-1,1,south))*dyu(jb-1)-
     *                ztd(i,jb))*dxt2r(i)*dytr(jb))*cstr(jb)
#endif
c
      elseif(ibdy.eq.north) then
c
c-----------------------------------------------------------------------
c  Northern boundary.
c-----------------------------------------------------------------------
c
        jb=jmtm1
        do 1670 i=2,imum1
          zunz(i,1,north)=c0
          zvnz(i,1,north)=c0
          do 1660 k = 1,km
            fxa=up(i,k)-ubp(i,k)
            fxb=vp(i,k)-vbp(i,k)
            zunz(i,1,north)=zunz(i,1,north)+fxa*dzvqz(i,k,0)
            zvnz(i,1,north)=zvnz(i,1,north)+fxb*dzvqz(i,k,0)       
1660      continue
          zunz(i,1,north)=zunz(i,1,north)*hvav(i,jb)*mu
          zvnz(i,1,north)=zvnz(i,1,north)*hvav(i,jb)*mu
c
c  Compute rate of change of barotropic vorticity at the boundary.
c
          ztd(i,jb)=(zunz(i  ,1,north)*dxu(i  )+
     *               zunz(i-1,1,north)*dxu(i-1))*cs(jb  )-
     *              (zunz(i  ,2,north)*dxu(i  )+
     *               zunz(i-1,2,north)*dxu(i-1))*cs(jb-1)
          ztd(i,jb)=(((zvnz(i,1,north)-zvnz(i-1,1,north))*dyu(jb  )+
     *                (zvnz(i,2,north)-zvnz(i-1,2,north))*dyu(jb-1)-
     *                ztd(i,jb))*dxt2r(i)*dytr(jb))*cstr(jb)
1670    continue
#ifndef cyclic
c
c  NW corner
c
        i = 1
        zunz(i,1,north)=c0
        zvnz(i,1,north)=c0
        do 1680 k = 1,km
          fxa=up(i,k)-ubp(i,k)
          fxb=vp(i,k)-vbp(i,k)
          zunz(i,1,north)=zunz(i,1,north)+fxa*dzvqz(i,k,0)
          zvnz(i,1,north)=zvnz(i,1,north)+fxb*dzvqz(i,k,0)       
1680    continue
        zunz(i,1,north)=zunz(i,1,north)*hvav(i,jb)*mu
        zvnz(i,1,north)=zvnz(i,1,north)*hvav(i,jb)*mu
c
        i = 2
        ztd(i,jb) = (zunz(i  ,1,north)*dxu(i  )+
     *               zunz(i-1,1,north)*dxu(i-1))*cs(jb  )-
     *              (zunz(i  ,2,north)*dxu(i  )+
     *               zunz(jb-1,1,west)*dxu(i-1))*cs(jb-1)
        ztd(i,jb) = (((zvnz(i,1,north)-zvnz(i-1,1,north))*dyu(jb  )+
     *                (zvnz(i,2,north)-zvnz(jb-1,1,west))*dyu(jb-1)-
     *                ztd(i,jb))*dxt2r(i)*dytr(jb))*cstr(jb)
c
c  NE corner
c
        i = imu
        zunz(i,1,north)=c0
        zvnz(i,1,north)=c0
        do 1690 k = 1,km
          fxa=up(i,k)-ubp(i,k)
          fxb=vp(i,k)-vbp(i,k)
          zunz(i,1,north)=zunz(i,1,north)+fxa*dzvqz(i,k,0)
          zvnz(i,1,north)=zvnz(i,1,north)+fxb*dzvqz(i,k,0)       
1690    continue
        zunz(i,1,north)=zunz(i,1,north)*hvav(i,jb)*mu
        zvnz(i,1,north)=zvnz(i,1,north)*hvav(i,jb)*mu
c
        ztd(i,jb) = (zunz(i  ,1,north)*dxu(i  )+
     *               zunz(i-1,1,north)*dxu(i-1))*cs(jb  )-
     *              (zunz(jb-1,1,east)*dxu(i  )+
     *               zunz(i-1,2,north)*dxu(i-1))*cs(jb-1)
        ztd(i,jb) = (((zvnz(i,1,north)-zvnz(i-1,1,north))*dyu(jb  )+
     *                (zvnz(jb-1,1,east)-zvnz(i-1,2,north))*dyu(jb-1)-
     *                ztd(i,jb))*dxt2r(i)*dytr(jb))*cstr(jb)
#endif
      endif
      return
      end
