scale_file_common_meshfield.F90

Go to the documentation of this file.

!> module FElib / File / Common
!!
!! @par Description
!!           A common module for outputting field data 
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_file_common_meshfield
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
 
  use scale_element_base, only: elementbase1d, elementbase2d, elementbase3d
  use scale_mesh_base, only: meshdiminfo
  use scale_mesh_base1d, only: meshbase1d, &
    meshbase1d_dimtypeid_x, meshbase1d_dimtypeid_xt, &
    meshbase1d_dimtype_num 
  use scale_mesh_base2d, only: meshbase2d, &
    meshbase2d_dimtypeid_x, meshbase2d_dimtypeid_y,    &
    meshbase2d_dimtypeid_xy, meshbase2d_dimtypeid_xyt, &
    meshbase2d_dimtype_num 
  use scale_mesh_base3d, only: meshbase3d, &
    meshbase3d_dimtypeid_x, meshbase3d_dimtypeid_y, meshbase3d_dimtypeid_z,  &
    meshbase3d_dimtypeid_zt,                                                 &
    meshbase3d_dimtypeid_xy, meshbase3d_dimtypeid_xyt,                       &
    meshbase3d_dimtypeid_xyz, meshbase3d_dimtypeid_xyzt,                     &
    meshbase3d_dimtype_num
 
  use scale_mesh_rectdom2d, only: meshrectdom2d
  use scale_mesh_cubedom3d, only: meshcubedom3d
  use scale_mesh_cubedspheredom2d, only: meshcubedspheredom2d
  use scale_mesh_cubedspheredom3d, only: meshcubedspheredom3d
  use scale_localmesh_1d, only: localmesh1d
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
 
  use scale_meshfield_base, only: meshfield1d, meshfield2d, meshfield3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
 
  public :: file_common_meshfield_get_dims1d
  public :: file_common_meshfield_get_dims2d
  public :: file_common_meshfield_get_dims3d
  interface file_common_meshfield_get_dims
    module procedure file_common_meshfield_get_dims1d
    module procedure file_common_meshfield_get_dims2d
    module procedure file_common_meshfield_get_dims2d_cubedsphere
    module procedure file_common_meshfield_get_dims3d
    module procedure file_common_meshfield_get_dims3d_cubedsphere
  end interface
  public :: file_common_meshfield_get_dims
 
 
  public :: file_common_meshfield_get_axis1d
  public :: file_common_meshfield_get_axis2d
  public :: file_common_meshfield_get_axis3d
  interface file_common_meshfield_get_axis
    module procedure file_common_meshfield_get_axis1d
    module procedure file_common_meshfield_get_axis2d
    module procedure file_common_meshfield_get_axis2d_cubedsphere
    module procedure file_common_meshfield_get_axis3d
    module procedure file_common_meshfield_get_axis3d_cubedsphere
  end interface
  public :: file_common_meshfield_get_axis
 
  public :: file_common_meshfield_put_field1d_cartesbuf
  public :: file_common_meshfield_put_field2d_cartesbuf
  public :: file_common_meshfield_put_field2d_cubedsphere_cartesbuf  
  public :: file_common_meshfield_put_field3d_cartesbuf
  public :: file_common_meshfield_put_field3d_cubedsphere_cartesbuf  
 
  public :: file_common_meshfield_set_cartesbuf_field1d
  public :: file_common_meshfield_set_cartesbuf_field2d
  public :: file_common_meshfield_set_cartesbuf_field2d_cubedsphere
  public :: file_common_meshfield_set_cartesbuf_field3d
  public :: file_common_meshfield_set_cartesbuf_field3d_cubedsphere
 
  public :: file_common_meshfield_set_cartesbuf_field1d_local
  public :: file_common_meshfield_set_cartesbuf_field2d_local
  public :: file_common_meshfield_set_cartesbuf_field3d_local
 
  type, public :: file_common_meshfield_diminfo
    character(len=H_SHORT) :: type
    integer :: ndim
    character(len=H_SHORT) :: dims(3)
    character(len=H_SHORT) :: name
    character(len=H_MID) :: desc
    character(len=H_SHORT) :: unit
    integer :: count(3)
    integer :: size
    logical :: positive_down(3)
  end type file_common_meshfield_diminfo
 
  public :: file_common_meshfield_get_dtype
  
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
 
  !--------------------
  !
  !++ Private procedures
  !
  !-------------------
 
  private :: get_uniform_grid1d
  private :: set_dimension
 
contains
 
  !- 1D ---------------
 
!OCL SERIAL
  subroutine file_common_meshfield_get_dims1d( mesh1D, dim_name_postfix, & ! (in)
    dimsinfo ) ! (out)
    implicit none
    class(meshbase1d), target, intent(in) :: mesh1D
    character(len=H_SHORT), intent(in) :: dim_name_postfix
    type(file_common_meshfield_diminfo), intent(out) :: dimsinfo(MeshBase1D_DIMTYPE_NUM)
 
    integer :: i_size
    type(meshdiminfo), pointer :: diminfo
    type(meshdiminfo), pointer :: diminfo_x
    !-------------------------------------------------
 
    i_size = mesh1d%NeG * mesh1d%refElem1D%Np
 
    diminfo_x => mesh1d%dimInfo(meshbase1d_dimtypeid_x)
    call set_dimension( dimsinfo(meshbase1d_dimtypeid_x),    &
      diminfo_x, "X", 1, (/ diminfo_x%name /), (/ i_size /), &
      dim_name_postfix )
 
    diminfo => mesh1d%dimInfo(meshbase1d_dimtypeid_xt)
    call set_dimension( dimsinfo(meshbase1d_dimtypeid_xt),   &
      diminfo, "XT", 1, (/ diminfo_x%name /), (/ i_size  /), &
      dim_name_postfix )
 
    return
  end subroutine file_common_meshfield_get_dims1d
 
!OCL SERIAL  
  subroutine file_common_meshfield_get_axis1d( mesh1D, dimsinfo, x, &
    force_uniform_grid )
    implicit none
 
    class(meshbase1d), target, intent(in) :: mesh1D
    type(file_common_meshfield_diminfo), intent(in) :: dimsinfo(MeshBase1D_DIMTYPE_NUM)
    real(RP), intent(out) :: x(dimsinfo(MeshBase1D_DIMTYPEID_X)%size)
    logical, intent(in), optional :: force_uniform_grid
 
    integer :: n
    integer :: i
    integer :: is, ie
    type(elementbase1d), pointer :: refElem
    type(localmesh1d), pointer :: lcmesh
 
    logical :: uniform_grid = .false.
    real(RP), allocatable :: x_local(:)
    !-------------------------------------------------
    
    if ( present(force_uniform_grid) ) uniform_grid = force_uniform_grid
 
    do n=1,mesh1d%LOCAL_MESH_NUM
      lcmesh => mesh1d%lcmesh_list(n)
      refelem => lcmesh%refElem1D
 
      allocate( x_local(refelem%Np) )
 
      do i=1,mesh1d%NeG
        x_local(:) = lcmesh%pos_en(:,i,1)
        if ( uniform_grid ) call get_uniform_grid1d( x_local, refelem%Nfp )
 
        is = 1 + (i-1)*refelem%Np + (n-1)*refelem%Np*lcmesh%Ne
        ie = is + refelem%Np -1 
        x(is:ie) = x_local(:)
      end do
 
      deallocate( x_local )
    end do
    
    return
  end subroutine file_common_meshfield_get_axis1d
 
!OCL SERIAL
  subroutine file_common_meshfield_put_field1d_cartesbuf( mesh1D, field1D, &
    buf, force_uniform_grid )
    use scale_polynomial, only: &
      polynomial_genlegendrepoly_sub
    implicit none
    class(meshbase1d), target, intent(in) :: mesh1d
    class(meshfield1d), intent(in) :: field1d
    real(rp), intent(inout) :: buf(:)
    logical, intent(in), optional :: force_uniform_grid
 
    integer :: n, kelem1, p
    integer :: i, i2
    type(localmesh1d), pointer :: lcmesh
    type(elementbase1d), pointer :: refelem
    integer :: i0_s
 
    logical :: uniform_grid = .false.
    integer :: np
    real(rp), allocatable :: x_local(:)
    real(rp) :: x_local0, delx
    real(rp) :: ox(1)
    real(rp), allocatable :: spectral_coef(:)
    real(rp), allocatable :: p1d_ori_x(:,:)
    !------------------------------------------------
 
    if ( present(force_uniform_grid) ) uniform_grid = force_uniform_grid
 
    i0_s = 0
    do n=1, mesh1d%LOCAL_MESH_NUM
      lcmesh => mesh1d%lcmesh_list(n)
      refelem => lcmesh%refElem1D
      np = refelem%Np
      
      if ( uniform_grid ) then
        allocate( x_local(np) )
        allocate( spectral_coef(np) )
        allocate( p1d_ori_x(1,np) )
      end if
 
      do kelem1=lcmesh%NeS, lcmesh%NeE
        if ( uniform_grid ) then
          x_local(:) = lcmesh%pos_en(:,kelem1,1)
          x_local0 = x_local(1); delx = x_local(np) - x_local0
          call get_uniform_grid1d( x_local, np )
  
          spectral_coef(:) = matmul(refelem%invV(:,:), field1d%local(n)%val(:,kelem1))
          do i2=1, np
            ox = - 1.0_rp + 2.0_rp * (x_local(i2) - x_local0) / delx  
            call polynomial_genlegendrepoly_sub( refelem%PolyOrder, ox, p1d_ori_x(:,:) )
  
            i = i0_s + i2 + (kelem1-1)*np 
            buf(i) = 0.0_rp 
            do p=1, np
              buf(i) = buf(i) + &
                  p1d_ori_x(1,p) * sqrt( dble(p-1) + 0.5_rp ) * spectral_coef(p)
            end do
          end do
        else
          do i2=1, np
            i = i0_s + i2 + (kelem1-1)*np 
            buf(i) = field1d%local(n)%val(i2,kelem1)
          end do
        end if
      end do
 
      i0_s = i0_s + lcmesh%Ne * refelem%Np
      if ( uniform_grid ) then
        deallocate( x_local )
        deallocate( spectral_coef )
        deallocate( p1d_ori_x )
      end if
    end do
 
    return
  end subroutine file_common_meshfield_put_field1d_cartesbuf
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field1d( mesh1D, buf, &
    field1D )
    implicit none
    class(meshbase1d), target, intent(in) :: mesh1d
    real(rp), intent(in) :: buf(:)
    class(meshfield1d), intent(inout) :: field1d
 
    integer :: n
    integer :: i0
    type(localmesh1d), pointer :: lcmesh
    type(elementbase1d), pointer :: refelem
    integer :: i0_s
    !----------------------------------------------------
 
    i0_s = 0
 
    do i0=1, mesh1d%LOCAL_MESH_NUM
      n = i0
      lcmesh => mesh1d%lcmesh_list(n)
      refelem => lcmesh%refElem1D
 
      call file_common_meshfield_set_cartesbuf_field1d_local(  &
        lcmesh, buf(:), i0_s,                                  &
        field1d%local(n)%val(:,:)                              )
 
      i0_s = i0_s + lcmesh%Ne * refelem%Np
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field1d
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field1d_local( &
    lcmesh, buf, i0_s,                                          &
    val )
    implicit none
    type(localmesh1d), intent(in) :: lcmesh
    real(rp), intent(in) :: buf(:)
    integer, intent(in) :: i0_s
    real(rp), intent(inout) :: val(lcmesh%refelem1d%np,lcmesh%nea)
 
    integer :: kelem1
    integer :: i, i1, i2
    type(elementbase1d), pointer :: refelem
    integer :: indx
    !----------------------------------------------------
 
    refelem => lcmesh%refElem1D
 
    do i1=1, lcmesh%Ne
      kelem1 = i1
      do i2=1, refelem%Np
        i = i0_s + i2 + (i1-1)*refelem%Np
        indx = i2
        val(indx,kelem1) = buf(i)
      end do
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field1d_local
 
  !- 2D ---------------
 
!OCL SERIAL
  subroutine file_common_meshfield_get_dims2d( mesh2D, dim_name_postfix, & ! (in)
    dimsinfo ) ! (out)
    implicit none
 
    class(meshrectdom2d), target, intent(in) :: mesh2D
    character(len=H_SHORT), intent(in) :: dim_name_postfix
    type(file_common_meshfield_diminfo), intent(out) :: dimsinfo(MeshBase2D_DIMTYPE_NUM)
 
    type(localmesh2d), pointer :: lcmesh
    integer :: i, j, n
 
    integer :: i_size, j_size
    type(meshdiminfo), pointer :: diminfo
    type(meshdiminfo), pointer :: diminfo_x
    type(meshdiminfo), pointer :: diminfo_y
    !-------------------------------------------------
    
    i_size = 0
    do i=1, size(mesh2d%rcdomIJ2LCMeshID,1)
      n = mesh2d%rcdomIJ2LCMeshID(i,1)
      lcmesh => mesh2d%lcmesh_list(n)
      i_size =i_size + lcmesh%NeX * lcmesh%refElem2D%Nfp
    end do
 
    j_size = 0
    do j=1, size(mesh2d%rcdomIJ2LCMeshID,2)
      n = mesh2d%rcdomIJ2LCMeshID(1,j)
      lcmesh => mesh2d%lcmesh_list(n)
      j_size = j_size + lcmesh%NeY * lcmesh%refElem2D%Nfp
    end do
  
    diminfo_x => mesh2d%dimInfo(meshbase2d_dimtypeid_x)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_x),    &
      diminfo_x, "X", 1, (/ diminfo_x%name /), (/ i_size /), &
      dim_name_postfix )
 
    diminfo_y => mesh2d%dimInfo(meshbase2d_dimtypeid_y)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_y),    &
      diminfo_y, "Y", 1, (/ diminfo_y%name /), (/ j_size /), &
      dim_name_postfix )
 
    diminfo => mesh2d%dimInfo(meshbase2d_dimtypeid_xy)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_xy),    &
      diminfo, "XY", 2, (/ diminfo_x%name, diminfo_y%name /), &
      (/ i_size, j_size /), dim_name_postfix )
 
    diminfo => mesh2d%dimInfo(meshbase2d_dimtypeid_xyt)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_xyt),    &
      diminfo, "XYT", 2, (/ diminfo_x%name, diminfo_y%name /), &
      (/ i_size, j_size /), dim_name_postfix )
  
    return
  end subroutine file_common_meshfield_get_dims2d
 
!OCL SERIAL
  subroutine file_common_meshfield_get_dims2d_cubedsphere( mesh2D, dim_name_postfix, & ! (in)
    dimsinfo ) ! (out)
    implicit none
 
    class(meshcubedspheredom2d), target, intent(in) :: mesh2D
    character(len=H_SHORT), intent(in) :: dim_name_postfix
    type(file_common_meshfield_diminfo), intent(out) :: dimsinfo(MeshBase2D_DIMTYPE_NUM)
 
    type(localmesh2d), pointer :: lcmesh
    integer :: i, j, n
 
    integer :: i_size, j_size
 
    type(meshdiminfo), pointer :: diminfo
    type(meshdiminfo), pointer :: diminfo_x
    type(meshdiminfo), pointer :: diminfo_y
    !-------------------------------------------------
    
    i_size = 0
    do i=1, size(mesh2d%rcdomIJP2LCMeshID,1)
      n = mesh2d%rcdomIJP2LCMeshID(i,1,1)
      lcmesh => mesh2d%lcmesh_list(n)
      i_size =i_size + lcmesh%NeX * lcmesh%refElem2D%Nfp
    end do
 
    j_size = 0
    do j=1, size(mesh2d%rcdomIJP2LCMeshID,2)
      n = mesh2d%rcdomIJP2LCMeshID(1,j,1)
      lcmesh => mesh2d%lcmesh_list(n)
      j_size = j_size + lcmesh%NeY * lcmesh%refElem2D%Nfp
    end do
 
    j_size = j_size * size(mesh2d%rcdomIJP2LCMeshID,3)
  
    diminfo_x => mesh2d%dimInfo(meshbase2d_dimtypeid_x)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_x),    &
      diminfo_x, "X", 1, (/ diminfo_x%name /), (/ i_size /), &
      dim_name_postfix )
 
    diminfo_y => mesh2d%dimInfo(meshbase2d_dimtypeid_y)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_y),   &
      diminfo_y, "Y", 1, (/ diminfo_y%name /), (/ j_size /), &
      dim_name_postfix )
 
    diminfo => mesh2d%dimInfo(meshbase2d_dimtypeid_xy)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_xy),    &
      diminfo, "XY", 2, (/ diminfo_x%name, diminfo_y%name /), &
      (/ i_size, j_size /), dim_name_postfix )
 
    diminfo => mesh2d%dimInfo(meshbase2d_dimtypeid_xyt)
    call set_dimension( dimsinfo(meshbase2d_dimtypeid_xyt),    &
      diminfo, "XYT", 2, (/ diminfo_x%name, diminfo_y%name /), &
      (/ i_size, j_size /), dim_name_postfix )
 
    return
  end subroutine file_common_meshfield_get_dims2d_cubedsphere
 
!OCL SERIAL
  subroutine file_common_meshfield_get_axis2d( mesh2D, dimsinfo, x, y, &
    force_uniform_grid  )
    implicit none
 
    class(meshrectdom2d), target, intent(in) :: mesh2D  
    type(file_common_meshfield_diminfo), intent(in) :: dimsinfo(MeshBase2D_DIMTYPE_NUM)
    real(RP), intent(out) :: x(dimsinfo(MeshBase2D_DIMTYPEID_X)%size)
    real(RP), intent(out) :: y(dimsinfo(MeshBase2D_DIMTYPEID_Y)%size)
    logical, intent(in), optional :: force_uniform_grid
 
    integer :: n
    integer :: ni, nj
    integer :: k
    integer :: i, j 
    type(elementbase2d), pointer :: refElem
    type(localmesh2d), pointer :: lcmesh
 
    integer :: is, js, ie, je, igs, jgs
 
    logical :: uniform_grid = .false.
    real(RP), allocatable :: x_local(:)
    real(RP), allocatable :: y_local(:)
    !-------------------------------------------------
    
    if ( present(force_uniform_grid) ) uniform_grid = force_uniform_grid
 
    igs = 0; jgs = 0
    do nj=1, size(mesh2d%rcdomIJ2LCMeshID,2) 
    do ni=1, size(mesh2d%rcdomIJ2LCMeshID,1)
      n = mesh2d%rcdomIJ2LCMeshID(ni,nj)
      lcmesh => mesh2d%lcmesh_list(n)
      refelem => lcmesh%refElem2D
 
      allocate( x_local(refelem%Nfp), y_local(refelem%Nfp) )
 
      do j=1, lcmesh%NeY
      do i=1, lcmesh%NeX
        k = i + (j-1) * lcmesh%NeX
        if ( j==1 .and. nj == 1 ) then
          x_local(:) = lcmesh%pos_en(refelem%Fmask(:,1),k,1)
          if ( uniform_grid ) call get_uniform_grid1d( x_local, refelem%Nfp )
 
          is = igs + 1 + (i-1)*refelem%Nfp
          ie = is + refelem%Nfp - 1
          x(is:ie) = x_local(:)
        end if
        if ( i==1 .and. ni == 1 ) then
          y_local(:) = lcmesh%pos_en(refelem%Fmask(:,4),k,2)
          if ( uniform_grid ) call get_uniform_grid1d( y_local, refelem%Nfp )
 
          js = jgs + 1 + (j-1)*refelem%Nfp
          je = js + refelem%Nfp - 1
          y(js:je) = y_local(:)
        end if
      end do
      end do
 
      igs = ie; jgs = je
      deallocate( x_local, y_local )
    end do
    end do
 
    return
  end subroutine file_common_meshfield_get_axis2d
 
!OCL SERIAL
  subroutine file_common_meshfield_get_axis2d_cubedsphere( mesh2D, dimsinfo, x, y )
    use scale_const, only: &
      pi => const_pi
    implicit none
 
    class(meshcubedspheredom2d), target, intent(in) :: mesh2D  
    type(file_common_meshfield_diminfo), intent(in) :: dimsinfo(MeshBase2D_DIMTYPE_NUM)
    real(RP), intent(out) :: x(dimsinfo(MeshBase2D_DIMTYPEID_X)%size)
    real(RP), intent(out) :: y(dimsinfo(MeshBase2D_DIMTYPEID_Y)%size)
 
    integer :: ni, nj, np, n
    integer :: k
    integer :: i, j
    type(elementbase2d), pointer :: refElem
    type(localmesh2d), pointer :: lcmesh
 
    integer :: is, js, ie, je, igs, jgs
 
    logical :: uniform_grid = .false.
    real(RP), allocatable :: x_local(:)
    real(RP), allocatable :: y_local(:)
    !-------------------------------------------------
    
    igs = 0; jgs = 0
    
    do np=1, size(mesh2d%rcdomIJP2LCMeshID,3) 
    do nj=1, size(mesh2d%rcdomIJP2LCMeshID,2) 
    do ni=1, size(mesh2d%rcdomIJP2LCMeshID,1)
      n = mesh2d%rcdomIJP2LCMeshID(ni,nj,np)
      lcmesh => mesh2d%lcmesh_list(n)
      refelem => lcmesh%refElem2D
 
      allocate( x_local(refelem%Nfp), y_local(refelem%Nfp) )
 
      do j=1, lcmesh%NeY
      do i=1, lcmesh%NeX
        k = i + (j-1) * lcmesh%NeX
        if ( j==1 .and. nj == 1 .and. np == 1) then
          x_local(:) = lcmesh%pos_en(refelem%Fmask(:,1),k,1)
 
          is = igs + 1 + (i-1)*refelem%Nfp
          ie = is + refelem%Nfp - 1
          x(is:ie) = x_local(:)
        end if
        if ( i==1 .and. ni == 1 ) then
          y_local(:) = lcmesh%pos_en(refelem%Fmask(:,4),k,2) &
                     + ( lcmesh%panelID - 1.0_rp ) * 0.5_rp * pi
 
          js = jgs + 1 + (j-1)*refelem%Nfp
          je = js + refelem%Nfp - 1
          y(js:je) = y_local(:)
        end if
      end do
      end do
 
      igs = ie; jgs = je
      deallocate( x_local, y_local )
    end do
    end do
    end do
 
    return
  end subroutine file_common_meshfield_get_axis2d_cubedsphere
 
!OCL SERIAL
  subroutine file_common_meshfield_put_field2d_cartesbuf( mesh2D, field2D, &
    buf, force_uniform_grid )
    use scale_polynomial, only: &
      polynomial_genlegendrepoly_sub
    implicit none
    class(meshrectdom2d), target, intent(in) :: mesh2d
    class(meshfield2d), intent(in) :: field2d
    real(rp), intent(inout) :: buf(:,:)
    logical, intent(in), optional :: force_uniform_grid
 
    integer :: n, kelem1
    integer :: i0, j0, i1, j1, i2, j2, i, j
    type(localmesh2d), pointer :: lcmesh
    type(elementbase2d), pointer :: refelem
    integer :: i0_s, j0_s
 
    logical :: uniform_grid = .false.
    integer :: nfp
    real(rp), allocatable :: x_local(:)
    real(rp) :: x_local0, delx
    real(rp), allocatable :: y_local(:)
    real(rp) :: y_local0, dely
    real(rp) :: ox(1), oy(1)
    real(rp), allocatable :: spectral_coef(:)
    real(rp), allocatable :: p1d_ori_x(:,:)
    real(rp), allocatable :: p1d_ori_y(:,:)
    integer :: l, p1, p2
    !------------------------------------------------
 
    if ( present(force_uniform_grid) ) uniform_grid = force_uniform_grid
 
    i0_s = 0; j0_s = 0
    
    do j0=1, size(mesh2d%rcdomIJ2LCMeshID,2)
      do i0=1, size(mesh2d%rcdomIJ2LCMeshID,1)
        n =  mesh2d%rcdomIJ2LCMeshID(i0,j0)
 
        lcmesh => mesh2d%lcmesh_list(n)
        refelem => lcmesh%refElem2D
        nfp = refelem%Nfp
        
        if ( uniform_grid ) then
          allocate( x_local(nfp), y_local(nfp) )
          allocate( spectral_coef(refelem%Np) )
          allocate( p1d_ori_x(1,nfp), p1d_ori_y(1,nfp) )
        end if
    
        do j1=1, lcmesh%NeY
        do i1=1, lcmesh%NeX
          kelem1 = i1 + (j1-1)*lcmesh%NeX
 
          if ( uniform_grid ) then
            x_local(:) = lcmesh%pos_en(refelem%Fmask(1:nfp,1),kelem1,1)
            x_local0 = x_local(1); delx = x_local(nfp) - x_local0
            y_local(:) = lcmesh%pos_en(refelem%Fmask(1:nfp,4),kelem1,2)
            y_local0 = y_local(1); dely = y_local(nfp) - y_local0
            call get_uniform_grid1d( x_local, nfp )
            call get_uniform_grid1d( y_local, nfp )
    
            spectral_coef(:) = matmul(refelem%invV(:,:), field2d%local(n)%val(:,kelem1))
            do j2=1, nfp
            do i2=1, nfp
              ox(1) = - 1.0_rp + 2.0_rp * (x_local(i2) - x_local0) / delx
              oy(1) = - 1.0_rp + 2.0_rp * (y_local(j2) - y_local0) / dely
    
              call polynomial_genlegendrepoly_sub( refelem%PolyOrder, ox, p1d_ori_x(:,:) ) 
              call polynomial_genlegendrepoly_sub( refelem%PolyOrder, oy, p1d_ori_y(:,:) ) 
    
              i = i0_s + i2 + (i1-1)*nfp
              j = j0_s + j2 + (j1-1)*nfp
              buf(i,j) = 0.0_rp 
              do p2=1, nfp
              do p1=1, nfp
                l = p1 + (p2-1)*nfp
                buf(i,j) = buf(i,j) + &
                    ( p1d_ori_x(1,p1) * p1d_ori_y(1,p2) )             &
                  * sqrt((dble(p1-1) + 0.5_rp)*(dble(p2-1) + 0.5_rp)) &
                  * spectral_coef(l)
              end do
              end do
            end do
            end do
          
          else
 
            do j2=1, nfp
            do i2=1, nfp
              i = i0_s + i2 + (i1-1)*nfp 
              j = j0_s + j2 + (j1-1)*nfp
              buf(i,j) = field2d%local(n)%val(i2+(j2-1)*nfp,kelem1)
            end do
            end do
          
          end if
        end do
        end do
        
        if ( uniform_grid ) then
          deallocate( x_local, y_local  )
          deallocate( spectral_coef )
          deallocate( p1d_ori_x, p1d_ori_y )
        end if
 
        i0_s = i0_s + lcmesh%NeX * refelem%Nfp
      end do
      j0_s = j0_s + lcmesh%NeY * refelem%Nfp
    end do
 
    return
  end subroutine file_common_meshfield_put_field2d_cartesbuf
 
!OCL SERIAL
  subroutine file_common_meshfield_put_field2d_cubedsphere_cartesbuf( mesh2D, field2D, &
    buf )
    use scale_polynomial, only: &
      polynomial_genlegendrepoly
    implicit none
    class(meshcubedspheredom2d), target, intent(in) :: mesh2d
    class(meshfield2d), intent(in) :: field2d
    real(rp), intent(inout) :: buf(:,:)
 
    integer :: n, kelem1
    integer :: i0, j0, p0, i1, j1, i2, j2, i, j
    type(localmesh2d), pointer :: lcmesh
    type(elementbase2d), pointer :: refelem
    integer :: i0_s, j0_s
    integer :: nfp
    !------------------------------------------------
 
    i0_s = 0; j0_s = 0
 
    do p0=1, size(mesh2d%rcdomIJP2LCMeshID,3)    
      do j0=1, size(mesh2d%rcdomIJP2LCMeshID,2)
        do i0=1, size(mesh2d%rcdomIJP2LCMeshID,1)
          n = mesh2d%rcdomIJP2LCMeshID(i0,j0,p0)
 
          lcmesh => mesh2d%lcmesh_list(n)
          refelem => lcmesh%refElem2D
          nfp = refelem%Nfp
              
          do j1=1, lcmesh%NeY
          do i1=1, lcmesh%NeX
            kelem1 = i1 + (j1-1)*lcmesh%NeX
 
            do j2=1, nfp
            do i2=1, nfp
              i = i0_s + i2 + (i1-1)*nfp
              j = j0_s + j2 + (j1-1)*nfp
              buf(i,j) = field2d%local(n)%val(i2+(j2-1)*nfp,kelem1)
            end do
            end do
            
          end do
          end do
          
          i0_s = i0_s + lcmesh%NeX * refelem%Nfp
        end do
        j0_s = j0_s + lcmesh%NeY * refelem%Nfp
      end do
      i0_s = 0
    end do
    
    return
  end subroutine file_common_meshfield_put_field2d_cubedsphere_cartesbuf
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field2d( mesh2D, buf, &
    field2D )
    implicit none
    class(meshrectdom2d), target, intent(in) :: mesh2d
    real(rp), intent(in) :: buf(:,:)
    class(meshfield2d), intent(inout) :: field2d
 
    integer :: n
    integer :: i0, j0
    type(localmesh2d), pointer :: lcmesh
    type(elementbase2d), pointer :: refelem
    integer :: i0_s, j0_s
    !----------------------------------------------------
 
    i0_s = 0; j0_s = 0
 
    do j0=1, size(mesh2d%rcdomIJ2LCMeshID,2)
      do i0=1, size(mesh2d%rcdomIJ2LCMeshID,1)
        n = mesh2d%rcdomIJ2LCMeshID(i0,j0)
        lcmesh => mesh2d%lcmesh_list(n)
        refelem => lcmesh%refElem2D
 
        call file_common_meshfield_set_cartesbuf_field2d_local(  &
          lcmesh, buf(:,:), i0_s, j0_s,                          &
          field2d%local(n)%val(:,:)                              )
        !$acc update device(field2d%local(n)%val)
 
        i0_s = i0_s + lcmesh%NeX * refelem%Nfp
      end do
      j0_s = j0_s + lcmesh%NeY * refelem%Nfp
      i0_s = 0
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field2d
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field2d_local( &
    lcmesh, buf, i0_s, j0_s,                                    &
    val )
    implicit none
    type(localmesh2d), intent(in) :: lcmesh
    real(rp), intent(in) :: buf(:,:)
    integer, intent(in) :: i0_s, j0_s
    real(rp), intent(inout) :: val(lcmesh%refelem2d%np,lcmesh%nea)
 
    integer :: kelem1
    integer :: i1, j1, i2, j2, i, j
    type(elementbase2d), pointer :: refelem
    integer :: indx
    !----------------------------------------------------
 
    refelem => lcmesh%refElem2D
 
    do j1=1, lcmesh%NeY
    do i1=1, lcmesh%NeX
      kelem1 = i1 + (j1-1)*lcmesh%NeX
      do j2=1, refelem%Nfp
      do i2=1, refelem%Nfp
        i = i0_s + i2 + (i1-1)*refelem%Nfp
        j = j0_s + j2 + (j1-1)*refelem%Nfp
        indx = i2 + (j2-1)*refelem%Nfp
        val(indx,kelem1) = buf(i,j)
      end do
      end do
    end do
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field2d_local
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field2d_cubedsphere( mesh2D, buf, &
    field2D )
    implicit none
    class(meshcubedspheredom2d), target, intent(in) :: mesh2d
    real(rp), intent(in) :: buf(:,:)
    class(meshfield2d), intent(inout) :: field2d
 
    integer :: n
    integer :: i0, j0, p0
    type(localmesh2d), pointer :: lcmesh
    type(elementbase2d), pointer :: refelem
    integer :: i0_s, j0_s, p0_s
    !----------------------------------------------------
 
    i0_s = 0; j0_s = 0; p0_s = 0
 
    do p0=1, size(mesh2d%rcdomIJP2LCMeshID,3)
      do j0=1, size(mesh2d%rcdomIJP2LCMeshID,2)
        do i0=1, size(mesh2d%rcdomIJP2LCMeshID,1)
          n = mesh2d%rcdomIJP2LCMeshID(i0,j0,p0)
          lcmesh => mesh2d%lcmesh_list(n)
          refelem => lcmesh%refElem2D
 
          call file_common_meshfield_set_cartesbuf_field2d_local(  &
            lcmesh, buf(:,:), i0_s, j0_s,                          &
            field2d%local(n)%val(:,:)                              )
          !$acc update device(field2d%local(n)%val)
 
          i0_s = i0_s + lcmesh%NeX * refelem%Nfp
        end do
        j0_s = j0_s + lcmesh%NeY * refelem%Nfp
        i0_s = 0
      end do
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field2d_cubedsphere
 
  !- 3D ------------
 
!OCL SERIAL
  subroutine file_common_meshfield_get_dims3d( mesh3D, dim_name_postfix, & ! (in)
    dimsinfo ) ! (out)
    implicit none
 
    class(meshcubedom3d), target, intent(in) :: mesh3D
    character(len=H_SHORT), intent(in) :: dim_name_postfix
    type(file_common_meshfield_diminfo), intent(out) :: dimsinfo(MESHBASE3D_DIMTYPE_NUM)
 
    type(localmesh3d), pointer :: lcmesh
    integer :: i, j, k, n
    integer :: i_size, j_size, k_size
 
    type(meshdiminfo), pointer :: dimInfo
    type(meshdiminfo), pointer :: dimInfo_x
    type(meshdiminfo), pointer :: dimInfo_y
    type(meshdiminfo), pointer :: dimInfo_z
    !-------------------------------------------------
    
    i_size = 0
    do i=1, size(mesh3d%rcdomIJK2LCMeshID,1)
      n = mesh3d%rcdomIJK2LCMeshID(i,1,1)
      lcmesh => mesh3d%lcmesh_list(n)
      i_size = i_size + lcmesh%NeX * lcmesh%refElem3D%Nnode_h1D
    end do
 
    j_size = 0
    do j=1, size(mesh3d%rcdomIJK2LCMeshID,2)
      n = mesh3d%rcdomIJK2LCMeshID(1,j,1)
      lcmesh => mesh3d%lcmesh_list(n)
      j_size = j_size + lcmesh%NeY * lcmesh%refElem3D%Nnode_h1D
    end do
 
    k_size = 0
    do k=1, size(mesh3d%rcdomIJK2LCMeshID,3)
      n = mesh3d%rcdomIJK2LCMeshID(1,1,k)
      lcmesh => mesh3d%lcmesh_list(n)
      k_size = k_size + lcmesh%NeZ * lcmesh%refElem3D%Nnode_v
    end do  
 
    diminfo_x => mesh3d%dimInfo(meshbase3d_dimtypeid_x)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_x),    &
      diminfo_x, "X", 1, (/ diminfo_x%name /), (/ i_size /), &
      dim_name_postfix )
 
    diminfo_y => mesh3d%dimInfo(meshbase3d_dimtypeid_y)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_y),    &
      diminfo_y, "Y", 1, (/ diminfo_y%name /), (/ j_size /), &
      dim_name_postfix )
 
    diminfo_z => mesh3d%dimInfo(meshbase3d_dimtypeid_z)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_z),           &
      diminfo_z, "Z", 1, (/ diminfo_z%name /), (/ k_size /),        &
      dim_name_postfix, positive_down=(/ diminfo_z%positive_down /) )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_zt)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_zt),  &
      diminfo, "ZT", 1, (/ diminfo_z%name /), (/ k_size /), &
      dim_name_postfix )
 
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xy)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xy),     &
      diminfo, "XY", 2, (/ diminfo_x%name, diminfo_y%name  /), &
      (/ i_size, j_size /), dim_name_postfix )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xyt)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xyt),     &
      diminfo, "XY", 2, (/ diminfo_x%name, diminfo_y%name  /),  &
      (/ i_size, j_size /), dim_name_postfix )
    
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xyz)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xyz),                    &
      diminfo, "XYZ", 3, (/ diminfo_x%name, diminfo_y%name, diminfo_z%name /), &
      (/ i_size, j_size, k_size /), dim_name_postfix )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xyzt)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xyzt),                    &
      diminfo, "XYZT", 3, (/ diminfo_x%name, diminfo_y%name, diminfo_z%name /), &
      (/ i_size, j_size, k_size /), dim_name_postfix )
 
    return
  end subroutine file_common_meshfield_get_dims3d
 
!OCL SERIAL
  subroutine file_common_meshfield_get_dims3d_cubedsphere( mesh3D, dim_name_postfix, & ! (in)
    dimsinfo ) ! (out)
    implicit none
 
    class(meshcubedspheredom3d), target, intent(in) :: mesh3D
    character(len=H_SHORT), intent(in) :: dim_name_postfix
    type(file_common_meshfield_diminfo), intent(out) :: dimsinfo(MESHBASE3D_DIMTYPE_NUM)
 
    type(localmesh3d), pointer :: lcmesh
    integer :: i, j, k, n
 
    integer :: i_size, j_size, k_size
 
    type(meshdiminfo), pointer :: diminfo
    type(meshdiminfo), pointer :: diminfo_x
    type(meshdiminfo), pointer :: diminfo_y
    type(meshdiminfo), pointer :: diminfo_z
    !-------------------------------------------------
    
    i_size = 0
    do i=1, size(mesh3d%rcdomIJKP2LCMeshID,1)
      n = mesh3d%rcdomIJKP2LCMeshID(i,1,1,1)
      lcmesh => mesh3d%lcmesh_list(n)
      i_size =i_size + lcmesh%NeX * lcmesh%refElem3D%Nnode_h1D
    end do
 
    j_size = 0
    do j=1, size(mesh3d%rcdomIJKP2LCMeshID,2)
      n = mesh3d%rcdomIJKP2LCMeshID(1,j,1,1)
      lcmesh => mesh3d%lcmesh_list(n)
      j_size = j_size + lcmesh%NeY * lcmesh%refElem3D%Nnode_h1D
    end do
 
    k_size = 0
    do k=1, size(mesh3d%rcdomIJKP2LCMeshID,3)
      n = mesh3d%rcdomIJKP2LCMeshID(1,1,k,1)
      lcmesh => mesh3d%lcmesh_list(n)
      k_size = k_size + lcmesh%NeZ * lcmesh%refElem3D%Nnode_v
    end do
 
    k_size = k_size * size(mesh3d%rcdomIJKP2LCMeshID,4)
  
    diminfo_x => mesh3d%dimInfo(meshbase3d_dimtypeid_x)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_x),    &
      diminfo_x, "X", 1, (/ diminfo_x%name /), (/ i_size /), &
      dim_name_postfix )
 
    diminfo_y => mesh3d%dimInfo(meshbase3d_dimtypeid_y)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_y),    &
      diminfo_y, "Y", 1, (/ diminfo_y%name /), (/ j_size /), &
      dim_name_postfix )
 
    diminfo_z => mesh3d%dimInfo(meshbase3d_dimtypeid_z)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_z),           &
      diminfo_z, "Z", 1, (/ diminfo_z%name /), (/ k_size /),        &
      dim_name_postfix, positive_down=(/ diminfo_z%positive_down /) )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_zt)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_zt),  &
      diminfo, "ZT", 1, (/ diminfo_z%name /), (/ k_size /), &
      dim_name_postfix )
 
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xy)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xy),     &
      diminfo, "XY", 2, (/ diminfo_x%name, diminfo_y%name  /), &
      (/ i_size, j_size /), &
      dim_name_postfix )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xyt)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xyt),     &
      diminfo, "XY", 2, (/ diminfo_x%name, diminfo_y%name  /),  &
      (/ i_size, j_size /), &
      dim_name_postfix )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xyz)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xyz),                    &
      diminfo, "XYZ", 3, (/ diminfo_x%name, diminfo_y%name, diminfo_z%name /), &
      (/ i_size, j_size, k_size /), &
      dim_name_postfix )
  
    diminfo => mesh3d%dimInfo(meshbase3d_dimtypeid_xyzt)
    call set_dimension( dimsinfo(meshbase3d_dimtypeid_xyzt),                    &
      diminfo, "XYZT", 3, (/ diminfo_x%name, diminfo_y%name, diminfo_z%name /), &
      (/ i_size, j_size, k_size /),                                             &
      dim_name_postfix )
 
    return
  end subroutine file_common_meshfield_get_dims3d_cubedsphere
 
!OCL SERIAL
  subroutine file_common_meshfield_get_axis3d( mesh3D, dimsinfo, x, y, z, &
    force_uniform_grid )
    implicit none
 
    class(meshcubedom3d), target, intent(in) :: mesh3D  
    type(file_common_meshfield_diminfo), intent(in) :: dimsinfo(MESHBASE3D_DIMTYPE_NUM)
    real(RP), intent(out) :: x(dimsinfo(MeshBase3D_DIMTYPEID_X)%size)
    real(RP), intent(out) :: y(dimsinfo(MeshBase3D_DIMTYPEID_Y)%size)
    real(RP), intent(out) :: z(dimsinfo(MeshBase3D_DIMTYPEID_Z)%size)
    logical, intent(in), optional :: force_uniform_grid
 
    integer :: n, kelem
    integer :: i, j, k
    type(elementbase3d), pointer :: refElem
    type(localmesh3d), pointer :: lcmesh
 
    integer :: is, js, ks, ie, je, ke, igs, jgs, kgs
    integer :: Nnode_h1D, Nnode_v
 
    logical :: uniform_grid = .false.
    real(RP), allocatable :: x_local(:)
    real(RP), allocatable :: y_local(:)
    real(RP), allocatable :: z_local(:)
    !------------------------------------------------------------------------------------------
    
    if ( present(force_uniform_grid) ) uniform_grid = force_uniform_grid
 
    igs = 0; jgs = 0; kgs = 0
 
    do n=1 ,mesh3d%LOCAL_MESH_NUM
      lcmesh => mesh3d%lcmesh_list(n)
      refelem => lcmesh%refElem3D
      nnode_h1d = refelem%Nnode_h1D
      nnode_v   = refelem%Nnode_v
 
      allocate( x_local(nnode_h1d), y_local(nnode_h1d) )
      allocate( z_local(nnode_v) )
 
      do k=1, lcmesh%NeZ
      do j=1, lcmesh%NeY
      do i=1, lcmesh%NeX
        kelem = i + (j-1)*lcmesh%NeX + (k-1)*lcmesh%NeX*lcmesh%NeY
        if ( j==1 .and. k==1) then
          x_local(:) = lcmesh%pos_en(refelem%Fmask_h(1:nnode_h1d,1),kelem,1)
          if ( uniform_grid ) call get_uniform_grid1d( x_local, nnode_h1d )
 
          is = igs + 1 + (i-1)*nnode_h1d
          ie = is + nnode_h1d - 1
          x(is:ie) = x_local(:)
        end if
        if ( i==1 .and. k==1) then
          y_local(:) = lcmesh%pos_en(refelem%Fmask_h(1:nnode_h1d,4),kelem,2)
          if ( uniform_grid ) call get_uniform_grid1d( y_local, nnode_h1d )
 
          js = jgs + 1 + (j-1)*nnode_h1d
          je = js + nnode_h1d - 1
          y(js:je) = y_local(:)
        end if
        if ( i==1 .and. j==1) then
          z_local(:) = lcmesh%pos_en(refelem%Colmask(:,1),kelem,3)
          if ( uniform_grid ) call get_uniform_grid1d( z_local, nnode_v )
 
          ks = kgs + 1 + (k-1)*nnode_v
          ke = ks + nnode_v - 1
          z(ks:ke) = z_local(:)
        end if
      end do
      end do
      end do
 
      igs = ie; jgs = je; kgs = ke
      deallocate( x_local, y_local )
      deallocate( z_local )
    end do
 
    return
  end subroutine file_common_meshfield_get_axis3d
 
!OCL SERIAL
  subroutine file_common_meshfield_get_axis3d_cubedsphere( mesh3D, dimsinfo, x, y, z )
    implicit none
 
    class(meshcubedspheredom3d), target, intent(in) :: mesh3D  
    type(file_common_meshfield_diminfo), intent(in) :: dimsinfo(MESHBASE3D_DIMTYPE_NUM)
    real(RP), intent(out) :: x(dimsinfo(MeshBase3D_DIMTYPEID_X)%size)
    real(RP), intent(out) :: y(dimsinfo(MeshBase3D_DIMTYPEID_Y)%size)
    real(RP), intent(out) :: z(dimsinfo(MeshBase3D_DIMTYPEID_Z)%size)
 
    integer :: n, ni, nj, nk, np
    integer :: kelem
    integer :: i, j, k
    type(elementbase3d), pointer :: refElem
    type(localmesh3d), pointer :: lcmesh
 
    integer :: is, js, ks, ie, je, ke, igs, jgs, kgs
 
    logical :: uniform_grid = .false.
    real(RP), allocatable :: x_local(:)
    real(RP), allocatable :: y_local(:)
    real(RP), allocatable :: z_local(:)
    !-------------------------------------------------
    
    igs = 0; jgs = 0; kgs = 0
    
    do np=1, size(mesh3d%rcdomIJKP2LCMeshID,4) 
    do nk=1, size(mesh3d%rcdomIJKP2LCMeshID,3) 
    do nj=1, size(mesh3d%rcdomIJKP2LCMeshID,2) 
    do ni=1, size(mesh3d%rcdomIJKP2LCMeshID,1)
      n = mesh3d%rcdomIJKP2LCMeshID(ni,nj,nk,np)
      lcmesh => mesh3d%lcmesh_list(n)
      refelem => lcmesh%refElem3D
 
      allocate( x_local(refelem%Nnode_h1D), y_local(refelem%Nnode_h1D), z_local(refelem%Nnode_v) )
 
      do k=1, lcmesh%NeZ
      do j=1, lcmesh%NeY
      do i=1, lcmesh%NeX
        kelem = i + (j-1) * lcmesh%NeX + (k-1) * lcmesh%NeX * lcmesh%NeY
        if ( j==1 .and. nj == 1 .and. k==1 .and. nk == 1 .and. np == 1) then
          x_local(:) = lcmesh%pos_en(refelem%Fmask_h(1:refelem%Nnode_h1D,1),kelem,1)
 
          is = igs + 1 + (i-1)*refelem%Nnode_h1D
          ie = is + refelem%Nnode_h1D - 1
          x(is:ie) = x_local(:)
        end if
        if ( i==1 .and. ni == 1 .and. k==1 .and. nk == 1 .and. np == 1 ) then
          y_local(:) = lcmesh%pos_en(refelem%Fmask_h(1:refelem%Nnode_h1D,4),kelem,2) 
 
          js = jgs + 1 + (j-1)*refelem%Nnode_h1D
          je = js + refelem%Nnode_h1D - 1
          y(js:je) = y_local(:)
        end if
        if ( i==1 .and. ni == 1 .and. j == 1 .and. nj == 1 ) then
          z_local(:) = lcmesh%pos_en(refelem%Colmask(:,1),kelem,3) &
                     + ( lcmesh%panelID - 1.0_rp ) * ( mesh3d%zmax_gl - mesh3d%zmin_gl )
 
          ks = kgs + 1 + (k-1)*refelem%Nnode_v
          ke = ks + refelem%Nnode_v - 1
          z(ks:ke) = z_local(:)
        end if
      end do
      end do
      end do
 
      igs = ie; jgs = je; kgs = ke
      deallocate( x_local, y_local, z_local )
    end do
    end do
    end do
    end do
 
    return
  end subroutine file_common_meshfield_get_axis3d_cubedsphere
 
!OCL_SERIAL
  subroutine file_common_meshfield_put_field3d_cartesbuf( mesh3D, field3D, &
    buf, force_uniform_grid )
    use scale_polynomial, only: &
      polynomial_genlegendrepoly_sub
    implicit none
    class(meshcubedom3d), target, intent(in) :: mesh3d
    class(meshfield3d), intent(in) :: field3d
    real(rp), intent(inout) :: buf(:,:,:)
    logical, intent(in), optional :: force_uniform_grid
 
    integer :: n, kelem1
    integer :: i0, j0, k0, i1, j1, k1, i2, j2, k2, i, j, k
    type(localmesh3d), pointer :: lcmesh
    type(elementbase3d), pointer :: refelem
    integer :: i0_s, j0_s, k0_s, indx
 
    logical :: uniform_grid = .false.
    integer :: nnode_h1d, nnode_v
    real(rp), allocatable :: x_local(:)
    real(rp) :: x_local0, delx
    real(rp), allocatable :: y_local(:)
    real(rp) :: y_local0, dely
    real(rp), allocatable :: z_local(:)
    real(rp) :: z_local0, delz
    real(rp) :: ox(1), oy(1), oz(1)
    real(rp), allocatable :: spectral_coef(:)
    real(rp), allocatable :: p1d_ori_x(:,:)
    real(rp), allocatable :: p1d_ori_y(:,:)
    real(rp), allocatable :: p1d_ori_z(:,:)
    integer :: l, p1, p2, p3
    !----------------------------------------------------
 
    if ( present(force_uniform_grid) ) uniform_grid = force_uniform_grid
 
    i0_s = 0; j0_s = 0; k0_s = 0
 
    do k0=1, size(mesh3d%rcdomIJK2LCMeshID,3)  
      do j0=1, size(mesh3d%rcdomIJK2LCMeshID,2)
        do i0=1, size(mesh3d%rcdomIJK2LCMeshID,1)
          n =  mesh3d%rcdomIJK2LCMeshID(i0,j0,k0)
 
          lcmesh => mesh3d%lcmesh_list(n)
          refelem => lcmesh%refElem3D
          nnode_h1d = refelem%Nnode_h1D
          nnode_v   = refelem%Nnode_v
 
          if ( uniform_grid ) then
            allocate( x_local(nnode_h1d), y_local(nnode_h1d) )
            allocate( z_local(nnode_v) ) 
            allocate( spectral_coef(refelem%Np) )
            allocate( p1d_ori_x(1,nnode_h1d), p1d_ori_y(1,nnode_h1d) )
            allocate( p1d_ori_z(1,nnode_v) )     
          end if
 
          !$omp parallel do collapse(2) private( kelem1, &
          !$omp i, i1, i2, j, j2, k, k2, indx,                               &
          !$omp x_local, x_local0, y_local, y_local0, z_local, z_local0,     &
          !$omp delx, dely, delz, ox, oy, oz,                                &
          !$omp spectral_coef, P1D_ori_x, P1D_ori_y, P1D_ori_z,              &
          !$omp p1, p2, p3, l                                                )        
          do k1=1, lcmesh%NeZ
          do j1=1, lcmesh%NeY
          do i1=1, lcmesh%NeX
            kelem1 = i1 + (j1-1)*lcmesh%NeX + (k1-1)*lcmesh%NeX*lcmesh%NeY
 
            if ( uniform_grid ) then
              x_local(:) = lcmesh%pos_en(refelem%Fmask_h(1:nnode_h1d,1),kelem1,1)
              x_local0 = x_local(1); delx = x_local(nnode_h1d) - x_local0
              y_local(:) = lcmesh%pos_en(refelem%Fmask_h(1:nnode_h1d,4),kelem1,2)
              y_local0 = y_local(1); dely = y_local(nnode_h1d) - y_local0
              z_local(:) = lcmesh%pos_en(refelem%Colmask(:,1),kelem1,3)
              z_local0 = z_local(1); delz = z_local(nnode_v  ) - z_local0
              call get_uniform_grid1d( x_local, nnode_h1d )
              call get_uniform_grid1d( y_local, nnode_h1d )
              call get_uniform_grid1d( z_local, nnode_v   )
 
              spectral_coef(:) = matmul(refelem%invV(:,:), field3d%local(n)%val(:,kelem1))
              do k2=1, nnode_v
              do j2=1, nnode_h1d
              do i2=1, nnode_h1d
                ox(1) = - 1.0_rp + 2.0_rp * (x_local(i2) - x_local0) / delx
                oy(1) = - 1.0_rp + 2.0_rp * (y_local(j2) - y_local0) / dely
                oz(1) = - 1.0_rp + 2.0_rp * (z_local(k2) - z_local0) / delz
      
                call polynomial_genlegendrepoly_sub( refelem%PolyOrder_h, ox, p1d_ori_x )
                call polynomial_genlegendrepoly_sub( refelem%PolyOrder_h, oy, p1d_ori_y )
                call polynomial_genlegendrepoly_sub( refelem%PolyOrder_v, oz, p1d_ori_z )
 
                i = i0_s + i2 + (i1-1)*nnode_h1d
                j = j0_s + j2 + (j1-1)*nnode_h1d
                k = k0_s + k2 + (k1-1)*nnode_v
                buf(i,j,k) = 0.0_rp 
                do p3=1, nnode_v
                do p2=1, nnode_h1d
                do p1=1, nnode_h1d
                  l = p1 + (p2-1)*nnode_h1d + (p3-1)*nnode_h1d**2
                  buf(i,j,k) = buf(i,j,k) + &
                      ( p1d_ori_x(1,p1) * p1d_ori_y(1,p2) * p1d_ori_z(1,p3) )                 &
                    * sqrt((dble(p1-1) + 0.5_rp)*(dble(p2-1) + 0.5_rp)*(dble(p3-1) + 0.5_rp)) &
                    * spectral_coef(l)
                end do
                end do
                end do            
              end do
              end do
              end do
            else
              do k2=1, nnode_v
              do j2=1, nnode_h1d
              do i2=1, nnode_h1d
                i = i0_s + i2 + (i1-1)*nnode_h1d
                j = j0_s + j2 + (j1-1)*nnode_h1d
                k = k0_s + k2 + (k1-1)*nnode_v
                indx = i2 + (j2-1)*nnode_h1d + (k2-1)*nnode_h1d**2
                buf(i,j,k) = field3d%local(n)%val(indx,kelem1)
              end do
              end do
              end do
            end if
          end do
          end do
          end do
 
          if ( uniform_grid ) then
            deallocate( x_local, y_local )
            deallocate( z_local )
            deallocate( spectral_coef )
          end if
 
          i0_s = i0_s + lcmesh%NeX * refelem%Nnode_h1D
        end do
        j0_s = j0_s + lcmesh%NeY * refelem%Nnode_h1D
        i0_s = 0
      end do
      k0_s = k0_s + lcmesh%NeZ * refelem%Nnode_v
      j0_s = 0
    end do
 
    return
  end subroutine file_common_meshfield_put_field3d_cartesbuf
 
!OCL SERIAL
  subroutine file_common_meshfield_put_field3d_cubedsphere_cartesbuf( mesh3D, field3D, &
    buf )
    implicit none
    class(meshcubedspheredom3d), target, intent(in) :: mesh3d
    class(meshfield3d), intent(in) :: field3d
    real(rp), intent(inout) :: buf(:,:,:)
 
    integer :: kelem1
    integer :: n, i0, j0, k0, p0
    integer :: i1, j1, k1, i2, j2, k2, i, j, k
    type(localmesh3d), pointer :: lcmesh
    type(elementbase3d), pointer :: refelem
    integer :: i0_s, j0_s, k0_s
    integer :: nnode_h1d
    integer :: nnode_v
    !------------------------------------------------
 
    i0_s = 0; j0_s = 0; k0_s = 0
 
    do p0=1, size(mesh3d%rcdomIJKP2LCMeshID,4)    
      do k0=1, size(mesh3d%rcdomIJKP2LCMeshID,3)
        do j0=1, size(mesh3d%rcdomIJKP2LCMeshID,2)
          do i0=1, size(mesh3d%rcdomIJKP2LCMeshID,1)
            n = mesh3d%rcdomIJKP2LCMeshID(i0,j0,k0,p0)
 
            lcmesh => mesh3d%lcmesh_list(n)
            refelem => lcmesh%refElem3D
            nnode_h1d = refelem%Nnode_h1D
            nnode_v = refelem%Nnode_v
                
            do k1=1, lcmesh%NeZ
            do j1=1, lcmesh%NeY
            do i1=1, lcmesh%NeX
              kelem1 = i1 + (j1-1)*lcmesh%NeX + (k1-1)*lcmesh%NeX*lcmesh%NeY
 
              do k2=1, nnode_v
              do j2=1, nnode_h1d
              do i2=1, nnode_h1d
                i = i0_s + i2 + (i1-1)*nnode_h1d
                j = j0_s + j2 + (j1-1)*nnode_h1d
                k = k0_s + k2 + (k1-1)*nnode_v
                buf(i,j,k) = field3d%local(n)%val(i2+(j2-1)*nnode_h1d+(k2-1)*nnode_h1d**2,kelem1)
              end do
              end do
              end do
            end do
            end do
            end do 
            
            i0_s = i0_s + lcmesh%NeX * refelem%Nnode_h1D
          end do
          j0_s = j0_s + lcmesh%NeY * refelem%Nnode_h1D
          i0_s = 0
        end do
        k0_s = k0_s + lcmesh%NeZ * refelem%Nnode_v
        j0_s = 0
      end do
    end do
    
    return
  end subroutine file_common_meshfield_put_field3d_cubedsphere_cartesbuf
 
  subroutine file_common_meshfield_set_cartesbuf_field3d( mesh3D, buf, &
    field3D )
    implicit none
    class(meshcubedom3d), target, intent(in) :: mesh3d
    real(rp), intent(in) :: buf(:,:,:)
    class(meshfield3d), intent(inout) :: field3d
 
    integer :: n
    integer :: i0, j0, k0
    type(localmesh3d), pointer :: lcmesh
    type(elementbase3d), pointer :: refelem
    integer :: i0_s, j0_s, k0_s
    !----------------------------------------------------
 
    i0_s = 0; j0_s = 0; k0_s = 0
 
    do k0=1, size(mesh3d%rcdomIJK2LCMeshID,3)  
      do j0=1, size(mesh3d%rcdomIJK2LCMeshID,2)
        do i0=1, size(mesh3d%rcdomIJK2LCMeshID,1)
          n = mesh3d%rcdomIJK2LCMeshID(i0,j0,k0)
          lcmesh => mesh3d%lcmesh_list(n)
          refelem => lcmesh%refElem3D
 
          call file_common_meshfield_set_cartesbuf_field3d_local(  &
            lcmesh, buf(:,:,:), i0_s, j0_s, k0_s,                  &
            field3d%local(n)%val(:,:)                              )
          !$acc update device(field3d%local(n)%val)
 
          i0_s = i0_s + lcmesh%NeX * refelem%Nnode_h1D
        end do
        j0_s = j0_s + lcmesh%NeY * refelem%Nnode_h1D
        i0_s = 0
      end do
      k0_s = k0_s + lcmesh%NeZ * refelem%Nnode_v
      j0_s = 0
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field3d
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field3d_cubedsphere( mesh3D, buf, &
    field3D )
    implicit none
    class(meshcubedspheredom3d), target, intent(in) :: mesh3d
    real(rp), intent(in) :: buf(:,:,:)
    class(meshfield3d), intent(inout) :: field3d
 
    integer :: n
    integer :: i0, j0, k0, p0
    type(localmesh3d), pointer :: lcmesh
    type(elementbase3d), pointer :: refelem
    integer :: i0_s, j0_s, k0_s, p0_s
    !----------------------------------------------------
 
    i0_s = 0; j0_s = 0; k0_s = 0; p0_s = 0
 
    do p0=1, size(mesh3d%rcdomIJKP2LCMeshID,4)
      do k0=1, size(mesh3d%rcdomIJKP2LCMeshID,3)
        do j0=1, size(mesh3d%rcdomIJKP2LCMeshID,2)
          do i0=1, size(mesh3d%rcdomIJKP2LCMeshID,1)
            n = mesh3d%rcdomIJKP2LCMeshID(i0,j0,k0,p0)
            lcmesh => mesh3d%lcmesh_list(n)
            refelem => lcmesh%refElem3D
 
            call file_common_meshfield_set_cartesbuf_field3d_local(  &
              lcmesh, buf(:,:,:), i0_s, j0_s, k0_s,                  &
              field3d%local(n)%val(:,:)                              )
            !$acc update device(field3d%local(n)%val)
 
            i0_s = i0_s + lcmesh%NeX * refelem%Nnode_h1D
          end do
          j0_s = j0_s + lcmesh%NeY * refelem%Nnode_h1D
          i0_s = 0
        end do
        k0_s = k0_s + lcmesh%NeZ * refelem%Nnode_v
        j0_s = 0
      end do
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field3d_cubedsphere
 
!OCL SERIAL
  subroutine file_common_meshfield_set_cartesbuf_field3d_local( &
    lcmesh, buf, i0_s, j0_s, k0_s,                              &
    val )
    implicit none
    type(localmesh3d), intent(in) :: lcmesh
    real(rp), intent(in) :: buf(:,:,:)
    integer, intent(in) :: i0_s, j0_s, k0_s
    real(rp), intent(inout) :: val(lcmesh%refelem3d%np,lcmesh%nea)
 
    integer :: kelem1
    integer :: i1, j1, k1, i2, j2, k2, i, j, k
    type(elementbase3d), pointer :: refelem
    integer :: indx
    !----------------------------------------------------
 
    refelem => lcmesh%refElem3D
 
    !$omp parallel do collapse(3) private( &
    !$omp kelem1, k2,j2,i2, i,j,k, indx    )
    do k1=1, lcmesh%NeZ
    do j1=1, lcmesh%NeY
    do i1=1, lcmesh%NeX
      kelem1 = i1 + (j1-1)*lcmesh%NeX + (k1-1)*lcmesh%NeX*lcmesh%NeY
      do k2=1, refelem%Nnode_v
      do j2=1, refelem%Nnode_h1D
      do i2=1, refelem%Nnode_h1D
        i = i0_s + i2 + (i1-1)*refelem%Nnode_h1D
        j = j0_s + j2 + (j1-1)*refelem%Nnode_h1D
        k = k0_s + k2 + (k1-1)*refelem%Nnode_v
        indx = i2 + (j2-1)*refelem%Nnode_h1D + (k2-1)*refelem%Nnode_h1D**2
        val(indx,kelem1) = buf(i,j,k)
      end do
      end do
      end do
    end do
    end do
    end do
 
    return
  end subroutine file_common_meshfield_set_cartesbuf_field3d_local
 
  function file_common_meshfield_get_dtype( datatype ) result( dtype )
    use scale_file_h, only: &
      file_real8, file_real4
    use scale_prc, &
      only: prc_abort
 
    implicit none
 
    character(*), intent(in) :: datatype
    integer :: dtype
    !--------------------------
 
    ! dtype is used to define the data type of axis variables in file
    if    ( datatype == 'REAL8' ) then
      dtype = file_real8
    elseif( datatype == 'REAL4' ) then
        dtype = file_real4
    else
      if    ( rp == 8 ) then
        dtype = file_real8
      elseif( rp == 4 ) then
        dtype = file_real4
      else
        log_error("file_restart_meshfield_get_dtype",*) 'unsupported data type. Check!', trim(datatype)
        call prc_abort
      endif
    endif
 
    return
  end function file_common_meshfield_get_dtype
 
  !- private -----------------------------------------------------------------------
 
  subroutine get_uniform_grid1d( pos1D, Np )
    implicit none
    integer, intent(in) :: np
    real(rp), intent(inout) :: pos1d(np)
 
    real(rp) :: del
    integer :: i
    !-----------------------------------------------
 
    del = ( pos1d(np) - pos1d(1) ) / dble(np)
    pos1d(1) = pos1d(1) + 0.5_rp * del
    do i=2, np
      pos1d(i) = pos1d(i-1) + del
    end do
 
    return
  end subroutine get_uniform_grid1d
 
  !> Set dimension information for file output
  subroutine set_dimension( dim,           & ! (out)
    diminfo, dim_type, ndims, dims, count, & ! (in)
    dim_postfix, positive_down             )  ! (in)
    implicit none
 
    type(file_common_meshfield_diminfo), intent(out) :: dim
    type(meshdiminfo), intent(in) :: diminfo
    character(*), intent(in) :: dim_type
    integer, intent(in) :: ndims
    character(len=*), intent(in) :: dims(ndims)
    integer, intent(in) :: count(ndims)
    character(*), intent(in) :: dim_postfix
    logical, intent(in), optional :: positive_down(ndims)
 
    integer :: d
    !----------------------------------------------------
 
    dim%name = trim(diminfo%name) // trim(dim_postfix)
    dim%unit = diminfo%unit
    dim%desc = diminfo%desc
    dim%type = trim(dim_type) // trim(dim_postfix)
    dim%ndim = ndims
    dim%size = 1
    do d=1, ndims
      dim%dims(d) = trim(dims(d)) // trim(dim_postfix)
      dim%count(d) = count(d)
      dim%size = dim%size * count(d)
      if ( present(positive_down) ) then
        dim%positive_down(d) = positive_down(d)
      else
        dim%positive_down(d) = .false.        
      end if
    end do
 
    return
  end subroutine set_dimension
 
end module scale_file_common_meshfield