program lev_extend

!  extend 3d field into land points

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
      implicit none

      integer :: imt, jmt, km, num_iters

      double precision, allocatable, dimension(:,:,:) :: FIELD
      integer, allocatable, dimension(:,:,:) :: MASK

      double precision :: land_val
 
      integer :: k, rec_length

      character*100 :: infile, outfile

!*******************************************************
!*******************************************************
      
      write(*,*)' enter nx, ny, nz'
      read(*,*) imt, jmt, km

      write(*,*)' enter name of real*8 input file'
      read(*,'(a100)') infile

      write(*,*)' enter name of real*8 output file'
      read(*,'(a100)') outfile

      write(*,*)' enter land value'
      read(*,*)land_val

      write(*,*)' enter number of iterations'
      read(*,*)num_iters

      allocate( FIELD(imt,jmt,km), MASK(imt,jmt,km))
      FIELD = 0.
      MASK = 1

      inquire(iolength=rec_length) FIELD(:,:,1)
      open(1,file=infile,access='direct',form='unformatted',
     &          recl=rec_length,status='unknown')
      do k = 1, km
        read(1,rec=k) FIELD(:,:,k)
        if(land_val > 0.) then
          where (FIELD(:,:,k) >= land_val) MASK(:,:,k) = 0
        elseif(land_val < 0.) then
          where (FIELD(:,:,k) <= land_val) MASK(:,:,k) = 0
        elseif(land_val == 0.) then
          where (FIELD(:,:,k) == land_val) MASK(:,:,k) = 0
        endif
      enddo
      close(1)

      call jacobi(FIELD,MASK,imt,jmt,km,num_iters)

      open(1,file=outfile,access='direct',form='unformatted',
     &          recl=rec_length,status='unknown')

      do k = 1, km
        write(1,rec=k) FIELD(:,:,k)
      enddo
      close(1)

      end

!----------------------------------------------------------------------

      subroutine jacobi(t,mask,imt,jmt,km,niter)

      implicit none

      integer :: imt, jmt, km, niter, init_val

      double precision, dimension(imt,jmt,km) :: t
      integer, dimension(imt,jmt,km) :: mask

      double precision, allocatable, dimension(:,:) :: tnew, told
      logical, allocatable, dimension(:,:) :: maskold

      integer :: kmin, kmax, ism, i, j, k, jj, iter
      double precision :: sum, average, eps, epsmax
!***********************************************************************

      allocate ( tnew(imt,jmt), told(imt,jmt), maskold(imt,jmt) )
      tnew = 0.
      told = 0.
      maskold = .false.

      kmin = 1
      kmax = km

      do k = kmin,kmax

      sum = 0.
      ism = 0
      do j = 1,jmt
         do i = 1,imt
            if(mask(i,j,k) .eq. 0) then
               maskold(i,j) = .false.
            else
               maskold(i,j) = .true.
               tnew(i,j) = t(i,j,k)
               sum = sum + tnew(i,j)
               ism = ism + 1
            endif
         enddo
      enddo

      average = sum/float(ism)
      write(6,77)average,k
 77   format('   JACOBI: average = ',f10.5,' at level ',i2)

!  set initial value
      init_val = 1.e-10
      do j = 1,jmt
         do i = 1,imt
            if(.not.maskold(i,j)) tnew(i,j) = init_val
         enddo
      enddo

      jj = jmt
 22   continue
      sum = 0.
      ism = 0
      do i = 1,imt
         if(maskold(i,jj)) then
            sum = sum + tnew(i,jj)
            ism = ism + 1
         endif
      enddo

      if(ism.eq.0) then
         jj = jj - 1
!        write(6,*)'**** NO LAND POINTS ****'
         go to 22
      endif
      average = sum/float(ism)

      jj = jmt
      do i = 1,imt
         if(.not.maskold(i,jj)) then
            tnew(i,jj) = init_val
            maskold(i,jj) = .true.
         endif
      enddo

      do iter = 1,niter

      if(mod(iter,50) .eq. 0) write(*,*)'iteration # ',iter

      epsmax = -1.e10

      do j = 1,jmt
         do i = 1,imt
            told(i,j) = tnew(i,j)
         enddo
      enddo

      do j = 1,jmt
         do i = 1,imt
            if(.not.maskold(i,j)) then
               sum = 0.

               if(i+1 .le. imt) then
                  sum = sum + told(i+1,j)
               else
                  sum = sum + told(3,j)
               endif

               if(i-1 .ge. 1) then
                  sum = sum + told(i-1,j)
               else
                  sum = sum + told(imt-2,j)
               endif

               if(j+1 .le. jmt) then
                  sum = sum + told(i,j+1)
               else
                  sum = sum + told(i,j)
               endif

               if(j-1 .ge. 1) then
                  sum = sum + told(i,j-1)
               else
                  sum = sum + told(i,j)
               endif

!-----------------------------------------------------------------------
!     fill in F at new ocean points with the smoothed field,
!     only fill in points which have neighboring old ocean points.
!-----------------------------------------------------------------------

               tnew(i,j) = 0.25*sum
               eps = abs((tnew(i,j) - told(i,j))/told(i,j))
               if(abs(tnew(i,j)).lt.0.1) eps = 0.
               epsmax = max(epsmax,eps)

            endif

         enddo
      enddo

      if(mod(iter,niter).eq.0) then
         print *, '       iteration ',iter,'; epsmax = ',epsmax
      endif

      enddo

      do j = 1,jmt
         do i = 1,imt
            t(i,j,k) = tnew(i,j)
         enddo
      enddo

      enddo

      return
      end