scale_filter.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_filter
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: filter_hyperdiff

  interface filter_hyperdiff
     module procedure filter_hyperdiff_2d
     module procedure filter_hyperdiff_3d
  end interface filter_hyperdiff

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine filter_hyperdiff_2d( &
       IA, IS, IE, JA, JS, JE, &
       data, order, nite, &
       limiter_sign )
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    use scale_prc, only: &
       prc_abort
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je

    real(RP), intent(inout) :: data(ia,ja)

    integer,  intent(in)    :: order
    integer,  intent(in)    :: nite

    real(RP), intent(in), optional :: limiter_sign(ia,ja)

    real(RP), pointer :: p1(:,:)
    real(RP), pointer :: p2(:,:)
    real(RP), target  :: work1(ia,ja)
    real(RP), target  :: work2(ia,ja)

    logical :: limiter

    integer :: i, j
    integer :: ite, n

    call prof_rapstart('FILTER', 3)

    if ( mod(order,2) .ne. 0 ) then
       log_error("FILTER_hyperdiff_2D", *) "order must be even"
       call prc_abort
    end if

    if ( is < 2 ) then
       log_error("FILTER_hyperdiff_2D", *) "IS must be >= 2"
       call prc_abort
    end if
    if ( ie > ia-1  ) then
       log_error("FILTER_hyperdiff_2D", *) "IS must be <= IA-1"
       call prc_abort
    end if
    if ( js < 2 ) then
       log_error("FILTER_hyperdiff_2D", *) "JS must be >= 2"
       call prc_abort
    end if
    if ( je > ja-1  ) then
       log_error("FILTER_hyperdiff_2D", *) "JS must be <= JA-1"
       call prc_abort
    end if


    limiter = present( limiter_sign )

    ! reduce grid-scale variation
    do ite = 1, nite

       call comm_vars8( data(:,:), 1 )
       call comm_wait ( data(:,:), 1, .true. )

       !$omp parallel do
       do j = 1, ja
       do i = 1, ia
          work2(i,j) = data(i,j)
       end do
       end do

       p1 => work2
       p2 => work1
       do n = 1, order
          !$omp parallel do
          do j = js, je
          do i = is, ie
             p2(i,j) = ( - p1(i+1,j) + p1(i,j)*2.0_rp - p1(i-1,j) &
                         - p1(i,j+1) + p1(i,j)*2.0_rp - p1(i,j-1) ) / 8.0_rp
          end do
          end do

          call comm_vars8( p2(:,:), 1 )
          call comm_wait ( p2(:,:), 1, .true. )

          if ( mod(n,2) == 0 ) then
             p1 => work2
             p2 => work1
          else
             p1 => work1
             p2 => work2
          end if
       end do

       !$omp parallel do
       do j = js, je
       do i = is, ie
          data(i,j) = data(i,j) - p1(i,j)
       end do
       end do

       if ( limiter ) then
          !$omp parallel do
          do j = js, je
          do i = is, ie
             data(i,j) = sign( max( data(i,j) * limiter_sign(i,j), 0.0_rp ), limiter_sign(i,j) )
          end do
          end do
       end if
    end do

    call prof_rapend('FILTER', 3)

    return
  end subroutine filter_hyperdiff_2d

  !-----------------------------------------------------------------------------
  subroutine filter_hyperdiff_3d( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       data, order, nite, &
       limiter_sign )
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    use scale_prc, only: &
       prc_abort
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(inout) :: data(ka,ia,ja)

    integer,  intent(in)    :: order
    integer,  intent(in)    :: nite

    real(RP), intent(in), optional :: limiter_sign(ka,ia,ja)

    real(RP) :: work_data(ia,ja)
    real(RP) :: work_sign(ia,ja)

    logical :: flag

    integer :: k, i, j

    if ( present(limiter_sign) ) then
       flag = .true.
    else
       flag = .false.
    end if

    do k = ks, ke

       !$omp parallel do
       do j = js, je
       do i = is, ie
          work_data(i,j) = data(k,i,j)
       end do
       end do
       if ( flag ) then
          !$omp parallel do
          do j = js, je
          do i = is, ie
             work_sign(i,j) = limiter_sign(k,i,j)
          end do
          end do

          call filter_hyperdiff_2d( ia, is, ie, ja, js, je, &
                                    work_data(:,:), order, nite, &
                                    limiter_sign = work_sign )

       else

          call filter_hyperdiff_2d( ia, is, ie, ja, js, je, &
                                    work_data(:,:), order, nite )
       end if


       !$omp parallel do
       do j = 1, ja
       do i = 1, ia
          data(k,i,j) = work_data(i,j)
       end do
       end do

    end do


    return
  end subroutine filter_hyperdiff_3d

end module scale_filter