program cube_topo

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

      integer,PARAMETER::M=64, NP=M+1  
     
!!!!!!     integer,PARAMETER::M=128, NP=M+1  !DRAGAN
    
    
      real,dimension(3,-1:NP, -1:NP)::xtab
      COMMON / QPAN / XTAB   
      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(0:im+1,0:jm+1,nm,4)::hsb,ndb,hsdum,nddum
      real,dimension(0:im+1,0:jm+1,nm,4,4)::hssb,ndsb
      real,dimension(:,:,:,:,:),allocatable::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,ks2
      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
      integer::ierr,mype,npes,nstart,nend,ndata,ndata2,hdl,irecv,isend, &
             istat(MPI_STATUS_SIZE),ks1,irtn

!-----------------------------------------------------------------------
!***
!***  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)

      call Infin3
!
! Read Table
!
      open (unit = 10, file = 'round.dat', form = 'unformatted')
      Read (10) XTAB  
      Close (10)  

      open(1,file="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)
        print *, "finish getting sea mask"
        do n=1,nm
        do i=1,im
        do j=1,jm
!          print *,i,j,sm(i,j,n)
        enddo
        enddo
        enddo

        call bocoh_cb(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)

        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

        hssb=0.
        ndsb=0.
        hsb=0.
        ndb=0.
        hsdum=0.
        nddum=0.
        do n=nstart,nend
          awlo=tln(n)
          aelo=tlx(n)
          asla=tfn(n)
          anla=tfx(n)

          if(mype.eq.0)then
            print *,n
          endif  

          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 lltocb(xlon,ylat,i2,j2,n2,k2,ks1,ks2)
!TEST#################################################################
            if(ks2==0) then
              print *,'xlon,ylat,i2,j2,n2,k2,ks1,ks2=', &
                     xlon,ylat,i2,j2,n2,k2,ks1,ks2
            end if
!TEST#################################################################
            hsb(i2,j2,n2,k2)=hsb(i2,j2,n2,k2)+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

          ndata=(im+2)*(jm+2)*nm*4
          ndata2=(im+2)*(jm+2)*nm*4*4

          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.

          CALL MPI_REDUCE(hsb,hsdum,ndata,MPI_REAL,MPI_SUM,0  &
      ,               MPI_COMM_WORLD,IRtn)
          CALL MPI_REDUCE(ndb,nddum,ndata,MPI_REAL,MPI_SUM,0  &
      ,               MPI_COMM_WORLD,IRtn)
          CALL MPI_REDUCE(hssb,hssdum,ndata2,MPI_REAL,MPI_SUM,0  &
      ,               MPI_COMM_WORLD,IRtn)
          CALL MPI_REDUCE(ndsb,ndsdum,ndata2,MPI_REAL,MPI_SUM,0  &
      ,               MPI_COMM_WORLD,IRtn)
 
          if(mype.eq.0)then
             hsb=hsdum
             ndb=nddum
             hssb=hssdum
             ndsb=ndsdum
          endif

          deallocate(hssdum,ndsdum)

!        print *,'ndsb,',ndsb(21,25,1,1,4)
          if(mype.eq.0)then
          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)
            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
        call bocohsub_cube(hsb,im,jm)
        call bocohsub_cube(ndb,im,jm)
        call bocohsub2_cube(hssb,im,jm)
        call bocohsub2_cube(ndsb,im,jm)

  write(0,*)' before zeff'
        call zeff(hssb,im,jm,nm)
  write(0,*)' after 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
  write(0,*)' before bocoh_cb'
      call bocoh_cb(hgt,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
  write(0,*)' before bocoh_cb 2'
     call bocoh_cb(agt,im,jm)

!***
!***  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

  write(0,*)' before bocoh_cb 3'
       call bocoh_cb(hgt,im,jm)


!    the silhouette/mean topography is done

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

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

      endif

      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
      call MPI_FINALIZE(ierr)              
      
      end program cube_topo