c---------------------------------------------------------
       subroutine river_init (npt,nx,ny,nz,nzi,xm,ym,iox,h)
c---------------------------------------------------------
       dimension xm(nx),ym(ny),iox(npt),h(npt,nz),nzi(npt)

c---------------------------------------------------------
c INPUT
c    xm           longitude, in degrees
c    ym           latitude , in degrees
c    nx           number of x points 
c    iox          compressed->uncompressed map

c OUTPUT - in common
c    n_river      number of rivers
c    x_river      river longitude in degrees
c    y_river      river latitude  in degrees
c    z_river      river depth in meters
c    i_river      compressed index of river input
c    t_river      monthly volume of river flow
c    v_river      volume fraction of river flow in each ocean layer
c    kmin_river   minimum vertical index of river input
c    kmax_river   maximum vertical index of river input
c---------------------------------------------------------
       include 'comm_new.h'
       parameter (idf_ri = 85)

       open (unit = idf_ri, file = fbriv(1:n_riv))
       n = 0

       do i= 1, max_river
          read(idf_ri,*,END=110) 
          read(idf_ri,*,END=110) y_river(n+1),x_river(n+1),z_river(n+1),
     *             (t_river(n+1,j),j=1,12)
c#          print*,n+1,x_river(n+1),y_river(n+1),z_river(n+1),t_river(n+1,6)
          if (y_river(n+1).gt.-99) n = n + 1
       enddo

  110  continue 

       close (idf_ri)

       n_river = n

       open (unit = idf_ri, file = finp(1:mlen(finp))//'.riv')
       write(idf_ri,*) 'river #       i,   j,   k           ipt'
       write(idf_ri,*) 

       do n = 1, n_river
          x = x_river(n)
          y = y_river(n)
          dmin = 180 + 360

          do i = 1, npt
             ixy = iox(i)
             ix = mod (ixy -1 ,nx) + 1
             iy = (ixy - ix)/nx + 1
             xdist = min(abs(xm(ix)-x),abs(xm(ix)-x-360),abs(xm(ix)-x+360))
             dist = xdist + abs(ym(iy)-y)
             if (dist.lt.dmin) then
                dmin = dist
                i_river(n) = i
             endif
          enddo

          if (dmin.gt.dmin_river) then
             print*,'suppressing river #',n,' at ',x,y
             do i=1,12
                t_river(n,i) = 0.
             enddo
          endif

          ipt = i_river(n)
          dep = 0.
          do k = 1, nzi(ipt)
             dep = dep + h(ipt,k)
             if (z_river(n).lt.dep) then
                k_river(n) = k
                goto 200
             endif
          enddo

          k_river(n) = nzi(ipt)

  200     continue 

          iz_river = 5    
          d_river  = -50    

          kr   = k_river(n)
          kmin = max(   1    , kr - iz_river)
          kmax = min(nzi(ipt), kr + iz_river)

          kmin_river(n) = kmin
          kmax_river(n) = kmax

          ik = 0
          dep = 0.
          do k = 1, nzi(ipt)
             dep = dep + h(ipt,k)
             if (k.ge.kmin) then
                v_river(n,k-kr) = exp( abs(dep - z_river(n)) / d_river)
             endif
             if (k.ge.kmax) goto 300
          enddo

  300     continue 
          v_tot = 0.
          do k = kmin-kr,kmax-kr
             v_tot = v_tot + v_river(n,k)
          enddo
          if (v_tot.ne.0) then
             do k = kmin-kr,kmax-kr
                v_river(n,k) = v_river(n,k)/v_tot
             enddo
          endif

          ixy = iox(ipt)
          ix = mod (ixy -1 ,nx) + 1
          iy = (ixy - ix)/nx + 1
          write(idf_ri,400) n,ix,iy,kmin_river(n),kmax_river(n),i_river(n)
  400     format(I8,4x,I4,','I4,','I4,8x,I6)

       enddo

       close (idf_ri)

       return
       end
             
       subroutine river_flux (npt,nz,nzi,tmonth,emx,emy,fsal,fsal_ri,flor)
       dimension nzi(npt), emx(npt), emy(npt)
       dimension fsal(npt,nz), fsal_ri(npt,nz)
c---------------------------------------------------------
c INPUT
c    emx          1/dx, in meters
c    emy          1/dy, in meters
c    fsal         salt flux, in psu*m/s

c OUTPUT
c    fsal         salt flux + river flux, in psu*m/s
c---------------------------------------------------------
       include 'comm_new.h'

c   note that t_river(*,it) is the monthly mean, and therefore is at 
c                  tmonth = it-.5   (e.g., january has it=1, tmonth=0.5)

       it2 = int(tmonth+1.5)
       tscl = (it2- .5) - tmonth
       tscp = 1.- tscl
       if (it2.eq.13) it2=1
       it1 = it2- 1
       if (it1.eq.0)  it1=12

       do nr = 1, n_river
          i = i_river(nr)
          kr = k_river(nr)
          do k = kmin_river(nr), kmax_river(nr)
c              river flux, in m^3/s
             t_riv= coef_river* (tscl*t_river(nr,it1)+ tscp*t_river(nr,it2))
             t_riv= t_riv * v_river(nr,k - kr)
             budget_ri = budget_ri - flor*t_riv
c              river flux into grid cell, in psu m/s
             t_riv = t_riv* emx(i)* emy(i)
             fsal_ri(i,k) = fsal_ri(i,k) - flor*t_riv
             fsal   (i,k) = fsal   (i,k) - t_riv* salt_fixed
          enddo
       enddo

       return
       end
