      subroutine add_topo(nlng,nlat,nlev,ulng,ulat,h,lmask)

      implicit none

      integer nlng,nlat,nlev
      double precision ulat(nlng,0:nlat),ulng(nlng,0:nlat)
      double precision h(nlng,0:nlat)
      integer ntlng,ntlat,lmask(nlng,0:nlat)
      parameter (ntlng=10800,ntlat=5400)
      real rtopo (ntlng,ntlat)

      double precision x(3,5),ullng(4),ullat(4),pi
      double precision real_lat,real_lng,hh,hland,hocean,xp,yp,zp
      double precision xv1,xv2,xpv,yv1,yv2,ypv,zv1,zv2,zpv
      integer i,j,k,i1,i2,im1,j1,l,kp1,ln,la
      integer lngmax,lngmin,latmax,latmin,nbland,nbpocean,nbpland
      logical incell, land
      character*80 fname

      pi = 3.141592653589793d0

 
c-----------------------------------------------------------------
c     reading topography...
c     The  integer values (depth at each point) are in meters, positive
c     above sealevel, neg. below.
c     The resolution is 2 min. x 2min.
c     The data is organized in full 180 deg. latitude circles
c     beginning at 90N (constant val.) and stepping south to 89 55.0S
c-----------------------------------------------------------------

      print *, " "
      print *, " reading topo data..."

c  dataset below archived at /project/oce/bryan/etopo_merge_bedmap.da
      fname = '/cgd/oce/yeager/obs/topo/etopo_merge_bedmap.da'
      open (55, file=fname, form='unformatted', status='old',
     & access='direct',action='read',recl=ntlat*ntlng*4)
      read (55,rec=1) rtopo
      close(55)
 
      print *, " "
      print *, "read rtopo. r(5000,2500) = ",rtopo(5000,2500) 
      print *, "calculating depth field..."
 
c----------------------------------------------------------
c     is there land at T-points?
c----------------------------------------------------------

      where(ulng.gt.pi) ulng=ulng-2.0*pi
 
      nbland = 0
      do j=1,nlat
         print *, " j = ",j
         do i=1,nlng
            if (i.eq.1) then
               im1 = nlng
            else
               im1 = i-1
            endif
 
c     what are the coordonates of the 4 edges of the cell?
            do k=1,4
               if (k.eq.1.or.k.eq.2) then
                  j1 = j-1
               else
                  j1 = j
               endif
               if (k.eq.1.or.k.eq.4) then
                  i1 =  im1
               else
                  i1 = i
               endif
               x (1,k) = cos(ulng(i1,j1))*cos(ulat(i1,j1))
               x (2,k) = sin(ulng(i1,j1))*cos(ulat(i1,j1))
               x (3,k) = sin(ulat(i1,j1))
               ullng (k) = ulng(i1,j1)
               ullat (k) = ulat(i1,j1)
          if (i.eq.81.and.j.eq.138) then
           write(6,*) 'ullng,ullat = ',ullng(k),ullat(k)
          endif
            enddo
 
 
c----------------------------------------------------------
c     what are the (long-lat) boundaries of the cell ?
c----------------------------------------------------------

            lngmax = (180.0*30+int(max(ulng(i,j),ulng(im1,j),
     &           ulng(i,j-1),ulng(im1,j-1))*180.*30./pi))+1
            lngmin = (180.0*30+int(min(ulng(i,j),ulng(im1,j),
     &           ulng(i,j-1),ulng(im1,j-1))*180.*30./pi))+1
            latmin = (90*30 - int(max(ulat(i,j),ulat(i,j-1),
     &           ulat(im1,j),ulat(im1,j-1))*30.*180./pi)) + 1
            latmax = (90*30 - int(min(ulat(i,j),ulat(i,j-1),
     &           ulat(im1,j),ulat(im1,j-1))*30.*180./pi)) + 1
            if (float(lngmax-lngmin)/float(ntlng).gt.0.25) then
               l = lngmin
               lngmin = lngmax
               lngmax = l+ntlng
            endif
         if (i.eq.81.and.j.eq.138) then
            write(6,*) 'lngmin,lngmax = ',lngmin,lngmax
            write(6,*) 'latmin,latmax = ',latmin,latmax
         endif
 
c----------------------------------------------------------
c     calculation of depth
c----------------------------------------------------------
 
            hland = 0.
            hocean = 0.
            incell = .true.
            land = .false.
            nbpocean =0
            nbpland =0
            do j1=latmin,latmax
               do i2=lngmin,lngmax
 
                  i1 = mod(i2,ntlng)
                  incell = .true.
                  real_lng = float(i1-1)*2.*pi/float(ntlng)-pi
                  real_lat = pi/2.-float(j1-1)*pi/float(ntlat)
                  xp = cos(real_lng)*cos(real_lat)
                  yp = sin(real_lng)*cos(real_lat)
                  zp = sin(real_lat)
 
                  do k=1,4
                     kp1 = mod(k,4)+1
                     xv1 = x(1,k) - xp
                     yv1 = x(2,k) - yp
                     zv1 = x(3,k) - zp
                     xv2 = x(1,kp1) - xp
                     yv2 = x(2,kp1) - yp
                     zv2 = x(3,kp1) - zp
                     xpv = yv1*zv2 - zv1*yv2
                     ypv =-xv1*zv2 + zv1*xv2
                     zpv = xv1*yv2 - yv1*xv2
                     incell = incell.and.(xpv*xp+ypv*yp+zpv*zp.ge.0)
                  enddo
 
                  if (incell) goto 20
                  goto 10
 20               continue
                  if (rtopo(i1,j1).ge.0.0) then
                     nbpland  = nbpland +1
                  else
                     hocean = hocean + rtopo(i1,j1)
                     nbpocean  = nbpocean +1
                  endif
 
 
 10               continue
 
               enddo
            enddo
         if (i.eq.81.and.j.eq.138) then
            write(6,*) 'nbpland = ',nbpland
            write(6,*) 'nbpocean = ',nbpocean
         endif
 
 
            if (nbpland.gt.(nbpocean)) then
               hh =  0.0
            else
               if (nbpocean+nbpland.eq.0) then
c                 print *,'nbpocean+nbpland.eq.0: ',i,j
                  hh=0.0
                  do k=1,4
                     ln = (-180.0*30-int(ullng(k)*180.*30./pi))+1
                     la = (90*12 - int(ullat(k)*12.*180./pi)) + 1
c                     print *,'     ln,la,ullng,ullat = ',ln, la,
c    &                        ullng(k),ullat(k)
                     hh = hh+ rtopo(ln,la)
c                    print *,'     hh = ',hh
                     nbpocean = 1
                  enddo
                  hh = hh/4.
c                 print *,i,j,tlng(i,j)*180./pi,tlat(i,j)*180./pi,hh
               else
C                 hh = (hocean+hland)/float(nbpland+nbpocean)
                  hh = (hocean)/float(nbpocean)
               endif
            endif
            if (hh.ge.0.0) then
               nbland = nbland +1
            endif
            h(i,j) = -hh*100.0
c	if there is any land in this grid cell, don't smooth later
            if (nbpland.gt.0) then
              lmask(i,j) = 1
            else
              lmask(i,j) = 0
            endif
         if (i.eq.81.and.j.eq.138) then
            write(6,*) 'h = ',h(i,j)
            write(6,*) 'lmask = ',lmask(i,j)
         endif
 
	 enddo
c        print *,'nb de points :',nbland,' pour j ',j
      enddo
 
      end
 
