subroutine zeff(hssb,im,jm,nm)
    
      implicit none

      integer::im,jm,nm
      real,dimension(0:im+1,0:jm+1,nm,4,4)::hssb
      real,dimension(0:im+1,0:jm+1,nm,4)::ds
      real,dimension(0:im+1,0:jm+1,nm)::rarea
      real::ao(4)
      real::exvl(6)
      real::exvlc(4)
      real::ha(4)
      real::hghmn
      real::hhadir(4)
      real::hhdir
      real::hi
      real::hmax2
      real::hmin1
      real::hmin3
      real::hob
      real::hodir
      real::htadir(4)
      real::htdir
      integer::i
      integer::i1
      integer::iex
      integer::ilob
      integer::imax
      integer::imin
      integer::imn
      integer::is
      integer::isrch
      integer::issb
      integer::j
      integer::j1
      integer::jmax
      integer::jmin
      integer::js
      integer::jssb
      integer::llsb
      integer::m
      integer::n
      integer::nadir(4)
      integer::ndir
      integer::ndirp
      integer::nexvlf
      integer::nhghmn
      integer::nhmm
      integer::nlim
      integer::nlob
      logical::isCorner
      real::saob
      real::shhob
      real::sl
      real::wtarea
      real::z0eff(0:im+1,0:jm+1,nm,4)
      real::zas(0:im+1,0:jm+1,nm,4)
      real::zeffl
      real::zefflc
      real::zha(0:im+1,0:jm+1,nm,4)
      real,parameter::dmin=1.0e-6,z0=0.15,rz0=1./z0,cd=0.6,vk=0.4,  &
                     rvk2=1./(vk**2),wtamin=dmin
      real,parameter::hmin=2.0*z0+dmin

      open(1,file='fzeff.dat', form="unformatted")
      read(1)rarea,ds
      close(1)

         
	  do n=1,nm
	  do i=1,im
	  do j=1,jm
	  do m=1,4
            hhdir=0.0
            htdir=0.0
            ndir=0.0
            imin=1
            if (mod(m,2).eq.1) then
!             cardinal direction, diagonal on grid
              imax=7
            else
!             diagonal direction, along grid
              imax=4
            endif
!
!           loop to select a line of subboxes
!
            do is=1,imax
              if (mod(m,2).eq.1) then
!               cardinal direction, diagonal on grid
                jmin=max(0,is-4)
                jmax=min(5,is+1)
              else
!               diagonal direction, along grid
                jmin=0
                jmax=5
              endif
              llsb=jmax-jmin+1
              nexvlf=0
!
!             loop to search for extrema along the line of subboxes
!
              do 4223 js=jmin,jmax
!               select subbox indices by direction indicator
                if (m.eq.3) then
                  issb=1+is-js
                  jssb=js
                endif
                if (m.eq.4) then
                  issb=js
                  jssb=is
                endif
                if (m.eq.1) then
                  issb=js
                  jssb=4-is+js
                endif
                if (m.eq.2) then
                  issb=is
                  jssb=js
                endif

                i1=i
		j1=j
!               reindex to the interior of a box
                if (issb.lt.1) then
                  issb=issb+4
		  i1=i1-1
                endif
                if (jssb.lt.1) then
                  jssb=jssb+4
		  j1=j1-1
                endif
                if (issb.gt.4) then
                  issb=issb-4
		  i1=i1+1
                endif
                if (jssb.gt.4) then
                  jssb=jssb-4
		  j1=j1+1
                endif
!               find minima and maxima
                if (js.eq.jmin) then
!                 first point, taken to be first minimum
                  hmin1=hssb(i1,j1,n,issb,jssb)
                  isrch=1
                else
!                 any other point
                  if (isrch.eq.1) then
!                   looking for the first minimum
                    if (hmin1.gt.hssb(i1,j1,n,issb,jssb)) then
!                     lower the first minimum
                      hmin1=hssb(i1,j1,n,issb,jssb)
                    elseif (hmin1.lt.hssb(i1,j1,n,issb,jssb)) then
!                     we have a minimum, now look for maximum
                      nexvlf=nexvlf+1
                      exvl(nexvlf)=hmin1
                      hmax2=hssb(i1,j1,n,issb,jssb)
!		      print *,hmax2,i1,j1,n,issb,jssb,hssb(i1,j1,n,issb,jssb)
                      isrch=2
                    endif
                  elseif (isrch.eq.2) then
!                   looking for a maximum
                    if (hmax2.lt.hssb(i1,j1,n,issb,jssb)) then
!                     raise the maximum
                      hmax2=hssb(i1,j1,n,issb,jssb)
                    elseif (hmax2.gt.hssb(i1,j1,n,issb,jssb)) then
!                     we have a maximum, now look for minimum
                      hmin3=hssb(i1,j1,n,issb,jssb)
                      isrch=3
                    endif
                  elseif (isrch.eq.3) then
!                   looking for a minimum
                    if (hmin3.gt.hssb(i1,j1,n,issb,jssb)) then
!                     lower the minimum
                      hmin3=hssb(i1,j1,n,issb,jssb)
                    endif
                    if ((hmin3.lt.hssb(i1,j1,n,issb,jssb)).or.  &
                        (js.eq.jmax)) then
                      nexvlf=nexvlf+1
                      exvl(nexvlf)=hmax2
!		      print *,i,j,n,m,nexvlf,hmax2,exvl(nexvlf)
                      nexvlf=nexvlf+1
                      exvl(nexvlf)=hmin3
!		      print *,i,j,n,m,nexvlf,exvl(nexvlf)
!                     move the minimum value to front of next mt
                      hmin1=hmin3
!                     now we are looking for another maximum
                      hmax2=hssb(i1,j1,n,issb,jssb)
                      isrch=2
                    endif
                  endif
                endif
 4223         continue
 
              if (nexvlf.ge.5) then
!
!  calculate z0eff of the line just processed,
!  assuming that the width of the obstacles is unity
!
!  sl = length of the line along which obstacles are calculated
!  nlob = n of the peak of the last obstacle along the line
!
               sl=llsb*ds(i,j,n,m)
!
                nlob=nexvlf-1
                saob=0.0
                shhob=0.0
                do 4235 ilob=2,nlob,2
                  hi=exvl (ilob)-0.5*(exvl (ilob-1)+exvl (ilob+1))
                  saob=saob+hi
                  shhob=shhob+hi**2
 4235           continue
                if (saob.gt.hmin) then
                  hob=shhob/saob
                else
                  hob=hmin
                  saob=hmin
                  shhob=hmin**2
                endif
!
!  calculate a tentative value of z0eff for this line
!
	if (0.5*rz0*hob.lt.0) write(6,*) 'neg alog', 0.5*rz0*hob

                if( hob==0.3 .and. 0.5*rz0*hob==1.0 ) then
                   hob = 0.30001 
                end if

                wtarea=0.5*cd*rvk2*saob/sl +1.0/((alog(0.5*rz0*hob))**2)
		if (wtarea.lt.wtamin) then
                  print *,' nonpositive tentative wtarea',wtarea,i,j,n,m
                  wtarea=wtamin
                endif
                zeffl =0.5*hob/exp(1.0/sqrt(wtarea))

!
!  of the one or more inside minima, find the highest
!
 4243           continue
                hghmn=exvl(3)
                nhghmn=3
                nlim=3
!
!  nlim = n of the last minimum
!
!  nlim = 5 if there is only one inside minimum; it is the highest
!
                if(nexvlf.gt.5) then
                  nlim=nexvlf-2
                  do 4245 imn=5,nlim,2
                    if (hghmn.lt.exvl(imn)) then
                      hghmn=exvl(imn)
                      nhghmn=imn
                    endif
 4245             continue
                endif
!
!  the highest inside minimum has been found.
!  does removing it and compressing exvl increase z0eff of the line?
!
                nhmm=nhghmn-1
                do 4251 iex=1,nhmm                  
		exvlc(iex)=exvl(iex)
 4251           continue
                if (exvl(nhghmn+1).gt.exvl(nhmm)) then
                  exvlc(nhmm)=exvl(nhghmn+1)
                endif
!
!  the higher of the two maxima has been selected and stored
!  now transfer the remainder of exvl to exvlc
!
                do 4253 iex=nhghmn,nlim
                  exvlc(iex)=exvl(iex+2)
 4253           continue
!
!  calculate z0eff of the compressed array, zefflc
!
                nlob=nlim-1
                saob=0.0
                shhob=0.0
                do 4255 ilob=2,nlob,2
                  hi=exvlc(ilob)-0.5*(exvlc(ilob-1)+exvlc(ilob+1))
                  saob=saob+hi
                  shhob=shhob+hi**2
 4255           continue
                if (saob.gt.hmin) then
                  hob=shhob/saob
                else
                  hob=hmin
                  saob=hmin
                  shhob=hmin**2
                endif
!
!  calculate a tentative value of z0eff for this line
!
                wtarea=0.5*cd*rvk2*saob/sl +1.0/((alog(0.5*rz0*hob))**2)
                  
                if (wtarea.lt.wtamin) then
                  print *,'nonpositive tentative wtarea',wtarea,i,j,n,m
                  wtarea=wtamin
                endif
                zefflc=0.5*hob/exp(1.0/sqrt(wtarea))
!
                if (zeffl.lt.zefflc) then
!
!  switch to compressed exvl and either exit the consolidation code
!  or search for another minimum to remove
!
                  do 4257 iex=nhmm,nlim
                    exvl(iex)=exvlc(iex)
 4257             continue
                  nexvlf=nlim
!
                  if (nexvlf.lt.5) go to 4263
!
                  zeffl =zefflc
                  go to 4243
                else
                  go to 4263
                endif
!
              endif
!
 4263         continue
!
!  There has been at most only one obstacle to start with,
!  or consolidation of obstacles is completed.
!  Add obstacles, if any, to the totals
!
      if (nexvlf.ge.3) then
          nlob=nexvlf-1
          do 4267 ilob=2,nlob,2
              hodir=exvl(ilob)-0.5*(exvl(ilob-1)+exvl(ilob+1))
!	      print *,ilob,hodir,exvl(ilob),exvl(ilob-1),exvl(ilob+1)
             if (hodir .eq. 0) then
              write(6,*) 'increasing hodir!!! '
              hodir=hmin
             endif
             htdir=htdir+hodir
             hhdir=hhdir+hodir**2
             ndir=ndir+1
 4267       continue
        endif

	enddo       !finished calculation of obstacles

!           in case we didn't find any obstacles
            if (ndir.eq.0) then
              ndir=1
              htdir=hmin
              hhdir=hmin**2
            endif
!           fill directional obstacle totals
            nadir(m)=ndir
            htadir(m)=htdir
            hhadir(m)=hhdir
            ndirp=ndir+1

          enddo


          do 4437 m=1,4
!
!         find the total area of obstacles in each direction
!
            ao(m)=htadir(m)*ds(i,j,n,m)*rarea(i,j,n)
!	    print *,ao(m),htadir(m),ds(i,j,n,m),rarea(i,j,n)
!
!         find the average heights of obstacles, use hmin to avoid zero
!           the average is weighted by the height of the obstacle
!
	if (htadir(m) .lt. 0) write(6,*) 'divide by zero!'
            ha(m)=hhadir(m)/htadir(m)
!            print *,ha(m),hhadir(m),htadir(m)
 4437     continue
!
!         calculate z0eff
!
          do 4552 m=1,4
            zha(i,j,n,m)=ha(m)
            zas(i,j,n,m)=ao(m)
!
! test for out-of-range values
!
            if (zha(i,j,n,m).le.hmin) then
              zha(i,j,n,m)=hmin
            endif

            if (zas(i,j,n,m).lt.0.0) then
              write(6,'('' negative zha(i,j,n,m,1)'',f20.10,2i8)')  &
                zas(i,j,n,m),i,j,n,m
              zas(i,j,n,m)=0.0
            endif
!
! calculate the effective z0
! from zas and zha
!
            wtarea=0.5*cd*rvk2*zas(i,j,n,m)  &
	         +1.0/((alog(0.5*rz0*zha(i,j,n,m)))**2)
              
            if (wtarea.lt.wtamin) then
              write(6,'('' nonpositive wtarea'',f20.10,2i8)')  &
                wtarea,i,j,n,m
              wtarea=wtamin
            endif
            z0eff(i,j,n,m)=0.5*zha(i,j,n,m)/exp(1.0/sqrt(wtarea))

 4552     continue
          enddo
          enddo
          enddo

	  do n=1,nm
	  do i=1,im
	  do j=1,jm
	    if(isCorner(i,j,n,im,jm,nm).ne.0)then
	      z0eff(i,j,n,:)=0.15
            endif
          enddo
          enddo
          enddo

	  open(1,file='zeff.dat', form="unformatted")
          write(1)z0eff
          close(1)


          end subroutine zeff