ccccccccccccccccccccccccccccccccccccccccc
      subroutine topolog(x,y)
ccccccccccccccccccccccccccccccccccccccccc
      include 'param.nml'
      include 'msg.inc'

......

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c-->  itopo=1   - read topo array from tape and redistribute it to pe's
c-->  ntop,mtop - dimension of topo array from tape
c-->  if (ntop.ne.n) or (mtop.ne.m) then we need interpolation
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      parameter(itopo=0,ntop=n,mtop=m)

      dimension zstmp(np,mp),    ! tmp array for transfer between pe's
     .         zsfull(n,m),      ! topo array on full domain (size n x m)
     .         zstape(ntop,mtop) ! topo araay on tape (size ntop x mtop)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                  MOUNTAIN SHAPE FUNCTIONS
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
construct fortran statement functions for the boundary shape
create time dependence of the lower boundary; if time dependent lower boundary
condition is used for initialisation make sure tend is multiple of dt
construct fortran statement functions for the boundary shape
c++++++++++++++++++++++++++
C Williamson et al. test case
      rd(xln,ylt)=sqrt(((xln-xml0)/xml)**2+((ylt-yml0)/yml)**2)
      rb(xln,ylt)=sqrt(((xln+xml0)/xml)**2+((ylt+yml0)/yml)**2)
      rc(xln,ylt)=sqrt(((xln     )/xml)**2+((ylt-xml0)/yml)**2)
c     profm(rad)=amax1(0.,1.-rad)
      profm(rad)=.5*(1.+cos(pi*rad))
      profb(rad)=rad
c++++++++++++++++++++++++++
c special for the conical mountain on the pole ccccc
c     rd(xln,ylt)=abs(ylt-y0)
c     profm(rad)=amax1(0., 1.-gamm*rad/rnot)
c     rnot = 2*acos(-1.)/128. * 10.
c     gamm=1.


c++++++++++++++++++++++++++
c     rd(xx,yy)=sqrt((xx/xml)**2+j3*(yy/yml)**2)
c     fi(rad)=amp*exp(-rad**2)
c     fi(rad)=.5*(1.+cos(pi*rad))
c     fi(rad)=amp/(1.+rad**2)
! klemp mountain
c      fi(rad)=amp*exp(-(rad/xml)**2)*(cos(pi*rad*.25e-3))**2
c membranas
c     fi(rad)=-.5*(1.+cos(pi*rad))
c     fid(rad)=.5*pi*sin(pi*rad)
! hyd. jump
c     fi(rad)=amp/(sqrt(1.+rad**2))**3
! cylindrical annulus
c     fi(rad,is,r1,r2)=(rad**2/r1**2)*
c    .  (1.-( ( ((r2/r1)**(is+2))-1.)*((rad/r1)**(is-2))
c    .       +( ((r2/r1)**(is-2))-1.)*((r2/rad)**(is+2)) )
c    .      /(((r2/r1)**(2*is))-1.) )
      fid(rad)=0.
 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
construct statement functions for the horizontal coordinate mapping

      parameter(syfact=2., iyorder=5,
     +          sycoef1=1./syfact,sycoef2=(syfact-1.)/syfact)
      parameter(sxfact=2., ixorder=5,
     +          sxcoef1=1./sxfact,sxcoef2=(sxfact-1.)/sxfact)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      domains are:    -1 < yln  < 1 and -1 < xln  < 1
c      and ranges are: -1 < yphy < 1 and -1 < xphy < 1
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c IDENTITY TRANSFORMATION (DEFAULT)
      xphy(xln,yln,tln)=xln                                   ! X0
      yphy(xln,yln,tln)=yln                                   ! Y0
c++++++++++++++++++++++++++
c --- Y GRID STRETCHING
c     syfact - stretching factor at y00
c     iyorder - order of polynomial function (odd order >= 3 required)
c     sycoef1 - inverse stretching factor at equator
c     sycoef2 - nonlinear coefficient
c     yphy(xln,yln,tln)=sycoef1*yln+sycoef2*yln**iyorder        ! Y_[-1,1]
C     yphy(xln,yln,tln)=yln*(sycoef1 + (1.-sycoef1)*yln*yln*    ! Y_cyc
C    +                           (10.-15.*yln+6.*yln*yln) )     ! Y_[0,1]
c++++++++++++++++++++++++++
c --- X GRID STRETCHING
c     sxfact - stretching factor at x00
c     ixorder - order of polynomial function (odd order >= 3 required)
c     sxcoef1 - inverse stretching factor at equator
c     sxcoef2 - nonlinear coefficient
c     xphy(xln,yln,tln)=sxcoef1*xln+sxcoef2*xln**ixorder        ! X_[-1,1]
cc    xphy(xln,yln,tln)=xln*(sxcoef1 + (1.-sxcoef1)*xln*xln*    ! X_cyc
cc   +                           (10.-15.*xln+6.*xln*xln) )     ! X_[0,1]
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


c--- compute horizontal grid point "physical" locations
      do j=1,mp
        do i=1,np
         ia=(npos-1)*np + i
         ja=(mpos-1)*mp + j
c         xcr(i,j)=xnorm*xphy(x(ia)*xnormi,y(ja)*ynormi,tt)
c         ycr(i,j)=ynorm*yphy(x(ia)*xnormi,y(ja)*ynormi,tt)
          xcr(i,j)=x(ia)*(rds*isphere+(1-isphere))
          ycr(i,j)=y(ja)*(rds*(isphere+icylind)+(1-(isphere+icylind)))
        enddo
      enddo


compute topografy
      angep=ang-1.e-10
      cang=cos(pi*angle/180.)
      sang=sin(pi*angle/180.)
      yj0=pi/180.*ang
      x0 = 1.5*pi*isphere+(1-isphere)*0.5*(x(1)+x(n))
      y0 = pi/6. *isphere+(1-isphere)*0.5*(y(1)+y(m))*j3
      dd=0.0133
      do 20 j=1,mp
      do 20 i=1,np
         ia=(npos-1)*np + i
         ja=(mpos-1)*mp + j

c --- compute trig and coriolis functions
         yj=ycr(i,j)*rdsi
         icoord=isphere+icylind
         cosa(i,j)=cos(yj)*icoord+(1-icoord)*1.
         sina(i,j)=sin(yj)*icoord+(1-icoord)*0.
         tnga(i,j)=tan(yj)*icoord+(1-icoord)*0.
         fcr2(i,j)=icoord*fcr0*cosa(i,j)
     *         +(1-icoord)*fcr0*(cos(yj0)-btpl*sin(yj0)*yj)
         fcr3(i,j)=icoord*fcr0*sina(i,j)
     *         +(1-icoord)*fcr0*(sin(yj0)+btpl*cos(yj0)*yj)

c --- compute lower boundary deflection

         xx= cang*(x(ia)-x0)+sang*(y(ja)-y0)
         yy=-sang*(x(ia)-x0)+cang*(y(ja)-y0)
         rr=rd(xx,yy)
      zs(i,j)= 0.
c     zs(i,j)=amp*profm(rr)
      zs(i,j)=amp*cvmgm(0.,profm(rr),1.-rr)
      zsd(i,j)= 0.
      zh(i,j) = zb
      zhd(i,j)= 0.

c read topography from external file

      if (itopo.eq.1) then  !read topografy from tape
      if (mype.eq.0) then
         call readtopo0(zs,zstmp,zstape,zsfull,ntop,mtop)
      else
         call readtopok(zs,zstmp)
      end if
      endif

filter topography if needed

      nitrf=-4
      do itrf=1,nitrf
      call filzs(zh)
      enddo

create imersed boundary objects

      if(imrsb.eq.1) then
c     ddh=0.5*dd
      do k=1,lib
       do j=1,mibp
        do i=1,nibp
         ia=(npos-1)*np + i
         ja=(mpos-1)*mp + j
c-------------------------------
c        zsib(i,j,k)=-1.
c         ri=sqrt((x(ia)-x0)**2+j3*(y(ja)-y0)**2)
c         zcrl=(k-1)*dz
c         if(ri.le.ddh.and.zb-zcrl.le.10.*dd) zsib(i,j,k)=1.
c-------------------------------
          gij=zb/(zh(i,j)-zs(i,j))           ! define terrain following transformation
          zloc=(k-1)*dz                      ! compute unstretched uniform height
          zcrl=zloc/gij+zs(i,j)              ! height in terrain following coordinates

          zsib(i,j,k)=-1.                    ! default, no IMB points
          xx= cang*(x(ia)-x0)+sang*(y(ja)-y0)
          yy=-sang*(x(ia)-x0)+cang*(y(ja)-y0)
          if(itopo.eq.0) rrb=rb(xx,yy)       ! object shifted in X 
          if(itopo.eq.1) rrb=rc(xx,yy)       ! object shifted in Y
          zsib2(i,j)=amp*cvmgm(0.,1.,1.-rrb)
          if(zcrl.lt.zsib2(i,j)) then        ! check if grid height below object
          zsib(i,j,k)=1.                     ! add IMB points below zsib2 surface
          endif