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
 
    limiter = present( limiter_sign )
 
    !$acc data copy(data) create(work1, work2)
    !$acc data copyin(limiter_sign) if(limiter)
 
    ! reduce grid-scale variation
    do ite = 1, nite
 
       call comm_vars8( data(:,:), 1 )
       call comm_wait ( data(:,:), 1, .false. )
 
       !$omp parallel do
       !$acc kernels
       do j = 1, ja
       do i = 1, ia
          work2(i,j) = data(i,j)
       end do
       end do
       !$acc end kernels
 
       p1 => work2
       p2 => work1
       do n = 1, order
          !$omp parallel do
          !$acc kernels
          !$acc loop collapse(2) independent
          do j = max(js,2), min(je,ja-1)
          do i = max(is,2), min(ie,ia-1)
             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
          !$acc end kernels
          if ( js == 1 ) then
             !$acc kernels
             !$acc loop independent
             do i = max(is,2), min(ie,ia-1)
                p2(i,js) = ( - p1(i+1,js) + p1(i,js)*2.0_rp - p1(i-1,js) &
                             - p1(i,js+1) + p1(i,js)                     ) / 6.0_rp
             end do
             !$acc end kernels
          else
             !$acc kernels
             !$acc loop independent
             do i = max(is,2), min(ie,ia-1)
                p2(i,js-1) = p2(i,js)
             end do
             !$acc end kernels
          end if
          if ( je == ja ) then
             !$acc kernels
             !$acc loop independent
             do i = max(is,2), min(ie,ia-1)
                p2(i,je) = ( - p1(i+1,je) + p1(i,je)*2.0_rp - p1(i-1,je) &
                                          + p1(i,je)        - p1(i,je-1) ) / 6.0_rp
             end do
             !$acc end kernels
          else
             !$acc kernels
             !$acc loop independent
             do i = max(is,2), min(ie,ia-1)
                p2(i,je+1) = p2(i,je)
             end do
             !$acc end kernels
          end if
          if ( is == 1 ) then
             !$acc kernels
             !$acc loop independent
             do j = max(js,2), min(je,ja-1)
                p2(is,j) = ( - p1(is+1,j) + p1(is,j)                     &
                             - p1(is,j+1) + p1(is,j)*2.0_rp - p1(is,j-1) ) / 6.0_rp
             end do
             !$acc end kernels
             if ( js == 1 ) then
                !$acc kernels
                p2(is,js) = ( - p1(is+1,js) + p1(is,js)               &
                              - p1(is,js+1) + p1(is,js)               ) / 4.0_rp
                !$acc end kernels
             end if
             if ( je == ja ) then
                !$acc kernels
                p2(is,je) = ( - p1(is+1,je) + p1(is,je)               &
                                            + p1(is,je) - p1(is,je-1) ) / 4.0_rp
                !$acc end kernels
             end if
          else
             !$acc kernels
             do j = max(js-1,1), min(je+1,ja)
                p2(is-1,j) = p2(is,j)
             end do
             !$acc end kernels
          end if
          if ( ie == ia ) then
             !$acc kernels
             !$acc loop independent
             do j = max(js,2), min(je,ja-1)
                p2(ie,j) = (              + p1(ie,j)        - p1(ie-1,j) &
                             - p1(ie,j+1) + p1(ie,j)*2.0_rp - p1(ie,j-1) ) / 6.0_rp
             end do
             !$acc end kernels
             if ( js == 1 ) then
                !$acc kernels
                p2(ie,js) = (               + p1(ie,js) - p1(ie-1,js) &
                              - p1(ie,js+1) + p1(ie,js)               ) / 4.0_rp
                !$acc end kernels
             end if
             if ( je == ja ) then
                !$acc kernels
                p2(ie,je) = (               + p1(ie,je) - p1(ie-1,je) &
                                            + p1(ie,je) - p1(ie,je-1) ) / 4.0_rp
                !$acc end kernels
             end if
          else
             !$acc kernels
             do j = max(js-1,1), min(je+1,ja)
                p2(ie+1,j) = p2(ie,j)
             end do
             !$acc end kernels
          end if
 
          call comm_vars8( p2(:,:), 1 )
          call comm_wait ( p2(:,:), 1, .false. )
 
          if ( mod(n,2) == 0 ) then
             p1 => work2
             p2 => work1
          else
             p1 => work1
             p2 => work2
          end if
       end do
 
       !$omp parallel do
       !$acc kernels
       do j = js, je
       do i = is, ie
          data(i,j) = data(i,j) - p1(i,j)
       end do
       end do
       !$acc end kernels
 
       if ( limiter ) then
          !$omp parallel do
          !$acc kernels
          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
          !$acc end kernels
       end if
    end do
 
    !$acc end data
    !$acc end data
 
    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
 
    !$acc data copy(data) create(work_data)
    !$acc data copyin(limiter_sign) create(work_sign) if(flag)
 
    do k = ks, ke
 
       !$omp parallel do
       !$acc kernels
       do j = js, je
       do i = is, ie
          work_data(i,j) = data(k,i,j)
       end do
       end do
       !$acc end kernels
       if ( flag ) then
          !$omp parallel do
          !$acc kernels
          do j = js, je
          do i = is, ie
             work_sign(i,j) = limiter_sign(k,i,j)
          end do
          end do
          !$acc end kernels
 
          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
       !$acc kernels
       do j = 1, ja
       do i = 1, ia
          data(k,i,j) = work_data(i,j)
       end do
       end do
       !$acc end kernels
 
    end do
 
    !$acc end data
    !$acc end data
 
    return
  end subroutine filter_hyperdiff_3d
 
end module scale_filter