scale_mesh_topography.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_mesh_topography
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_localmesh_3d, only: localmesh3d
  use scale_localmesh_2d, only: localmesh2d
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
  use scale_element_base, only: elementbase3d
  use scale_meshfield_base, only: &
    meshfield2d, meshfield3d
  use scale_meshfieldcomm_base, only: &
    meshfieldcommbase, meshfieldcontainer
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  type, public :: meshtopography
    type(meshfield2d) :: topo 
  contains
    procedure :: init => meshtopography_init
    procedure :: final => meshtopography_final
    procedure :: setvcoordinate => meshtopography_set_vcoordinate
  end type meshtopography
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
!OCL SERIAL
  subroutine meshtopography_init( this, varname, mesh )
    implicit none
 
    class(meshtopography), intent(inout) :: this
    character(len=*), intent(in) :: varname
    class(meshbase2d), intent(in), target :: mesh
 
    integer :: n
    !-----------------------------------------------------------------------------
 
    call this%topo%Init( varname, "m", mesh )
 
    do n=1, mesh%LOCAL_MESH_NUM
!$omp parallel workshare
      this%topo%local(n)%val(:,:) = 0.0_rp
!$omp end parallel workshare
    end do
 
    return
  end subroutine meshtopography_init
 
!OCL SERIAL
  subroutine meshtopography_final( this )
    implicit none
 
    class(meshtopography), intent(inout) :: this
    !-----------------------------------------------------------------------------
          
    call this%topo%Final()
 
    return
  end subroutine meshtopography_final
  
!OCL SERIAL
  subroutine meshtopography_set_vcoordinate( this, mesh3D, &
      vcoord_id, zTop, comm3D, comm2D                      )
    
    use scale_sparsemat, only: sparsemat
    use scale_meshutil_vcoord, only: meshutil_vcoord_getmetric
    use scale_meshfieldcomm_cubedspheredom3d, only: meshfieldcommcubedspheredom3d
    implicit none
 
    class(meshtopography), target, intent(inout) :: this
    class(meshbase3d), intent(inout), target :: mesh3D
    integer, intent(in) :: vcoord_id
    real(RP), intent(in) :: zTop
    class(meshfieldcommbase), intent(inout) :: comm3D
    class(meshfieldcommbase), intent(inout) :: comm2D
 
    type(sparsemat) :: Dx2D, Dy2D, Lift2D
    class(localmesh3d), pointer :: lcmesh
    class(localmesh2d), pointer :: lcmesh2D
    class(elementbase3d), pointer :: elem3D
 
    integer :: n
    integer :: ke, ke2D
 
    type(meshfieldcontainer) :: comm2d_varlist(1)
    type(meshfieldcontainer) :: comm3d_varlist(4)
 
    type(meshfield3d), target :: zlev, GsqrtV, G13, G23
    type(meshfield3d), target :: tmp_G13, tmp_G23
 
    logical :: flag_covariantvec
    !-------------------------------------------------------------
 
    lcmesh => mesh3d%lcmesh_list(1)
    lcmesh2d => lcmesh%lcmesh2D
    call dx2d  %Init( lcmesh2d%refElem2D%Dx1  )
    call dy2d  %Init( lcmesh2d%refElem2D%Dx2  )
    call lift2d%Init( lcmesh2d%refElem2D%Lift )
 
    ! Exchange topography data to fill halo
 
    comm2d_varlist(1)%field2d => this%topo
    call comm2d%Put(comm2d_varlist, 1)
    call comm2d%Exchange()
    call comm2d%Get(comm2d_varlist, 1)
 
    ! Calculate metric factors associated with general vertical coordinates
 
    call zlev%Init( "zlev", "", mesh3d )
    call gsqrtv%Init( "GsqrtV", "", mesh3d )
    call g13%Init( "G13", "", mesh3d )
    call g23%Init( "G23", "", mesh3d )
 
    call tmp_g13%Init( "tmp_G13", "", mesh3d )
    call tmp_g23%Init( "tmp_G23", "", mesh3d )
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh => mesh3d%lcmesh_list(n)
      lcmesh2d => lcmesh%lcmesh2D
      elem3d => lcmesh%refElem3D
 
      call meshutil_vcoord_getmetric( &
        g13%local(n)%val, g23%local(n)%val, zlev%local(n)%val,  & ! (out)
        gsqrtv%local(n)%val,                                    & ! (out)
        this%topo%local(n)%val(:,:), ztop, vcoord_id,           & ! (in)
        lcmesh, lcmesh%refElem3D, lcmesh2d, lcmesh2d%refElem2D, & ! (in)
        dx2d, dy2d, lift2d                                      ) ! (in)
 
      !$omp parallel do private(ke2D)
      do ke=lcmesh%NeS, lcmesh%NeE
        ke2d = lcmesh%EMap3Dto2D(ke)
        tmp_g13%local(n)%val(:,ke) = &
          lcmesh%GIJ(elem3d%IndexH2Dto3D(:),ke2d,1,1) * g13%local(n)%val(:,ke) &
        + lcmesh%GIJ(elem3d%IndexH2Dto3D(:),ke2d,1,2) * g23%local(n)%val(:,ke)
        tmp_g23%local(n)%val(:,ke) = &
          lcmesh%GIJ(elem3d%IndexH2Dto3D(:),ke2d,2,1) * g13%local(n)%val(:,ke) &
        + lcmesh%GIJ(elem3d%IndexH2Dto3D(:),ke2d,2,2) * g23%local(n)%val(:,ke)
      end do        
    end do
 
    ! Exchange metric data to fill halo
 
    comm3d_varlist(1)%field3d => zlev
    comm3d_varlist(2)%field3d => gsqrtv
    comm3d_varlist(3)%field3d => tmp_g13
    comm3d_varlist(4)%field3d => tmp_g23
 
    flag_covariantvec = .false.
    select type(comm3d)
    type is (meshfieldcommcubedspheredom3d)
      call comm3d%SetCovariantVec( 1, g13, g23 )
      flag_covariantvec = .true.
    end select
 
    call comm3d%Put(comm3d_varlist, 1)
    call comm3d%Exchange()
    call comm3d%Get(comm3d_varlist, 1)
 
    ! Update metric data managed by local mesh
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh => mesh3d%lcmesh_list(n)
      elem3d => lcmesh%refElem3D
 
      if ( flag_covariantvec ) then
        call mesh3d%Set_geometric_with_vcoord( n, &
          gsqrtv%local(n)%val, zlev%local(n)%val, g13%local(n)%val, g23%local(n)%val )
      else
        call mesh3d%Set_geometric_with_vcoord( n, &
          gsqrtv%local(n)%val, zlev%local(n)%val, tmp_g13%local(n)%val, tmp_g23%local(n)%val )
      end if
    end do
 
    !---
    call zlev%Final()
    call gsqrtv%Final()
    call g13%Final()
    call g23%Final()
    call tmp_g13%Final()
    call tmp_g23%Final()
 
    call dx2d%Final()
    call dy2d%Final()
    call lift2d%Final()
 
    return
  end subroutine meshtopography_set_vcoordinate
 
end module scale_mesh_topography