c     The followings [search for "cgomc27:rivers"] indicate [simple] 
c     changes to be made to pom to effect rivers as described in 
c
c     Oey, 1996, JPO, pp.145-175; [river description is on p.154]
c
c     It's most convenient to declare a new array WRIV(IM,JM), which 
c     =0 except at grid points (ir(nr),jr(nr)) [nr=#rivers] where 
c     WRIV(ir(n),jr(n)) = -discharge_of_nth_river (m^3/s) divided by
c                         [dx(ir(n),jr(n))*dy(ir(n),jr(n))]
c     i.e., WRIV = downward velocity due to rivers [note negative sign]
c
c     A sample code follows:
c
c     Specify river discharge as downward vertical velocity
      wriv(:,:)=0.0
      do n=1,nr                    !nr=total_# of rivers
      j=jr(n)                      !3<jr<jm-2 to be safe
      i=ir(n)                      !3<ir<im-2
      wriv(i,j)=-qriv(n)/art(i,j)  !qriv=discharge_of_nth_river
      enddo
c
c
c     An older version [when pom still had S=S-35 etc and also i simply 
c     used w(*,*,1) w/o bothering to declare a new WRIV] was distributed 
c     to a number of pom users [who requested] since my 1996 paper.
c
c                          l.oey (oct/29/2001)
c
c---------------------------------------------------------------------
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      .
      .
      .
      .

C********** BEGIN EXTERNAL MODE ***************************************
                      DO 8000 IEXT=1,ISPLIT
      DO 405 J=2,JM
      DO 405 I=2,IM
      FLUXUA(I,J)=.25E0*(D(I,J)+D(I-1,J))*(DY(I,J)+DY(I-1,J))*UA(I,J)
 405  FLUXVA(I,J)=.25E0*(D(I,J)+D(I,J-1))*(DX(I,J)+DX(I,J-1))*VA(I,J)
C
      DO 410 J=2,JMM1
      DO 410 I=2,IMM1
  410 ELF(I,J)=ELB(I,J)
cgomc27:rivers
c--**1    -DTE2*(FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J))
c--**2                    /ART(I,J)
     1   -DTE2*((FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J))
c--  2                    /ART(I,J) + W(I,J,1) )
     2                    /ART(I,J) + WRIV(I,J))
C
C
      CALL BCOND(1)
C

      .
      .
      .
      .




C----------------------------------------------------------------
C     COMPUTE TF AN SF USING UF AND VF AS TEMPORARY VARIABLES
C----------------------------------------------------------------
      IF(MODE.EQ.4) GO TO 360
cgomc27:rivers:
c--   CALL ADVT(TB,T,TCLIM,DTI2,UF)
c--   CALL ADVT(SB,S,SCLIM,DTI2,VF)
c     Note that TSURF at grid points where WRIV.ne.0 is = river's temperature
c     This is convenient if one knows the temp [say cold and fresh] of river
c     One may also set it to some climatological temperature at the river
c     mouth
      CALL ADVT(TB,T,TCLIM,DTI2,UF,TSURF)
c     Put fresh water S=TPS=0.0 at grid points where WRIV.ne.0 (i.e. rivers)
c     [can be modified to put non-zero salinity]
      TPS(:,:)=SSURF(:,:)*(1.-WRIV(:,:)/(WRIV(:,:)+1.E-28))
      CALL ADVT(SB,S,SCLIM,DTI2,VF,TPS)

      .
      .
      .
      .


cgomc27:rivers
c--   SUBROUTINE ADVT(FB,F,FCLIM,DTI2,FF)
      SUBROUTINE ADVT(FB,F,FCLIM,DTI2,FF,FAT)
C
C     THIS SUBROUTINE INTEGRATES CONSERVATIVE SCALAR EQUATIONS
C
      .
      .
      .
      .

C
C****** DO VERTICAL ADVECTION *************************************
cgomc27:
      DO 505 J=2,JMM1
      DO 505 I=2,IMM1
cgomc27:rivers
c505  FF(I,J,1)=-.5*(F(I,J,1)+F(I,J,2))*W(I,J,2)*ART(I,J)/DZ(1)
 505  FF(I,J,1)=ART(I,J)*(          FAT(I,J)      *WRIV(I,J)
     1                   -.5E0*(F(I,J,1)+F(I,J,2))*W(I,J,2))/DZ(1)
c
      DO 520 K=2,KBM1
      DO 520 J=2,JMM1
      DO 520 I=2,IMM1

      .
      .
      .
      .

      RETURN
      END
c
      .
      .
      .
      .

c---------------------------------------------------------------------
c---------------------------------------------------------------------
c---------------------------------------------------------------------