!! $$  main program documentation block
!
! MODULE  :  SMOOTH
!   PRGMMR: BIJU THOMAS, URI/GSO, 2014-09-03 
!   LANG: FORTRAN 90/95
! 

MODULE smooth2d_md
  CONTAINS
  
  SUBROUTINE smooth2dave(maskval,msk,fld,nx,ny,npt)
  IMPLICIT NONE

   INTEGER,             INTENT(IN)           :: nx, ny, maskval
   INTEGER, DIMENSION(nx,ny), INTENT(IN)     :: msk
   INTEGER,                     INTENT(IN)   :: npt
   REAL, DIMENSION(nx,ny),   INTENT(INOUT)   :: fld
   REAL, DIMENSION(nx,ny)                    :: tmp
   INTEGER                                   :: i, j, i1, j1, n
   INTEGER                                   :: imn, imx, jmn, jmx

   DO j = 1, ny
   DO i = 1, nx 
      IF ( msk(i,j) /= maskval ) THEN
         n = 0
         tmp(i,j) = 0.0
         imn = max0(i - npt , 1)
         imx = min0(i + npt , nx)
         jmn = max0(j - npt , 1)
         jmx = min0(j + npt , ny)
         DO j1 = jmn, jmx
         DO i1 = imn, imx
            IF ( msk(i1, j1) /= maskval ) THEN
               tmp(i,j) = tmp(i,j) + fld(i1,j1)  
               n = n + 1
            END IF
         ENDDO
         ENDDO
         tmp(i,j) = tmp(i,j)/FLOAT(n)
      END IF
   ENDDO
   ENDDO

   DO j = 1, ny
   DO i = 1, nx
      IF ( msk(i,j) /= maskval) THEN
         fld(i,j) = tmp(i,j)
      END IF
   ENDDO
   ENDDO

  END SUBROUTINE smooth2dave

  SUBROUTINE smooth2d9ptwave(maskval,msk,fld,nx,ny)
  IMPLICIT NONE

   INTEGER,             INTENT(IN)           :: nx, ny, maskval
   INTEGER, DIMENSION(nx,ny), INTENT(IN)     :: msk
   INTEGER, PARAMETER                        :: npt = 1
   REAL, DIMENSION(nx,ny),   INTENT(INOUT)   :: fld
   REAL, DIMENSION(nx,ny)                    :: tmp
   INTEGER                                   :: i, j, i1, j1, n
   INTEGER                                   :: imn, imx, jmn, jmx
   REAL, PARAMETER                           :: w4 = 1./4., w8 = 1./8.
   REAL, PARAMETER                           :: w16 = 1./16.

! Spatial Layout of Weights  

! 1./16. 1./8.  1./16.
! 1./8.  1./4.  1./8.
! 1./16. 1./8.  1./16.

   DO j = 1, ny
   DO i = 1, nx
      IF ( msk(i,j) /= maskval ) THEN
         IF(fld(i,j) .gt. 200000.) PRINT*,n,i,j,fld(i,j) 
         n = 0
         tmp(i,j) = 0.0
         imn = max0(i - npt , 1)
         imx = min0(i + npt , nx)
         jmn = max0(j - npt , 1)
         jmx = min0(j + npt , ny)
         DO j1 = jmn, jmx
         DO i1 = imn, imx
            IF ( msk(i1, j1) /= maskval ) THEN
               tmp(i,j) = tmp(i,j) + fld(i1,j1)
               n = n + 1
!               if ( i==1003 .and. j == 28) then
!                 write(9,*)n,fld(i1,j1),tmp(i,j),i1,j1
!               endif
            END IF
         ENDDO
         ENDDO
         IF ( n == 9 ) THEN
            tmp(i,j) = w16*fld(i-1,j-1) + w8*fld(i-1,j) + w16*fld(i-1,j+1) &
                     + w8*fld(i,j-1)    + w4*fld(i,j)   + w8*fld(i,j+1)    &
                     + w16*fld(i+1,j-1) + w8*fld(i+1,j) + w16*fld(i+1,j+1)
            if ( i==1003 .and. j == 28) then
                 write(9,*)fld(i-1,j-1),fld(i-1,j),fld(i-1,j+1),fld(i,j-1), &
                 fld(i,j),fld(i,j+1),fld(i+1,j-1),fld(i+1,j),fld(i+1,j+1)
               endif

         ELSE
            tmp(i,j) = tmp(i,j)/FLOAT(n)
         END IF
      END IF
   ENDDO
   ENDDO

   DO j = 1, ny
   DO i = 1, nx
      IF ( msk(i,j) /= maskval) THEN
         fld(i,j) = tmp(i,j)
      END IF
   ENDDO
   ENDDO

  END SUBROUTINE smooth2d9ptwave

  SUBROUTINE smooth2d25ptwave(maskval,msk,fld,nx,ny)
  IMPLICIT NONE

   INTEGER,             INTENT(IN)           :: nx, ny, maskval
   INTEGER, DIMENSION(nx,ny), INTENT(IN)     :: msk
   INTEGER, PARAMETER                        :: npt = 2
   REAL, DIMENSION(nx,ny),   INTENT(INOUT)   :: fld
   REAL, DIMENSION(nx,ny)                    :: tmp
   INTEGER                                   :: i, j, i1, j1, n
   INTEGER                                   :: imn, imx, jmn, jmx
   REAL, PARAMETER                           :: w1 = 0.003, w2 = 0.013
   REAL, PARAMETER                           :: w3 = 0.022, w4 = 0.059
   REAL, PARAMETER                           :: w5 = 0.097, w6 = 0.159
 
! Spatial Layout of Weights

!  0.003 0.013 0.022 0.013 0.003
!  0.013 0.059 0.097 0.059 0.013
!  0.022 0.097 0.159 0.097 0.022
!  0.013 0.059 0.097 0.059 0.013
!  0.003 0.013 0.022 0.013 0.003 

   DO j = 1, ny
   DO i = 1, nx
      IF ( msk(i,j) /= maskval ) THEN
         n = 0
         tmp(i,j) = 0.0
         imn = max0(i - npt , 1)
         imx = min0(i + npt , nx)
         jmn = max0(j - npt , 1)
         jmx = min0(j + npt , ny)
         DO j1 = jmn, jmx
         DO i1 = imn, imx
            IF ( msk(i1, j1) /= maskval ) THEN
               tmp(i,j) = tmp(i,j) + fld(i1,j1)
               n = n + 1
            END IF
         ENDDO
         ENDDO
         IF ( n == 25 ) THEN
            tmp(i,j) = w1*fld(i-2,j-2) + w2*fld(i-2,j-1) + w3*fld(i-2,j) + &
                       w2*fld(i-2,j+1) + w1*fld(i-2,j+2) +                 &
                       w2*fld(i-1,j-2) + w4*fld(i-1,j-1) + w5*fld(i-1,j) + &
                       w4*fld(i-1,j+1) + w2*fld(i-1,j+2) +                 &
                       w3*fld(i,j-2)   + w5*fld(i,j-1)   + w6*fld(i,j)   + &
                       w5*fld(i,j+1)   + w3*fld(i,j+2)   +                 &
                       w2*fld(i+1,j-2) + w4*fld(i+1,j-1) + w5*fld(i+1,j) + &
                       w4*fld(i+1,j+1) + w2*fld(i+1,j+2) +                 &
                       w1*fld(i+2,j-2) + w2*fld(i+2,j-1) + w3*fld(i+2,j) + &
                       w2*fld(i+2,j+1) + w1*fld(i+2,j+2)
         ELSE
            tmp(i,j) = tmp(i,j)/FLOAT(n)
         END IF
      END IF
   ENDDO
   ENDDO

   DO j = 1, ny
   DO i = 1, nx
      IF ( msk(i,j) /= maskval) THEN
         fld(i,j) = tmp(i,j)
      END IF
   ENDDO
   ENDDO

  END SUBROUTINE smooth2d25ptwave

  SUBROUTINE smooth2d(msk,fld,nx,ny,coe)
! 
! Original (f77) version is from POM Code
! Algorithm by Dr. Ginis
! Abstract: Smoothing(coe>0) and de-smoothing(coe<0).
!
  IMPLICIT NONE
    
   INTEGER,             INTENT(IN)           :: nx, ny
   INTEGER, DIMENSION(nx,ny), INTENT(IN)     :: msk
   REAL,                     INTENT(IN)      :: coe
   REAL, DIMENSION(nx,ny),   INTENT(INOUT)   :: fld
   REAL, DIMENSION(nx,ny)                    :: tmp
   INTEGER                                   :: i, j

   DO j = 1, ny
   DO i = 2, nx -1
      IF ( msk(i,j) .ne. 0 )  THEN
         IF ( msk(i-1,j) .eq. 0 )  THEN  
            IF ( msk(i+1,j) .eq. 0 )  THEN
               tmp(i,j)  = fld(i,j)
            ELSE
               tmp(i,j)  = fld(i,j) + coe * ( fld(i+1,j) + fld(i,j) -   &
                                              2.0 * fld(i,j) )
            END IF
         ELSE
            IF ( msk(i+1,j) .eq. 0 )  THEN
               tmp(i,j)  = fld(i,j) + coe * ( fld(i,j) + fld(i-1,j) -   &
                                              2.0 * fld(i,j) )
            ELSE
               tmp(i,j)  = fld(i,j) + coe * ( fld(i+1,j) + fld(i-1,j) - &
                                              2.0 * fld(i,j) )
            END IF
         END IF
      END IF
   END DO
   END DO

   DO j = 1, ny
   DO i = 2, nx -1
     IF ( msk(i,j) .ne. 0 )  THEN
        fld(i,j) = tmp(i,j)
     END IF
   END DO
   END DO  

   DO j = 2, ny - 1 
   DO i = 1, nx 
      IF ( msk(i,j) .ne. 0 )  THEN
         IF ( msk(i,j-1) .eq. 0 )  THEN
            IF ( msk(i,j+1) .eq. 0 )  THEN
               tmp(i,j)  = fld(i,j)
            ELSE
               tmp(i,j)  = fld(i,j) + coe * ( fld(i,j+1) + fld(i,j) -   &
                                              2.0 * fld(i,j) )
            END IF
         ELSE
            IF ( msk(i,j+1) .eq. 0 )  THEN
               tmp(i,j)  = fld(i,j) + coe * ( fld(i,j) + fld(i,j-1) -   &
                                              2.0 * fld(i,j) )
            ELSE
               tmp(i,j)  = fld(i,j) + coe * ( fld(i,j+1) + fld(i,j-1) - &
                                              2.0 * fld(i,j) )
            END IF
         END IF
      END IF
   END DO
   END DO      

   DO j = 2, ny - 1
   DO i = 1, nx 
      IF ( msk(i,j) .ne. 0 )  THEN
         fld(i,j) = tmp(i,j)
      END IF
   END DO
   END DO

  END SUBROUTINE smooth2d
END MODULE smooth2d_md