soften Subroutine

   subroutine soften(tabin, mask, tabout, long, larg)
   !================================================================================================
   !< @note Function to smooth out a 2D array: each point is replaced by a weighted mean of its neighbors.
   !<
   !< \[
   !<   h_{i,j} = \frac{1}{16} \left( 4 h_{i, j} + 2 h_{i + 1, j    } + 2 h_{i - 1, j    } + 2 h_{i    , j + 1} + 2 h_{i    , j - 1}
   !<                                              + h_{i + 1, j - 1} +   h_{i - 1, j + 1} +   h_{i - 1, j + 1} +   h_{i + 1, j - 1} \right)
   !< \]
   !<
   !< @endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                                      :: long       !! *2D array width*
   integer(kind=I4), intent(in )                                      :: larg       !! *2D array height*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg)           :: tabin      !! *2D array in*
   real   (kind=R8), intent(out), dimension(1:long, 1:larg)           :: tabout     !! *2D array out*
   integer(kind=I4), intent(in ), dimension(1:long, 1:larg), optional :: mask       !! *mask*

      integer(kind=I4) :: i, j

      tabout(1:long, 1:larg) = tabin(1:long, 1:larg)

      if ( present(mask) ) then

         do j = 1 +1, larg -1
         do i = 1 +1, long -1

            if ( sum(mask(i-1:i+1, j-1:j+1)) < 9 ) then
               cycle
            else
               tabout(i, j) = ( 2*tabin(i, j) +tabin(i +1, j   ) +tabin(i -1, j   ) +                                          &  !
                                               tabin(i   , j +1) +tabin(i   , j -1) + ( tabin(i +1, j -1) +tabin(i -1, j -1) + &  !
                                                                                        tabin(i -1, j +1) +tabin(i +1, j +1) ) / 2._R8 ) / 8
            endif

         enddo
         enddo

      else

         do j = 1 +1, larg -1
         do i = 1 +1, long -1
            tabout(i, j) = ( 2*tabin(i, j) +tabin(i +1, j   ) +tabin(i -1, j   ) +                                          &  !
                                            tabin(i   , j +1) +tabin(i   , j -1) + ( tabin(i +1, j -1) +tabin(i -1, j -1) + &  !
                                                                                     tabin(i -1, j +1) +tabin(i +1, j +1) ) / 2._R8 ) / 8
         enddo
         enddo

      endif

   return
   endsubroutine soften