scale_mesh_cubedspheredom3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_mesh_cubedspheredom3d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_io
  use scale_precision
 
  use scale_mesh_base3d, only: &
    meshbase3d, meshbase3d_init, meshbase3d_final,                               &
    meshbase3d_setgeometricinfo,                                                 &
    meshbase3d_dimtype_num,                                                      &
    meshbase3d_dimtypeid_x, meshbase3d_dimtypeid_y, meshbase3d_dimtypeid_z,      &
    meshbase3d_dimtypeid_xyz, meshbase3d_dimtypeid_zt, meshbase3d_dimtypeid_xyzt
  use scale_localmesh_3d, only: localmesh3d
  use scale_element_hexahedral, only: hexahedralelement
 
  use scale_mesh_base2d, only: &
    meshbase2d, meshbase2d_init, meshbase2d_final, &
    meshbase2d_setgeometricinfo
  use scale_mesh_cubedspheredom2d, only: &
    meshcubedspheredom2d, meshcubedspheredom2d_setuplocaldom
  use scale_localmesh_2d, only: localmesh2d
  use scale_element_quadrilateral, only: quadrilateralelement
 
  use scale_element_base, only: elementbase2d, elementbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  type, extends(meshbase3d), public :: meshcubedspheredom3d
    integer :: negx
    integer :: negy
    integer :: negz
    
    real(rp), public :: xmin_gl, xmax_gl
    real(rp), public :: ymin_gl, ymax_gl
    real(rp), public :: zmin_gl, zmax_gl
 
    real(rp), allocatable :: fz(:)
 
    integer, allocatable :: rcdomijkp2lcmeshid(:,:,:,:)
 
    real(rp) :: rplanet
 
    type(meshcubedspheredom2d) :: mesh2d
    type(quadrilateralelement) :: refelem2d
 
    logical :: shallow_approx
  contains
    procedure :: init => meshcubedspheredom3d_init
    procedure :: final => meshcubedspheredom3d_final
    procedure :: generate => meshcubedspheredom3d_generate
    procedure :: assigndomid => messhcubedspheredom3d_assigndomid
    procedure :: getmesh2d => meshcubedspheredom3d_getmesh2d
    procedure :: set_geometric_with_vcoord => meshcubedom3d_set_geometric_with_vcoord
  end type meshcubedspheredom3d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
 
  private :: meshcubedspheredom3d_calc_normal
  private :: meshcubedspheredom3d_coord_conv
  private :: meshcubedspheredom3d_set_metric
  private :: fill_halo_metric
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
!OCL SERIAL
  subroutine meshcubedspheredom3d_init( this, &
    NeGX, NeGY, NeGZ, RPlanet,                &
    dom_zmin, dom_zmax,                       &
    refElem, NLocalMeshPerPrc,                &
    nproc, myrank,                            &
    FZ, shallow_approx                        )
    
    use scale_const, only: &
      pi => const_pi
    implicit none
 
    class(meshcubedspheredom3d), intent(inout) :: this
    integer, intent(in) :: NeGX
    integer, intent(in) :: NeGY
    integer, intent(in) :: NeGZ
    real(RP), intent(in) :: RPlanet
    real(RP), intent(in) :: dom_Zmin
    real(RP), intent(in) :: dom_zmax
    type(hexahedralelement), intent(in), target :: refElem
    integer, intent(in) :: NLocalMeshPerPrc
    integer, intent(in), optional :: nproc
    integer, intent(in), optional :: myrank
    real(RP), intent(in), optional :: FZ(NeGZ+1)
    logical, intent(in), optional :: shallow_approx
 
    integer :: k
    real(RP) :: dz
    !-----------------------------------------------------------------------------
 
    this%NeGX = negx
    this%NeGY = negy
    this%NeGZ = negz
 
    this%xmin_gl = - 0.25_rp * pi 
    this%xmax_gl = + 0.25_rp * pi 
    this%ymin_gl = - 0.25_rp * pi 
    this%ymax_gl = + 0.25_rp * pi 
    this%zmin_gl = dom_zmin
    this%zmax_gl = dom_zmax
    this%RPlanet = rplanet
    this%dom_vol = 4.0_rp / 3.0_rp * pi * ( ( dom_zmax + rplanet )**3 - ( dom_zmin + rplanet )**3 )
 
 
    !- Fz
    allocate( this%FZ(this%NeGZ+1) )
    if ( present(fz) ) then
      this%FZ(:) = fz(:)
    else
      this%FZ(1          ) = dom_zmin
      this%FZ(this%NeGZ+1) = dom_zmax
      dz = (dom_zmax - dom_zmin) / dble(this%NeGZ)
      do k=2, this%NeGZ
        this%FZ(k) = this%FZ(k-1) + dz
      end do
    end if
 
    !-
    if ( present(shallow_approx) ) then
      this%shallow_approx = shallow_approx
    else
      this%shallow_approx = .true.
    end if
 
    !--
    call meshbase3d_init( this, refelem, nlocalmeshperprc, 6, &
      nproc, myrank )
 
    !---
    call this%refElem2D%Init( this%refElem3D%PolyOrder_h, refelem%IsLumpedMatrix() )
    call this%mesh2D%Init( negx, negy, rplanet, this%refElem2D, nlocalmeshperprc, &
      nproc, myrank )
      
    !-- Modify the information of dimension for the cubed sphere mesh
    call this%SetDimInfo( meshbase3d_dimtypeid_x, "x", "1", "X-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_y, "y", "1", "Y-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_z, "z", "m", "Z-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_xyz, "xyz", "1", "XYZ-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_zt, "zt", "1", "XYZ-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_xyzt, "xyzt", "1", "XYZ-coordinate" )
  
    return
  end subroutine meshcubedspheredom3d_init
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_final( this )
    implicit none
    class(meshcubedspheredom3d), intent(inout) :: this
    !-----------------------------------------------------------------------------
    
    if (this%isGenerated) then
      if ( allocated(this%rcdomIJKP2LCMeshID) ) then
        deallocate( this%rcdomIJKP2LCMeshID )
      end if
    else
      if ( allocated( this%FZ ) ) deallocate( this%FZ )
    end if
 
    call this%mesh2D%Final()
    call this%refElem2D%Final()
 
    call meshbase3d_final( this )
 
    return
  end subroutine meshcubedspheredom3d_final
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_getmesh2d( this, ptr_mesh2D )
    implicit none
    class(meshcubedspheredom3d), intent(in), target :: this
    class(meshbase2d), pointer, intent(out) :: ptr_mesh2D
    !-------------------------------------------------------
 
    ptr_mesh2d => this%mesh2D
    return
  end subroutine meshcubedspheredom3d_getmesh2d
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_generate( this )
    use scale_mesh_cubedspheredom2d, only: &
      meshcubedspheredom2d_check_division_params
    implicit none
 
    class(meshcubedspheredom3d), intent(inout), target :: this
 
    integer :: n
    type(localmesh3d), pointer :: mesh
 
    integer :: tileID_table(this%LOCAL_MESH_NUM, this%PRC_NUM)
    integer :: panelID_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pi_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pj_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pk_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
 
    integer :: NprcX_lc, NprcY_lc, NprcZ_lc
    integer :: tileID   
    !-----------------------------------------------------------------------------
 
    call meshcubedspheredom2d_check_division_params( &
      nprcx_lc, nprcy_lc,                            &
      this%PRC_NUM, this%LOCAL_MESH_NUM_global,      &
      .true. )
 
    nprcz_lc = 1
 
    !--- Construct the connectivity of patches  (only master node)
 
    call this%AssignDomID( & 
      nprcx_lc, nprcy_lc, nprcz_lc,  & ! (in)
      tileid_table, panelid_table,   & ! (out)
      pi_table, pj_table, pk_table   ) ! (out)
 
    do n=1, this%LOCAL_MESH_NUM
      mesh => this%lcmesh_list(n)
      tileid = tileid_table(n, mesh%PRC_myrank+1)
 
      call meshcubedspheredom3d_setuplocaldom( mesh, &
        tileid,  panelid_table(tileid),                             &
        pi_table(tileid), pj_table(tileid), pk_table(tileid),       &
        nprcx_lc, nprcy_lc, nprcz_lc,                               &
        this%xmin_gl, this%xmax_gl, this%ymin_gl, this%ymax_gl,     &
        this%zmin_gl, this%zmax_gl, this%RPlanet,                   &
        this%NeGX/nprcx_lc, this%NeGY/nprcy_lc, this%NeGZ/nprcz_lc, &
        this%FZ(:) )
 
      call meshcubedspheredom2d_setuplocaldom( this%mesh2D%lcmesh_list(n),   &
        tileid,  panelid_table(tileid),                                      &
        pi_table(tileid), pj_table(tileid), nprcx_lc, nprcy_lc,              &
        this%xmin_gl, this%xmax_gl, this%ymin_gl, this%ymax_gl,              &
        this%RPlanet, this%NeGX/nprcx_lc, this%NeGY/nprcy_lc                 )
 
      call mesh%SetLocalMesh2D( this%mesh2D%lcmesh_list(n) )
 
      !---
      ! write(*,*) "** my_rank=", mesh%PRC_myrank
      ! write(*,*) " tileID:", mesh%tileID
      ! write(*,*) " pnlID:", mesh%panelID, "-- i,j (within a panel)=", pi_table(tileID), pj_table(tileID)
      ! write(*,*) " local mesh:", n, "( total", this%LOCAL_MESH_NUM, ")"
      ! write(*,*) " panel_connect:", this%tilePanelID_globalMap(:,mesh%tileID)
      ! write(*,*) " tile_connect:", this%tileID_globalMap(:,mesh%tileID)
      ! write(*,*) " face_connect:", this%tileFaceID_globalMap(:,mesh%tileID)
      ! write(*,*) " domain size"
      ! write(*,*) "   NeX, NeY:", mesh%NeX, mesh%NeY
      ! write(*,*) "   [X], [Y]:",  mesh%xmin, mesh%xmax, ":", mesh%ymin, mesh%ymax
    end do
 
    ! Set lon&lat position and metrics with the cubed sphere mesh
    call meshcubedspheredom3d_set_metric( this )
 
    ! To set rcdomIJP2LCMeshID, call AssignDomID for 2D mesh
    call this%mesh2D%AssignDomID( & 
      nprcx_lc, nprcy_lc,            & ! (in)
      tileid_table, panelid_table,   & ! (out)
      pi_table, pj_table             ) ! (out)
 
    this%isGenerated = .true.
    this%mesh2D%isGenerated = .true.
    
    deallocate( this%FZ )
    
    return
  end subroutine meshcubedspheredom3d_generate
 
!OCL SERIAL
  subroutine meshcubedom3d_set_geometric_with_vcoord(this, lcdomID, GsqrtV_lc, zlev_lc, G13_lc, G23_lc)
    implicit none
    class(meshcubedspheredom3d), intent(inout), target :: this
    integer, intent(in) :: lcdomID
    real(RP), intent(in) :: GsqrtV_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    real(RP), intent(in) :: zlev_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    real(RP), intent(in) :: G13_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    real(RP), intent(in) :: G23_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
 
    integer :: ke
    class(localmesh3d), pointer :: lcmesh
    !-------------------------------------------------------
 
    lcmesh => this%lcmesh_list(lcdomid)
 
    !$omp parallel private(ke)
    !$omp do
    do ke=lcmesh%NeS, lcmesh%NeE
      lcmesh%zlev(:,ke) = zlev_lc(:,ke)
    end do
    !$omp do
    do ke=lcmesh%NeS, lcmesh%NeA   
      if ( this%shallow_approx ) then
        lcmesh%gam(:,ke) = 1.0_rp
      else
        lcmesh%gam(:,ke) = 1.0_rp + zlev_lc(:,ke) / this%RPlanet
      end if
 
      lcmesh%Gsqrt(:,ke) = gsqrtv_lc(:,ke) * lcmesh%gam(:,ke)**2 * lcmesh%Gsqrt(:,ke)
      lcmesh%GI3(:,ke,1) = g13_lc(:,ke)
      lcmesh%GI3(:,ke,2) = g23_lc(:,ke)
    end do
    !$omp end parallel
 
    return
  end subroutine meshcubedom3d_set_geometric_with_vcoord
 
  !- private ------------------------------
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_setuplocaldom( lcmesh,        &
    tileID, panelID,                                            &
    i, j, k, NprcX, NprcY, NprcZ,                               &
    dom_xmin, dom_xmax, dom_ymin, dom_ymax, dom_zmin, dom_zmax, &
    planet_radius, NeX, NeY, NeZ,                               &
    FZ )
 
    use scale_meshutil_cubedsphere3d, only: &
      meshutilcubedsphere3d_genconnectivity,   &
      meshutilcubedsphere3d_gencubedomain,     &
      meshutilcubedsphere3d_buildinteriormap,  &
      meshutilcubedsphere3d_genpatchboundarymap
    
    use scale_localmesh_base, only: bctype_interior
 
    implicit none
      
    type(localmesh3d), intent(inout) :: lcmesh
    integer, intent(in) :: tileID
    integer, intent(in) :: panelID
    integer, intent(in) :: i, j, k
    integer, intent(in) :: NprcX, NprcY, NprcZ
    real(RP), intent(in) :: dom_xmin, dom_xmax
    real(RP), intent(in) :: dom_ymin, dom_ymax
    real(RP), intent(in) :: dom_zmin, dom_zmax
    real(RP), intent(in) :: planet_radius
    integer, intent(in) :: NeX, NeY, NeZ
    real(RP), intent(in) :: FZ(NeZ*NprcZ+1)
 
    class(elementbase3d), pointer :: elem
    real(RP) :: delx, dely
    real(RP) :: FZ_lc(NeZ+1)
 
    integer :: ii, jj, kk
    integer :: ke
    !-----------------------------------------------------------------------------
 
    elem => lcmesh%refElem3D
    lcmesh%tileID = tileid
    lcmesh%panelID = panelid
        
    !--
    lcmesh%Ne   = nex * ney * nez
    lcmesh%Nv  = (nex + 1)*(ney + 1)*(nez + 1)
    lcmesh%NeS = 1
    lcmesh%NeE = lcmesh%Ne
    lcmesh%NeA = lcmesh%Ne + 2*(nex + ney)*nez + 2*nex*ney
 
    lcmesh%NeX = nex
    lcmesh%NeY = ney
    lcmesh%NeZ = nez
 
    lcmesh%Ne2D  = nex * ney
    lcmesh%Ne2DA = nex * ney + 2*(nex + ney)
 
    !--
    delx = ( dom_xmax - dom_xmin ) / dble(nprcx)
    dely = ( dom_ymax - dom_ymin ) / dble(nprcy)
    fz_lc(:) = fz((k-1)*nez+1:k*nez+1)
    lcmesh%xmin = dom_xmin + (i-1)*delx
    lcmesh%xmax = dom_xmin +  i   *delx
    lcmesh%ymin = dom_ymin + (j-1)*dely
    lcmesh%ymax = dom_ymin +  j   *dely
    lcmesh%zmin = fz_lc(1)
    lcmesh%zmax = fz_lc(nez+1)
 
    !--
    allocate( lcmesh%pos_ev(lcmesh%Nv,3) )
    allocate( lcmesh%EToV(lcmesh%Ne,elem%Nv) )
    allocate( lcmesh%EToE(lcmesh%Ne,elem%Nfaces) )
    allocate( lcmesh%EToF(lcmesh%Ne,elem%Nfaces) )
    allocate( lcmesh%BCType(lcmesh%refElem%Nfaces,lcmesh%Ne) )
    allocate( lcmesh%VMapM(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%VMapP(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%MapM(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%MapP(elem%NfpTot, lcmesh%Ne) )
 
    allocate( lcmesh%EMap3Dto2D(lcmesh%Ne) )
 
    lcmesh%BCType(:,:) = bctype_interior
 
    !----
    call meshutilcubedsphere3d_gencubedomain( lcmesh%pos_ev, lcmesh%EToV, & ! (out)
      lcmesh%NeX, lcmesh%xmin, lcmesh%xmax,                               & ! (in)
      lcmesh%NeY, lcmesh%ymin, lcmesh%ymax,                               & ! (in) 
      lcmesh%NeZ, lcmesh%zmin, lcmesh%zmax, fz=fz_lc                      ) ! (in) 
    
    !---
    call meshbase3d_setgeometricinfo(lcmesh, meshcubedspheredom3d_coord_conv, meshcubedspheredom3d_calc_normal )
    
    !---
    call meshutilcubedsphere3d_genconnectivity( lcmesh%EToE, lcmesh%EToF, & ! (out)
      lcmesh%EToV, lcmesh%Ne, elem%Nfaces )                                 ! (in)
    
    !---
    call meshutilcubedsphere3d_buildinteriormap( lcmesh%VmapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP,  & ! (out)
      lcmesh%pos_en, lcmesh%pos_ev, lcmesh%EToE, lcmesh%EtoF, lcmesh%EtoV,                              & ! (in)
      elem%Fmask_h, elem%Fmask_v, lcmesh%Ne, lcmesh%Nv, elem%Np, elem%Nfp_h, elem%Nfp_v, elem%NfpTot,   & ! (in)
      elem%Nfaces_h, elem%Nfaces_v, elem%Nfaces ) 
 
    
    call meshutilcubedsphere3d_genpatchboundarymap( lcmesh%VMapB, lcmesh%MapB, lcmesh%VMapP,            & !(out)
      lcmesh%pos_en, lcmesh%xmin, lcmesh%xmax, lcmesh%ymin, lcmesh%ymax, lcmesh%zmin, lcmesh%zmax,      & ! (in)
      elem%Fmask_h, elem%Fmask_v, lcmesh%Ne, lcmesh%Nv, elem%Np, elem%Nfp_h, elem%Nfp_v, elem%NfpTot,   & ! (in)
      elem%Nfaces_h, elem%Nfaces_v, elem%Nfaces )  
    
    !---
    !$omp parallel do collapse(2) private(ii,ke)
    do kk=1, lcmesh%NeZ
    do jj=1, lcmesh%NeY
    do ii=1, lcmesh%NeX
      ke = ii + (jj-1) * lcmesh%NeX + (kk-1) * lcmesh%NeX * lcmesh%NeY
      lcmesh%EMap3Dto2D(ke) = ii + (jj-1) * lcmesh%NeX
    end do
    end do
    end do
    
    return
  end subroutine meshcubedspheredom3d_setuplocaldom
 
!OCL SERIAL
  subroutine messhcubedspheredom3d_assigndomid( this, &
    NprcX_lc, NprcY_lc, NprcZ_lc,                     &
    tileID_table, panelID_table,                      &
    pi_table, pj_table, pk_table )
  
    use scale_meshutil_cubedsphere3d, only: &       
      meshutilcubedsphere3d_buildglobalmap
    
    implicit none
 
    class(meshcubedspheredom3d), target, intent(inout) :: this    
    integer, intent(in) :: NprcX_lc
    integer, intent(in) :: NprcY_lc
    integer, intent(in) :: NprcZ_lc
    integer, intent(out) :: tileID_table(this%LOCAL_MESH_NUM, this%PRC_NUM)
    integer, intent(out) :: panelID_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer, intent(out) :: pi_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer, intent(out) :: pj_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer, intent(out) :: pk_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    
    integer :: n
    integer :: prc
    integer :: tileID
    integer :: is_lc, js_lc, ks_lc, ps_lc
    integer :: ilc_count, jlc_count, klc_count, plc_count
    integer :: ilc, jlc, klc, plc
    
    type(localmesh3d), pointer :: lcmesh
    !-----------------------------------------------------------------------------
 
    call meshutilcubedsphere3d_buildglobalmap( &
      panelid_table, pi_table, pj_table, pk_table,                                  & ! (out)
      this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap, & ! (out)
      this%LOCAL_MESH_NUM_global, nprcz_lc )                                          ! (in)
 
    !----
    
    do prc=1, this%PRC_NUM
    do n=1, this%LOCAL_MESH_NUM
      tileid = n + (prc-1)*this%LOCAL_MESH_NUM
      lcmesh => this%lcmesh_list(n)
      
      !-
      tileid_table(n,prc)                 = tileid
      this%tileID_global2localMap(tileid) = n
      this%PRCRank_globalMap(tileid)      = prc - 1
 
      !-
      if ( this%PRCRank_globalMap(tileid) == lcmesh%PRC_myrank ) then
        if (n==1) then
          is_lc = pi_table(tileid); ilc_count = 1
          js_lc = pj_table(tileid); jlc_count = 1
          ks_lc = pk_table(tileid); klc_count = 1
          ps_lc = panelid_table(tileid); plc_count = 1
        end if
        if(is_lc < pi_table(tileid)) ilc_count = ilc_count + 1
        if(js_lc < pj_table(tileid)) jlc_count = jlc_count + 1
        if(ks_lc < pk_table(tileid)) klc_count = klc_count + 1
        if(ps_lc < panelid_table(tileid)) plc_count = plc_count + 1
      end if 
    end do
    end do
 
    allocate( this%rcdomIJKP2LCMeshID(ilc_count,jlc_count,klc_count,plc_count) )
    do plc=1, plc_count
    do klc=1, klc_count
    do jlc=1, jlc_count
    do ilc=1, ilc_count
      this%rcdomIJKP2LCMeshID(ilc,jlc,klc,plc) = ilc + (jlc - 1)*ilc_count + (klc - 1)*ilc_count*jlc_count &
                                               + (plc-1)*ilc_count*jlc_count*klc_count
    end do
    end do
    end do
    end do
 
    return
  end subroutine messhcubedspheredom3d_assigndomid
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_coord_conv( x, y, z, xX, xY, xZ, yX, yY, yZ, zX, zY, zZ, &
    vx, vy, vz, elem )
 
    implicit none
 
    type(elementbase3d), intent(in) :: elem
    real(RP), intent(out) :: x(elem%Np), y(elem%Np), z(elem%Np)
    real(RP), intent(out) :: xX(elem%Np), xY(elem%Np), xZ(elem%Np)
    real(RP), intent(out) :: yX(elem%Np), yY(elem%Np), yZ(elem%Np)
    real(RP), intent(out) :: zX(elem%Np), zY(elem%Np), zZ(elem%Np)
    real(RP), intent(in) :: vx(elem%Nv), vy(elem%Nv), vz(elem%Nv)
 
    !-------------------------------------------------
    
    x(:) = vx(1) + 0.5_rp*(elem%x1(:) + 1.0_rp)*(vx(2) - vx(1))
    y(:) = vy(1) + 0.5_rp*(elem%x2(:) + 1.0_rp)*(vy(3) - vy(1))
    z(:) = vz(1) + 0.5_rp*(elem%x3(:) + 1.0_rp)*(vz(5) - vz(1))
 
    xx(:) = 0.5_rp*(vx(2) - vx(1)) !matmul(refElem%Dx1,mesh%x1(:,n))
    xy(:) = 0.0_rp                 !matmul(refElem%Dx2,mesh%x1(:,n))
    xz(:) = 0.0_rp                 !matmul(refElem%Dx3,mesh%x1(:,n))
    yx(:) = 0.0_rp                 !matmul(refElem%Dx1,mesh%x2(:,n))
    yy(:) = 0.5_rp*(vy(3) - vy(1)) !matmul(refElem%Dx2,mesh%x2(:,n))
    yz(:) = 0.0_rp                 !matmul(refElem%Dx3,mesh%x2(:,n))
    zx(:) = 0.0_rp                 !matmul(refElem%Dx1,mesh%x3(:,n))
    zy(:) = 0.0_rp                 !matmul(refElem%Dx2,mesh%x3(:,n))
    zz(:) = 0.5_rp*(vz(5) - vz(1)) !matmul(refElem%Dx3,mesh%x3(:,n))
 
    return
  end subroutine meshcubedspheredom3d_coord_conv  
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_calc_normal( normal_fn, &
    Escale_f, fid_h, fid_v, elem )
 
    implicit none
 
    type(elementbase3d), intent(in) :: elem
    real(RP), intent(out) :: normal_fn(elem%NfpTot,3)
    integer, intent(in) :: fid_h(elem%Nfp_h,elem%Nfaces_h)
    integer, intent(in) :: fid_v(elem%Nfp_v,elem%Nfaces_v)        
    real(RP), intent(in) :: Escale_f(elem%NfpTot,3,3)
 
    integer :: d
    !-------------------------------------------------
 
    do d=1, 3
      normal_fn(fid_h(:,1),d) = - escale_f(fid_h(:,1),2,d)
      normal_fn(fid_h(:,2),d) = + escale_f(fid_h(:,2),1,d)
      normal_fn(fid_h(:,3),d) = + escale_f(fid_h(:,3),2,d)
      normal_fn(fid_h(:,4),d) = - escale_f(fid_h(:,4),1,d)
 
      normal_fn(fid_v(:,1),d) = - escale_f(fid_v(:,1),3,d)
      normal_fn(fid_v(:,2),d) = + escale_f(fid_v(:,2),3,d)    
    end do
 
    return
  end subroutine meshcubedspheredom3d_calc_normal 
 
  !--
 
!OCL SERIAL
  subroutine meshcubedspheredom3d_set_metric( this )
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_cs2lonlatpos, &
      cubedspherecoordcnv_getmetric
        
    implicit none
    class(meshcubedspheredom3d), intent(inout), target :: this
 
    integer :: n
    integer :: ke, ke2D
 
    class(localmesh3d), pointer :: lcmesh
    class(localmesh2d), pointer :: lcmesh2D
    class(elementbase2d), pointer :: elem2D
 
    real(RP), allocatable :: gam2D(:,:)
    !----------------------------------------------------
 
    do n=1, this%mesh2D%LOCAL_MESH_NUM
      lcmesh => this%lcmesh_list(n)
      lcmesh2d => this%mesh2D%lcmesh_list(n)
      elem2d => lcmesh2d%refElem2D
 
      allocate( gam2d(elem2d%Np,lcmesh2d%Ne) )
      gam2d(:,:) = 1.0_rp
 
      call cubedspherecoordcnv_cs2lonlatpos( &
        lcmesh2d%panelID, lcmesh2d%pos_en(:,:,1), lcmesh2d%pos_en(:,:,2), gam2d(:,:), & ! (in)
        lcmesh2d%Ne * elem2d%Np,                                                      & ! (in)
        lcmesh%lon2D(:,:), lcmesh%lat2D(:,:)                                          ) ! (out)
 
      call cubedspherecoordcnv_getmetric( &
        lcmesh2d%pos_en(:,:,1), lcmesh2d%pos_en(:,:,2), elem2d%Np * lcmesh2d%Ne, this%RPlanet, & ! (in)
        lcmesh%G_ij, lcmesh%GIJ, lcmesh%GsqrtH                                                 ) ! (out)
 
      !$omp parallel do private(ke2D)
      do ke=lcmesh%NeS, lcmesh%NeE
        ke2d = lcmesh%EMap3Dto2D(ke)
 
        if ( this%shallow_approx ) then
          lcmesh%gam(:,ke) = 1.0_rp
        else
          lcmesh%gam(:,ke) = 1.0_rp + lcmesh%pos_en(:,ke,3) / this%RPlanet
        end if
        lcmesh%Gsqrt(:,ke) = lcmesh%GsqrtH(lcmesh%refElem3D%IndexH2Dto3D(:),ke2d)        
      end do
 
      call fill_halo_metric( lcmesh%Gsqrt, lcmesh%gam, lcmesh%VMapM, lcmesh%VMapP, lcmesh, lcmesh%refElem3D )
 
      deallocate( gam2d )
      !--
    end do
 
 
    return
  end subroutine meshcubedspheredom3d_set_metric
 
!OCL SERIAL
  subroutine fill_halo_metric( Gsqrt, gam, vmapM, vmapP, lmesh, elem )
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: vmapM(elem%NfpTot*lmesh%Ne)
    integer, intent(in) :: vmapP(elem%NfpTot*lmesh%Ne)
    real(RP), intent(inout) :: Gsqrt(elem%Np*lmesh%NeA)
    real(RP), intent(inout) :: gam(elem%Np*lmesh%NeA)
 
    integer :: i, iM, iP
    !------------------------------------------------
 
    !$omp parallel do private(i, iM, iP)
    do i=1, elem%NfpTot*lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
      if ( ip > elem%Np * lmesh%Ne ) then
        gsqrt(ip) = gsqrt(im)
        gam(ip) = gam(im)
      end if
    end do  
    return
  end subroutine fill_halo_metric
 
end module scale_mesh_cubedspheredom3d