program octa_topo

      implicit none
      include 'param_o.h'
      include 'mpif.h'

      real,dimension(0:im+1,0:jm+1,nm)::hlam,hphi
      real,dimension(0:im+1,0:jm+1,nm)::sm,hgt,agt
      integer,parameter::ntiles=18*36
      real,dimension(ntiles)::tfn,tfx,tln,tlx
      character*8,dimension(ntiles)::ctile
      real,dimension(:,:,:,:),allocatable::hsb,hsb2,ndb,hsdum,nddum
      real,dimension(:,:,:),allocatable::hgtsub,hgts1
      real,dimension(:,:,:,:,:),allocatable::hssb,ndsb,hssdum,ndsdum
      real::awlo,aelo,asla,anla,xlon,ylat,shs,alah,sumtx,sumtxp,sumty  &
           ,sumtyp
      integer::n,nlon,nlat,irecl,nrrd,i,j,i2,j2,n2,k2,kqsb,nphd,nsbhd  &
              ,r,k
      real,parameter::dlo=1./120.,rad=3.1415926/180.
      character*255::fname
      integer*2,dimension(1200,1200)::iht
      real,dimension(4)::hs,ndpsb
       real,parameter::flon0=270.  &
      , flat0=-90.,scen=1
      CHARACTER (len=39) :: tfile=& 
          '../data_in/grid/SOCT0004.DAT'
      integer::ierr,mype,npes,nstart,nend,ndata,ndata2,hdl,irecv,isend, &
             istat(MPI_STATUS_SIZE),ks1,ks2

!-----------------------------------------------------------------------
!***
!***  Initial MPI calls
!***
      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,mype,ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,npes,ierr)


      allocate(hssb(0:im+1,0:jm+1,nm,4,4),ndsb(0:im+1,0:jm+1,nm,4,4))
      allocate(hsb(0:im+1,0:jm+1,nm,4),ndb(0:im+1,0:jm+1,nm,4))
      allocate(hsb2(0:im+1,0:jm+1,nm,4))
      allocate(hgtsub(0:im+1,0:jm+1,nm))
      allocate(hgts1(0:im+1,0:jm+1,nm))
      allocate(hsdum(0:im+1,0:jm+1,nm,4),nddum(0:im+1,0:jm+1,nm,4))
      CALL insoct2(flon0, flat0, scen, tfile)

      open(1,file="../data_in/grid/hpos.dat",  &
          form="unformatted")
      read(1)hlam,hphi
      close(1)

      hlam=hlam*180./3.1415926
      hphi=hphi*180./3.1415926
      do n=1,nm
      do i=1,im
      do j=1,jm
        if(hlam(i,j,n).gt.180)then
	  hlam(i,j,n)=hlam(i,j,n)-360
        endif
      enddo
      enddo
      enddo


!
!
!    get the land/sea mask
!
!    5400,2700 for 4 minute global sea mask file
!    
!
        call get_sea(sm,5400,2700,hphi,hlam,im,jm,nm)

	call bocoh(sm,im,jm)

!
!
!    get the file name, max and min longlitude and latitude
!    of the topo tile files
!

        call getTopoTileInfo(ntiles,tfn,tfx,tln,tlx,ctile)

	hsb=0.
	ndb=0
	print *,ntiles

        k=ntiles/npes
	if(mod(ntiles,npes).ne.0)then
	  k=k+1
        endif

	nstart=1+mype*k
	nend=(mype+1)*k
	if(nend.gt.ntiles)then
	  nend=ntiles
        endif

	do n=nstart,nend
	  if(mype.eq.0)then
	    print *,n,nstart,nend
          endif

	  awlo=tln(n)
	  aelo=tlx(n)
	  asla=tfn(n)
	  anla=tfx(n)

	  nlon=int((aelo-awlo)/dlo)+1
	  nlat=int((anla-asla)/dlo)+1

	  irecl=nlon*2

          fname='../data_in/static/'//ctile(n)
	  open(1,file=fname,status='old',form='unformatted', &
	    access='direct',recl=irecl)
          do nrrd=1,nlat
	    read(1,rec=nrrd) (iht(i,nrrd),i=1,nlon)
	    do i=1,nlon
	      call swapval(iht(i,nrrd))
            enddo
          enddo
	  close(1)

	  do j=1,nlat
	  do i=1,nlon
	    xlon=(awlo+(i-1)/120.)*rad
	    ylat=(anla-(j-1)/120.)*rad
	    call lltooc(xlon,ylat,i2,j2,n2,k2,ks1,ks2)
	    hsb(i2,j2,n2,k2)=hsb(i2,j2,n2,k2)+iht(i,j)
	    hsb2(i2,j2,n2,k2)=hsb2(i2,j2,n2,k2)+iht(i,j)*iht(i,j)
	    ndb(i2,j2,n2,k2)=ndb(i2,j2,n2,k2)+1
	    hssb(i2,j2,n2,ks1,ks2)=hssb(i2,j2,n2,ks1,ks2)+iht(i,j)
	    ndsb(i2,j2,n2,ks1,ks2)=ndsb(i2,j2,n2,ks1,ks2)+1
          enddo
	  enddo
	  enddo

          allocate(hssdum(0:im+1,0:jm+1,nm,4,4),    &
	         ndsdum(0:im+1,0:jm+1,nm,4,4))
	  hssdum=0.
	  ndsdum=0.


          ndata=(im+2)*(jm+2)*nm*4
          ndata2=(im+2)*(jm+2)*nm*4*4
	  if(mype.eq.0)then
	    do n=1,npes-1
	       print *,'n=',n
	        call MPI_IRECV(hsdum, ndata, MPI_REAL, n, n, &
                      MPI_COMM_WORLD, hdl, irecv)
                call MPI_WAIT( hdl, istat, ierr )
 	        call MPI_IRECV(nddum, ndata, MPI_REAL, n, n, &
                      MPI_COMM_WORLD, hdl, irecv)
                call MPI_WAIT( hdl, istat, ierr )
		hsb=hsb+hsdum
		ndb=ndb+nddum
	        call MPI_IRECV(hsdum, ndata, MPI_REAL, n, n, &
                      MPI_COMM_WORLD, hdl, irecv)
                call MPI_WAIT( hdl, istat, ierr )
		hsb2=hsb2+hsdum

	        call MPI_IRECV(hssdum, ndata2, MPI_REAL, n, n, &
                      MPI_COMM_WORLD, hdl, irecv)
                call MPI_WAIT( hdl, istat, ierr )
 	        call MPI_IRECV(ndsdum, ndata2, MPI_REAL, n, n, &
                      MPI_COMM_WORLD, hdl, irecv)
                call MPI_WAIT( hdl, istat, ierr )
		hssb=hssb+hssdum
		ndsb=ndsb+ndsdum
	       print *,'n=',n
             enddo
          else
	     call MPI_ISEND(hsb, ndata,MPI_REAL, 0, mype,  &
                     MPI_COMM_WORLD, hdl, isend)
             call MPI_WAIT( hdl, istat, ierr )
 	     call MPI_ISEND(ndb, ndata,MPI_REAL, 0, mype,  &
                     MPI_COMM_WORLD, hdl, isend)
             call MPI_WAIT( hdl, istat, ierr )
	     call MPI_ISEND(hsb2, ndata,MPI_REAL, 0, mype,  &
                     MPI_COMM_WORLD, hdl, isend)
             call MPI_WAIT( hdl, istat, ierr )

	     call MPI_ISEND(hssb, ndata2,MPI_REAL, 0, mype,  &
                     MPI_COMM_WORLD, hdl, isend)
             call MPI_WAIT( hdl, istat, ierr )
 	     call MPI_ISEND(ndsb, ndata2,MPI_REAL, 0, mype,  &
                     MPI_COMM_WORLD, hdl, isend)
             call MPI_WAIT( hdl, istat, ierr )

          endif
	  deallocate(hssdum,ndsdum)
   

          if(mype.eq.0)then
	  print *,'begin processing sfc ...'
          do n=1,nm
	  do j=1,jm
	  do i=1,im
	  do kqsb=1,4
	    if(ndb(i,j,n,kqsb).ne.0)then
	      hsb(i,j,n,kqsb)=hsb(i,j,n,kqsb)/ndb(i,j,n,kqsb)
	      hsb2(i,j,n,kqsb)=hsb2(i,j,n,kqsb)/ndb(i,j,n,kqsb)
	      hsb2(i,j,n,kqsb)=hsb2(i,j,n,kqsb)-hsb(i,j,n,kqsb)**2
	      hsb2(i,j,n,kqsb)=sqrt(hsb2(i,j,n,kqsb))
            endif
	    do ks2=1,4
	    if(ndsb(i,j,n,kqsb,ks2).ne.0)then
	      hssb(i,j,n,kqsb,ks2)=hssb(i,j,n,kqsb,ks2)/ndsb(i,j,n,kqsb,ks2)
            endif
            enddo
          enddo
          enddo
          enddo
          enddo
	  print *,'finish processing sfc ...'
        call bocohsub(hsb,im,jm)
        call bocohsub(hsb2,im,jm)
        call bocohsub(ndb,im,jm)
        call bocohsub2(hssb,im,jm)
        call bocohsub2(ndsb,im,jm)

        print *,'zeff ...'
        call zeff(hssb,im,jm,nm)
        print *,'finished zeff ...'

!
!   Initialize the height 
!

       hgt=0.

!
!  Calculate silhouette height
!  
!  It is the average of the max(hsb(4),hsb(1)), max(hsb(3),hsb(2))
!  max(hsb(4),hsb(3)) and max(hsb(1),hsb(2)). hsb(i) is the average
!  height of the subbox i.
! 
!
        do n=1,nm
        do j=0,jm+1
        do i=0,im+1
            hs(1)   =amax1(hsb(i,j,n,4),hsb(i,j,n,1))
            NDPSB(1)=      NDB(i,j,n,4)+NDB(i,j,n,1)
            HS   (2)=AMAX1(HSB(i,j,n,3),HSB(i,j,n,2))
            NDPSB(2)=      NDB(i,j,n,3)+NDB(i,j,n,2)
            HS   (3)=AMAX1(HSB(i,j,n,4),HSB(i,j,n,2))
            NDPSB(3)=      NDB(i,j,n,4)+NDB(i,j,n,2)
            HS   (4)=AMAX1(HSB(i,j,n,1),HSB(i,j,n,3))
            NDPSB(4)=      NDB(i,j,n,1)+NDB(i,j,n,3)   
          
	    NPHD=4
            DO r=1,4
              IF(NDPSB(r).EQ.0) NPHD=NPHD-1
            enddo
     
        if(nphd.ne.0) then
           SHS=0.
           DO r=1,4
            SHS=SHS+HS(r)
           enddo
           HGT(i,j,n)=SHS/NPHD
        endif
      enddo
      enddo
      enddo

      do n=1,nm
      do j=0,jm+1
      do i=0,im+1
        hgts1(i,j,n)=0.25*(hsb2(i,j,n,1)+hsb2(i,j,n,2)+hsb2(i,j,n,3)  &
	              +hsb2(i,j,n,4))
      enddo
      enddo
      enddo

      call bocoh(hgts1,im,jm)

      do n=1,nm
      do j=1,jm-1
      do i=1,im-1
        hgtsub(i,j,n)=0.25*(hgts1(i,j,n)+hgts1(i,j+1,n)+hgts1(i+1,j,n)  &
	              +hgts1(i+1,j+1,n))
      enddo
      enddo
      enddo


      call bocoh(hgt,im,jm)
      call bocoh(hgtsub,im,jm)
    

!-----------------------------------------------------------
!***  CALCULATE AND SAVE AVERAGE ELEVATIONS
!

      AGT=0.

      do r=1,nm
      DO j=1,jm
      do i=1,im
         NSBHD=4
         DO n=1,4
           IF(NDB(i,j,r,n).EQ.0)NSBHD=NSBHD-1
         enddo

         IF(NSBHD.ne.0)then
            AGT(i,j,r)=(HSB(i,j,r,1)+HSB(i,j,r,2)+HSB(i,j,r,3)   &
	       +HSB(i,j,r,4))/NSBHD
         endif
     enddo
     enddo
     enddo
     call bocoh(agt,im,jm)
     print *,'silhouette/mean mountains...'

!***
!***  THE silhouette/mean MOUNTAINS:  AT LAND POINTS WITH CONCAVE
!***  ACTUAL (AVERAGE) TOPOGRAPHY, REPLACE THE
!***  SILHOUETTE ELEVATION WITH THE AVERAGE ELEVATION
!***  
      do r=1,nm
      DO j=1,jm
      do i=1,im
        IF(SM(i,j,r).le.0.5)then
           ALAH=AGT(i+1,j+1,r)+AGT(i-1,j+1,r)+AGT(i+1,j-1,r)  &
            +AGT(i-1,j-1,r)+2.*(AGT(i+1,j,r)                 &
            +AGT(i-1,j,r)+AGT(i,j+1,r)+AGT(i,j-1,r))         &
            -12.*AGT(i,j,r)
          IF(ALAH.GE.0.)HGT(i,j,r)=AGT(i,j,r)
        endif
      enddo
      enddo
      enddo
!
!     If K is a point seen as a "valley" when looking at three-point
!     averages in any one of the four directions, and it is at the
!     same time not higher than each of its four nearest neighbors,
!     and also not higher than each of its four second-nearest
!     neighbors, choose always the mean elevation (23 June 97)
!
!     (Saddle points are expected to be frequently declared mean as
!     a result)
!
      do r=1,nm
      DO j=1,jm
      do i=1,im
       IF(SM(i,j,r).le.0.5)then
       if (.not.(agt(i,j,r).gt.agt(i+1,j,r) .and. agt(i,j,r).gt.agt(i-1,j,r) .and. &
         agt(i,j,r).gt.agt(i,j+1,r) .and. agt(i,j,r).gt.agt(i,j-1,r))) then 
       if (.not.(agt(i,j,r).gt.agt(i+1,j+1,r) .and. agt(i,j,r).gt.agt(i-1,j+1,r) .and. &
         agt(i,j,r).gt.agt(i+1,j-1,r) .and. agt(i,j,r).gt.agt(i-1,j-1,r)))then   

           sumtx =agt(i-1,j+1,r)+agt(i,j,r)+agt(i+1,j-1,r)
           sumtxp=agt(i-1,j,r)+agt(i,j,r)+agt(i+1,j,r)
           sumty =agt(i+1,j+1,r)+agt(i,j,r)+agt(i-1,j-1,r)
           sumtyp=agt(i,j+1,r)+agt(i,j,r)+agt(i,j-1,r)

          if (sumtx .lt.agt(i+1,j,r)+agt(i+1,j+1,r)+agt(i,j+1,r) .and.  &
             sumtx .lt.agt(i-1,j,r)+agt(i-1,j-1,r)+agt(i,j-1,r))        &
               hgt(i,j,r)=agt(i,j,r)
          if (sumtxp.lt.agt(i-1,j+1,r)+agt(i,j+1,r)+agt(i+1,j+1,r) .and. &
              sumtxp.lt.agt(i+1,j-1,r)+agt(i,j-1,r)+agt(i-1,j-1,r))      &
               hgt(i,j,r)=agt(i,j,r)
          if (sumty .lt.agt(i,j+1,r)+agt(i-1,j+1,r)+agt(i-1,j,r) .and.   &
              sumty .lt.agt(i+1,j-1,r)+agt(i,j-1,r)+agt(i+1,j,r))        &
               hgt(i,j,r)=agt(i,j,r)
          if (sumtyp.lt.agt(i+1,j+1,r)+agt(i+1,j,r)+agt(i+1,j-1,r) .and. &
              sumtyp.lt.agt(i-1,j,r)+agt(i-1,j-1,r)+agt(i-1,j+1,r))      &
               hgt(i,j,r)=agt(i,j,r)
       endif
       endif
       endif
      enddo
      enddo
      enddo

       call bocoh(hgt,im,jm)

      if(mype.eq.0)then
      do i=1,im
      do j=1,jm
      do n=1,nm
!        print *,i,j,n,hgt(i,j,n)
      enddo
      enddo
      enddo
      endif


!    the silhouette/mean topography is done

      write(6,'('' silhouette/mean topography is done'')')

      open(1,file='../data_in/grid/topo.dat',  &
           status='unknown',form='unformatted')
      write(1) hgt,sm
      close(1)

      endif

      open(1,file='../data_in/grid/topo_sub.dat',form='unformatted')
      write(1)hgtsub
      close(1)

      deallocate(hssb,ndsb)
      deallocate(hsb,ndb)
      deallocate(hsb2,hgtsub,hgts1)
      deallocate(hsdum,nddum)
 
      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
      call MPI_FINALIZE(ierr)              
      
      end program octa_topo