subroutine wad_real
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Sets up for real topodata problem w/wad.            *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      include 'pom08.c'
C
      real delx,elejmid,elwjmid,ra,vel,hland,dely  !lyo:!wad:define hland
      integer nvar                                 !lyo:!wad:define nvar
      integer i,j,k,m
      real rad,re,dlat,dlon,cff

      rad=0.01745329
      re=6371.E3

!lyo:!wad:Define variable or uniform grid option:
      nvar=1  !=1 for variable grid; =0 otherwise
!lyo:!wad:Special test case for smaller delx=4km:
C
C     Current velocity:
C
      vel=0.0e0
C
C     Read Ocean Topography dataset
C     
      Open(40,file='POM08_Topo.dat',form='formatted')
      do j=jm,1,-1
      do i=1,im
        read(40,*) east_e(i,j),north_e(i,j),h(i,j)
      enddo
      enddo
      close(40)

!
!lyo:_20080415:Move "calc. Coriolis etc" through "call areas_masks" 
!     from below (after call dens) to here; otherwise, fsm will not 
!     be defined but is used in subroutine dens. (This is a pom2k_bug).
C
C --- calc. Coriolis Parameter
C
        do j=1,jm
          do i=1,im
            cor(i,j)=2.*7.29E-5*sin(north_e(i,j)*rad)
!           aam2d(i,j)=aam(i,j,1) !lyo:_20080415:pom2k_bug:aam not defined
            elb(i,j)=0.
            etb(i,j)=0.
            dt(i,j)=h(i,j)
          end do
        end do
C
        do j=1,jm
          do i=2,im-1
            dx(i,j)=0.5*rad*re*sqrt(((east_e(i+1,j)-east_e(i-1,j))
     1 *cos(north_e(i,j)*rad))**2+(north_e(i+1,j)-north_e(i-1,j))**2)
          end do
            dx(1,j)=dx(2,j)
            dx(im,j)=dx(im-1,j)
        end do

        do i=1,im
          do j=2,jm-1
            dy(i,j)=0.5*rad*re*sqrt(((east_e(i,j+1)-east_e(i,j-1))
     1 *cos(north_e(i,j)*rad))**2+(north_e(i,j+1)-north_e(i,j-1))**2)
          end do
            dy(i,1)=dy(i,2)
            dy(i,jm)=dy(i,jm-1)
        end do

C
C --- the following grids are needed only for netcdf plotting
C
C     Corner of cell points:
C
      do j=2,jm
        do i=2,im
          east_c(i,j)=(east_e(i,j)+east_e(i-1,j)
     $                  +east_e(i,j-1)+east_e(i-1,j-1))/4.e0
          north_c(i,j)=(north_e(i,j)+north_e(i-1,j)
     $                   +north_e(i,j-1)+north_e(i-1,j-1))/4.e0
        end do
      end do
C
C
C     Extrapolate ends (approx.):
C
      do i=2,im
        east_c(i,1)=2.*east_c(i,2)-east_c(i,3)
        north_c(i,1)=2.*north_c(i,2)-north_c(i,3)
      end do
        east_c(1,1)=2.*east_c(2,1)-east_c(3,1)
C
      do j=2,jm
        east_c(1,j)=2.*east_c(2,j)-east_c(3,j)
        north_c(1,j)=2.*north_c(2,j)-north_c(3,j)
      end do
        north_c(1,1)=2.*north_c(1,2)-north_c(1,3)
C
C     u-points:
C
      do j=1,jm-1
        do i=1,im
          east_u(i,j)=(east_c(i,j)+east_c(i,j+1))/2.e0
          north_u(i,j)=(north_c(i,j)+north_c(i,j+1))/2.e0
        end do
      end do
C
C     Extrapolate ends:
C
      do i=1,im
        east_u(i,jm)=(east_c(i,jm)*3.e0-east_c(i,jm-1))/2.e0
        north_u(i,jm)=(north_c(i,jm)*3.e0-north_c(i,jm-1))/2.e0
      end do
C
C     v-points:
C
      do j=1,jm
        do i=1,im-1
          east_v(i,j)=(east_c(i,j)+east_c(i+1,j))/2.e0
          north_v(i,j)=(north_c(i,j)+north_c(i+1,j))/2.e0
        end do
      end do
C
C     Extrapolate ends:
C
      do j=1,jm
        east_v(im,j)=(east_c(im,j)*3.e0-east_c(im-1,j))/2.e0
        north_v(im,j)=(north_c(im,j)*3.e0-north_c(im-1,j))/2.e0
      end do
C
C     rot is the angle (radians, anticlockwise) of the i-axis relative
C     to east, averaged to a cell centre:
C
C     (NOTE that the following calculation of rot only works properly
C     for this particular rectangular grid)
C
      do j=1,jm
        do i=1,im
          rot(i,j)=0.e0
        end do
      end do
C
C     Define LANDSEA COVER
C
!tne:!wad:!lyo:!wad:
!     Note that unlike original seamount, h<0 is now water below MSL, and 
!     h>0 is above the MSL and is either land (if h>hland, see below)
!     or is potential wet-and-dry region
!
      hland=5.e0 !note all pnts are potential WAD if we set
!                             !hland >= -hmax*(1.e0-delh) (=7.5m);
!                             !i.e. "if" below is NOT satisfied
      do j=1,jm
      do i=1,im
        if(h(i,j).gt.hland) h(i,j) = hhi + 1.e0
      enddo
      enddo
C
C     Calculate areas and masks:
C
!lyo:!wad:Define "h" consistent with WAD and calculate wetmask, cell areas
!         and fsm etc.
      call wadh
C
C     Adjust bottom topography so that cell to cell variations
C     in h do not exceed parameter slmax:
C
!     if(slmax.lt.1.e0) call slpmax  !lyo:wad:now done in wadh above
C
C     Set initial conditions:
C
      do k=1,kbm1
        do j=1,jm
          do i=1,im
!lyo:wad:
             tb(i,j,k)=5.e0+15.e0*exp(zz(k)*h(i,j)/1000.e0)-tbias
             if (h(i,j).gt.hhi) then !lyo:wad:stratified water below MSL:
                tb(i,j,k)=5.e0+15.e0*exp(zz(k)*(h(i,j)-hhi)/1000.e0)
             else                    !lyo:wad:well-mixed water above MSL:
                tb(i,j,k)=5.e0+15.e0
             endif
            tb(i,j,k)=tb(i,j,k)-tbias 
!
            sb(i,j,k)=35.e0-sbias
            tclim(i,j,k)=tb(i,j,k)
            sclim(i,j,k)=sb(i,j,k)
            ub(i,j,k)=vel*dum(i,j)
          end do
        end do
      end do
C
C     Initialise uab and vab as necessary
C     (NOTE that these have already been initialised to zero in the
C     main program):
C
      do j=1,jm
        do i=1,im
          uab(i,j)=vel*dum(i,j)
        end do
      end do
C
C     Set surface boundary conditions, e_atmos, vflux, wusurf,
C     wvsurf, wtsurf, wssurf and swrad, as necessary
C     (NOTE:
C      1. These have all been initialised to zero in the main program.
C      2. The temperature and salinity of inflowing water must be
C         defined relative to tbias and sbias.):
C
      do j=1,jm
        do i=1,im
C     No conditions necessary for this problem
!
!lyo:!wad:!pom2k_bug:tsurf and ssurf were never defined, but should be:
             tsurf(i,j)=tb(i,j,1)
             ssurf(i,j)=sb(i,j,1)
!
        end do
      end do
C
C     Initialise elb, etb, dt and aam2d:
C
      do j=1,jm
        do i=1,im
          elb(i,j)=(-e_atmos(i,j)-hhi)*fsm(i,j)      !lyo:!wad:
!lyo:!wad:note:The following "if" is not satisfied if nwad=0 since     !
!     then fsm=wetmask at all cells                                    !
!     if (fsm(i,j).ne.0.0.and.wetmask(i,j).eq.0.0) ! dry but not land
      if (fsm(i,j).ne.wetmask(i,j))                ! dry but not land
     &elb(i,j)=-h(i,j)+hco !note slightly smaller hco<hc to ensure that
                           !initially, cell is wet with elb+h=hco < hc
          etb(i,j)=elb(i,j)               !lyo:!wad:
          dt(i,j)=h(i,j)+elb(i,j)         !lyo:!wad:
          aam2d(i,j)=aam(i,j,1)
        end do
      end do
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Make sure that initial wetmask remains compatible with       !
!     initial elb.  Also define marsh evaporation, +ve upward or       !
!     evaporate.  Also initialize wriv for rivers.                     !
!lyo:!wad:note:wriv is added here but the implementation following     !
!     Oey [1996; JPO, 26, 145-175; but see p.154 in particular] is not !
!     complete in this code.  E-mail me if want to implement it.       !
!                                                                      !
      if (nwad.eq.1) then
      do j=1,jm; do i=1,im
      wetmask(i,j)=fsm(i,j)
      if ((h(i,j)+elb(i,j)).le.hc) wetmask(i,j)=0.0
      enddo; enddo
      endif
!                                                                      !
      do j=1,jm; do i=1,im
      wmarsh(i,j)=0.0; wriv(i,j)=0.0
      if (fsm(i,j).ne.wetmask(i,j)) then           ! dry but not land
!     2 examples of wmarsh:
!--   wmarsh(i,j)=fsm(i,j)*10./art(i,j)
!--   wmarsh(i,j)=fsm(i,j)*(4.e-7)  !Gill's [1982] value
      endif
      enddo; enddo
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad: Set up pdens before 1st call dens; used also in profq:      !
      do k=1,kbm1; do j=1,jm; do i=1,im
         pdens(i,j,k)=grav*rhoref*(-zz(k)*max(h(i,j)-hhi,0.e0))*1.e-5
         enddo; enddo; enddo
!                                                                      !
      call dens(sb,tb,rho)
!                                                                      !
!----------------------------------------------------------------------!
C
C     Generated horizontally averaged density field (in this
C     application, the initial condition for density is a function
C     of z (the vertical cartesian coordinate) -- when this is not
C     so, make sure that rmean has been area averaged BEFORE transfer
C     to sigma coordinates):
C
      do k=1,kbm1
        do j=1,jm
          do i=1,im
            rmean(i,j,k)=rho(i,j,k)
          end do
        end do
      end do
C
C     Set lateral boundary conditions, for use in subroutine bcond
C     (in the seamount problem, the east and west boundaries are open,
C     while the south and north boundaries are closed through the
C     specification of the masks fsm, dum and dvm):
C
      rfe=1.e0
      rfw=1.e0
      rfn=1.e0
      rfs=1.e0
C
      do j=2,jmm1
        uabw(j)=uab(2,j)
        uabe(j)=uab(imm1,j)
C
C     Set geostrophically conditioned elevations at the boundaries:
C
!lyo:!wad:Note keep (temporarily) all el* defined wrt MSL - ele, elw,
!         eln & els were initialized to be =0 in MAIN; then adjust
!         them later (below) with hhi:
!
        ele(j)=ele(j-1)-cor(imm1,j)*uab(imm1,j)/grav*dy(imm1,j-1)
        elw(j)=elw(j-1)-cor(2,j)*uab(2,j)/grav*dy(2,j-1)
      end do
C
C     Adjust boundary elevations so that they are zero in the middle
C     of the channel:
C
      elejmid=ele(jmm1/2)
      elwjmid=elw(jmm1/2)
      do j=2,jmm1
        ele(j)=(ele(j)-elejmid)*fsm(im,j)
        elw(j)=(elw(j)-elwjmid)*fsm(2,j)
      end do
!
!tne:!wad:!lyo:wad:
      do i=1,im
         eln(i)=(0.-hhi)*fsm(i,jm)
         els(i)=(0.-hhi)*fsm(i,2)
      enddo
      do j=1,jm
         ele(j)=(ele(j)-hhi)*fsm(im,j)
         elw(j)=(elw(j)-hhi)*fsm(2,j)
      enddo
C
C     Set thermodynamic boundary conditions (for the seamount
C     problem, and other possible applications, lateral thermodynamic
C     boundary conditions are set equal to the initial conditions and
C     are held constant thereafter - users may, of course, create
C     variable boundary conditions):
C
      do k=1,kbm1
C
        do j=1,jm
          tbe(j,k)=tb(im,j,k)
          tbw(j,k)=tb(1,j,k)
          sbe(j,k)=sb(im,j,k)
          sbw(j,k)=sb(1,j,k)
        end do
C
        do i=1,im
          tbn(i,k)=tb(i,jm,k)
          tbs(i,k)=tb(i,1,k)
          sbn(i,k)=sb(i,jm,k)
          sbs(i,k)=sb(i,1,k)
        end do
C
      end do
C
      return
C
      end