scale_meshutil_cubedsphere2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshutil_cubedsphere2d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_const, only: &
    pi => const_pi,      &
    eps => const_eps
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prc
 
  use scale_meshutil_2d, only: &
    meshutilcubedsphere2d_genrectdomain    => meshutil2d_genrectdomain,    &
    meshutilcubedsphere2d_genconnectivity  => meshutil2d_genconnectivity,  &
    meshutilcubedsphere2d_buildinteriormap => meshutil2d_buildinteriormap, &
    meshutilcubedsphere2d_genpatchboundarymap => meshutil2d_genpatchboundarymap
  !-----------------------------------------------------------------------------
  implicit none
  private
  
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: meshutilcubedsphere2d_genrectdomain
  public :: meshutilcubedsphere2d_genconnectivity  
  public :: meshutilcubedsphere2d_buildinteriormap
  public :: meshutilcubedsphere2d_buildglobalmap
  public :: meshutilcubedsphere2d_genpatchboundarymap
  public :: meshutilcubedsphere2d_modifyconnectivity
  public :: meshutilcubedsphere2d_getpanelconnectivity
  public :: meshutilcubedsphere2d_getpanelid
  
contains
!OCL SERIAL
  subroutine meshutilcubedsphere2d_buildglobalmap( &
    panelID_table, pi_table, pj_table,              &
    tileID_map, tileFaceID_map, tilePanelID_map,    &
    Ntile  )
    
    use scale_meshutil_2d, only: &
      meshutil2d_genconnectivity
    implicit none
 
    integer, intent(in) :: ntile
    integer, intent(out) :: panelid_table(ntile)
    integer, intent(out) :: pi_table(ntile)
    integer, intent(out) :: pj_table(ntile)
    integer, intent(out) :: tileid_map(4,ntile)
    integer, intent(out) :: tilefaceid_map(4,ntile)
    integer, intent(out) :: tilepanelid_map(4,ntile)
 
    integer :: ntileperpanel
    integer :: nex, ney, nvx, nvy
    integer, allocatable :: nodesid_2d(:,:,:)
    integer, allocatable :: etov(:,:)
    integer, allocatable :: etoe(:,:)
    integer, allocatable :: etof(:,:)
    integer :: i, j, k, f
    integer :: panelid
    integer :: tileid, tileid_r
    integer :: counter
    !-----------------------------------------------------------------------------
 
    ntileperpanel = ntile / 6
    ney = int( sqrt(dble(ntileperpanel)) )
    nex = ntileperpanel/ney
    nvx = nex + 1
    nvy = ney + 1
    allocate( nodesid_2d(nvx,nvy,6) )
    allocate( etov(ntile,4), etoe(ntile,4), etof(ntile,4) )
 
    counter = 0
    do panelid = 1, 6
    do j = 1, nvy
    do i = 1, nvx
      counter = counter + 1
      nodesid_2d(i,j,panelid) = counter
    end do
    end do
    end do
 
    !----
 
    tileid = 0
    do panelid = 1, 6
    do j = 1, ney
    do i = 1, nex
      tileid = tileid + 1
      panelid_table(tileid) = panelid
      pi_table(tileid) = i; pj_table(tileid) = j;
      etov(tileid,:) = (/ nodesid_2d(i,j  ,panelid), nodesid_2d(i+1,j  ,panelid),   &
                          nodesid_2d(i,j+1,panelid), nodesid_2d(i+1,j+1,panelid)   /)
    end do
    end do
    end do
 
    call meshutil2d_genconnectivity( etoe, etof, &
      etov, ntile, 4 )
    tileid_map(:,:) = transpose(etoe)
    tilefaceid_map(:,:) = transpose(etof)
    
    do tileid=1, ntile
    do f=1, 4
      tileid_r = tileid_map(f,tileid)
      tilepanelid_map(f,tileid) = panelid_table(tileid_r)
    end do
    end do
 
    call meshutilcubedsphere2d_modifyconnectivity( &
      tilepanelid_map, tileid_map, tilefaceid_map,          & ! (inout)
      panelid_table, pi_table, pj_table, nex, ney, ntile, 4 ) ! (in)
 
    return
  end subroutine meshutilcubedsphere2d_buildglobalmap
 
  !----
!OCL SERIAL
  subroutine meshutilcubedsphere2d_modifyconnectivity( tilePanelID_map, tileID_map, tileFaceID_map, &
    panelID_table, pi_table, pj_table, NeX, NeY, Ntile, Nface )
 
    integer, intent(in) :: ntile
    integer, intent(in) :: nface
    integer, intent(out) :: tileid_map(nface,ntile)
    integer, intent(out) :: tilefaceid_map(nface,ntile)
    integer, intent(out) :: tilepanelid_map(nface,ntile)
    integer, intent(in) :: panelid_table(ntile)
    integer, intent(in) :: pi_table(ntile)
    integer, intent(in) :: pj_table(ntile)
    integer, intent(in) :: nex, ney
 
    integer :: panel_connectivity(4,6)
    integer :: face_connectivity (4,6)
 
    integer :: tileid
    integer :: panelid
    integer :: f
    integer :: pi_, pj_    
    !-----------------------------------------------------------------------------
 
    call meshutilcubedsphere2d_getpanelconnectivity( &
      panel_connectivity, face_connectivity )
 
    
    do tileid=1, ntile
      panelid = panelid_table(tileid)
 
      do f=1, 4
        if ( tilefaceid_map(f,tileid) /= f ) cycle ! Does the face correspond the boundary of panel of cubed sphere?
 
        pi_ = pi_table(tileid)          
        pj_ = pj_table(tileid)  
        
        select case( panelid )
        case ( 1 )
          if ( mod(f,2) == 0 ) then ! West / East
            if ( pi_table(tileid) == 1   ) pi_ = nex
            if ( pi_table(tileid) == nex ) pi_ = 1
          else                      ! North / South
            if ( pj_table(tileid) == 1   ) pj_ = ney
            if ( pj_table(tileid) == ney ) pj_ = 1  
          end if
        case ( 2 )
          if ( mod(f,2) == 0 ) then ! West / East
            if ( pi_table(tileid) == 1   ) pi_ = nex
            if ( pi_table(tileid) == nex ) pi_ = 1
          else                      ! North / South
            if ( pj_table(tileid) == 1   ) then
              pi_ = nex; pj_ = ney - pi_table(tileid) + 1
            end if
            if ( pj_table(tileid) == ney ) then
              pi_ = nex ; pj_ = pi_table(tileid)
            end if
          end if
        case ( 3 )
          if ( mod(f,2) == 0 ) then ! West / East
            if ( pi_table(tileid) == 1   ) pi_ = nex
            if ( pi_table(tileid) == nex ) pi_ = 1
          else                      ! North / South
            if ( pj_table(tileid) == 1   ) then
              pi_ = nex - pi_table(tileid) + 1; pj_ = 1
            end if
            if ( pj_table(tileid) == ney ) then
              pi_ = nex - pi_table(tileid) + 1; pj_ = ney
            end if
          end if
        case ( 4 )
          if ( mod(f,2) == 0 ) then ! West / East
            if ( pi_table(tileid) == 1   ) pi_ = nex
            if ( pi_table(tileid) == nex ) pi_ = 1
          else                      ! North / South 
            if ( pj_table(tileid) == 1   ) then
              pi_ = 1; pj_ = pi_table(tileid)
            end if
            if ( pj_table(tileid) == ney ) then
              pi_ = 1; pj_ = ney - pi_table(tileid) + 1; 
            end if
          end if          
        case ( 5 )
          pj_ = ney
          if ( mod(f,2) == 0 ) then 
            if ( pi_table(tileid) == 1   ) pi_ = ney - pj_table(tileid) + 1  ! West            
            if ( pi_table(tileid) == nex ) pi_ = pj_table(tileid)            ! East
          else
            if ( pj_table(tileid) == 1   ) pi_ = pi_table(tileid)            ! South            
            if ( pj_table(tileid) == ney ) pi_ = nex - pi_table(tileid) + 1  ! North
          end if
        case ( 6 )
          pj_ = 1
          if ( mod(f,2) == 0 ) then 
            if ( pi_table(tileid) == 1   ) pi_ = pj_table(tileid)            ! West            
            if ( pi_table(tileid) == nex ) pi_ = ney - pj_table(tileid) + 1  ! East
          else
            if ( pj_table(tileid) == 1   ) pi_ = nex - pi_table(tileid) + 1  ! South            
            if ( pj_table(tileid) == ney ) pi_ = pi_table(tileid)            ! North
          end if
        end select
 
        tilepanelid_map(f,tileid) = panel_connectivity(f,panelid)
        tileid_map(f,tileid) = pi_ + (pj_ - 1) * nex + (tilepanelid_map(f,tileid) - 1) * nex * ney
        tilefaceid_map(f,tileid) = face_connectivity(f,panelid)
      end do ! loop for f
    end do ! loop for tile
 
    ! if (PRC_myrank==0) then
    !   write(*,*) " MeshUtilCubedSphere2D_modifyConnectivity"
    !   do tileID=1, Ntile
    !     write(*,'(a,i3,a,2i3)') "tileID=", tileID, ":", pi_table(tileID), pj_table(tileID)
    !     write(*,'(a,6i4)') "map_tile:", tileID_map(:,tileID)
    !     write(*,'(a,6i4)') "map_panel:", tilePanelID_map(:,tileID)
    !     write(*,'(a,6i4)') "mac_face:", tileFaceID_map(:,tileID)
    !   end do
    ! end if
 
    return
  end subroutine meshutilcubedsphere2d_modifyconnectivity
 
!OCL SERIAL
  subroutine meshutilcubedsphere2d_getpanelconnectivity( panel_connectivity, face_connectivity )
 
    implicit none
 
    integer, intent(out) :: panel_connectivity(4,6)
    integer, intent(out) :: face_connectivity(4,6)
 
    integer :: n
    integer :: zonal_panelid_list(0:5)
 
    zonal_panelid_list(:) = (/ 4, 1, 2, 3, 4, 1/)
    do n=1, 4
      panel_connectivity(1,n) = 6
      if (zonal_panelid_list(4-n) > 2) then
        face_connectivity(1,n)  = +zonal_panelid_list(4-n)
      else
        face_connectivity(1,n) =  -zonal_panelid_list(4-n)
      end if
 
      panel_connectivity(2,n) = zonal_panelid_list(n+1)
      face_connectivity(2,n)  = 4
 
      panel_connectivity(3,n) = 5
      face_connectivity(3,n)  = n
      if (n > 2) then
        face_connectivity(3,n)  = -n
      else
        face_connectivity(3,n) = n
      end if
 
      panel_connectivity(4,n) = zonal_panelid_list(n-1)
      face_connectivity(4,n)  = 2
    end do
 
    panel_connectivity(:,5) = (/ 1, 2, 3, 4 /)
    face_connectivity(:,5) = (/ 3, 3, -3, -3 /)
    panel_connectivity(:,6) = (/ 3, 2, 1, 4 /)
    face_connectivity(:,6) = (/ -1, -1, 1, 1 /)
 
  end subroutine meshutilcubedsphere2d_getpanelconnectivity
 
!OCL SERIAL
  subroutine meshutilcubedsphere2d_getpanelid( panelID, lon, lat, Np )
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_lonlat2cspos
    implicit none
 
    integer, intent(in) :: np
    integer, intent(out) :: panelid(np)
    real(rp), intent(in) :: lon(np)
    real(rp), intent(in) :: lat(np)
 
    integer :: p
    integer :: pnl
    real(rp) :: alph(np), beta(np)
    real(rp) :: lon_(np), lat_(np)
 
    real(rp), parameter :: eps = 1.0e-64_rp
    !------------------------------------------
 
    !$omp parallel do
    do p=1, np
      if ( abs(cos(lon(p))) < eps ) then
        lon_(p) = lon(p) + eps
      else
        lon_(p) = lon(p)
      end if
      if ( abs( lat(p) - 0.5_rp * pi ) < eps ) then
        lat_(p) = lat(p) - sign(eps, lat(p))
      else
        lat_(p) = lat(p)
      end if
    end do
 
    panelid(:) = -1
    do pnl=1, 6
      call cubedspherecoordcnv_lonlat2cspos( pnl, lon_, lat_, np, &
        alph, beta )
    
      select case(pnl)
      case (5)
        where ( lat(:) > 0.0_rp .and. abs(alph(:)) <= 0.25_rp * pi .and. abs(beta(:)) <= 0.25_rp * pi  )
          panelid(:) = pnl
        end where
      case (6)
        where ( lat(:) < 0.0_rp .and. abs(alph(:)) <= 0.25_rp * pi .and. abs(beta(:)) <= 0.25_rp * pi  )
          panelid(:) = pnl
        end where
      case default 
        where ( abs(alph(:)) <= 0.25_rp * pi .and. abs(beta(:)) <= 0.25_rp * pi  )
          panelid(:) = pnl
        end where
      end select     
    end do
 
    do p=1, np
      if (panelid(p) <  0) then
        log_error("MeshUtilCubedSphere2D_getPanelID",*) 'Fail to search a panel ID of cubed sphere grid!'
        write(*,*) "p=", p, ": (lon,lat)=", lon(p), lat(p), "(alpha,beta)=", alph(p), beta(p) 
        call prc_abort
      end if 
    end do
 
    return
  end subroutine meshutilcubedsphere2d_getpanelid
 
end module scale_meshutil_cubedsphere2d