mod_ocean_mesh_gm.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Ocean / Mesh
!!
!! @par Description
!!          Module for mesh with ocean global model
!!
!! @author Yuta kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_ocean_mesh_gm
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
 
  use scale_meshfield_base, only: meshfield3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_mesh_cubedspheredom3d, only: meshcubedspheredom3d
  use scale_element_base, only: elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_3d, only: localmesh3d
  use scale_meshfieldcomm_cubedspheredom3d, only: meshfieldcommcubedspheredom3d
  use scale_sparsemat, only: sparsemat
  use scale_file_restart_meshfield, only: file_restart_meshfield_component
 
  use scale_model_var_manager, only: modelvarmanager
  use scale_model_mesh_manager, only: modelmesh3d
 
  use mod_ocean_mesh, only: &
    oceanmesh, ocn_mesh_max_commnuicator_num
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedures
  !
 
  !> Derived type to manage a computational mesh of global ocean model
  !!
  type, extends(oceanmesh), public :: oceanmeshgm
    type(meshcubedspheredom3d) :: mesh                                               !< Object for 3D cubed-sphere mesh
    type(meshfieldcommcubedspheredom3d) :: comm_list(ocn_mesh_max_commnuicator_num)
  contains
    procedure :: init => oceanmeshgm_init
    procedure :: final => oceanmeshgm_final
    procedure :: create_communicator => oceanmeshgm_create_communicator
    procedure :: setup_restartfile1 => oceanmeshgm_setup_restartfile1
    procedure :: setup_restartfile2 => oceanmeshgm_setup_restartfile2
    procedure :: calc_uvmet => oceanmeshgm_calc_uvmet
    procedure :: setup_vcoordinate => oceanmeshgm_setup_vcoordinate
  end type oceanmeshgm
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-------------------
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
 
  !> Initialize an object to manage computational mesh
  !!
!OCL SERIAL
  subroutine oceanmeshgm_init( this )    
    use scale_const, only: &
      rplanet => const_radius
    use scale_mesh_base2d, only: meshbase2d
    use scale_meshutil_vcoord, only: &
      meshutil_get_vcoord_typeid
    
    implicit none
    class(oceanmeshgm), target, intent(inout) :: this
 
    real(RP) :: dom_zmin          = 0.0_rp      !< Minimum vertical coordinate value of the computational domain
    real(RP) :: dom_zmax          = 10.0e3_rp   !< Maximum vertical coordinate value of the computational domain
    logical  :: isPeriodicZ       = .false.     !< Flag whether a periodic boundary condition is applied in the vertical direction
  
    integer, parameter :: FZ_nmax = 1000
    real(RP) :: FZ(FZ_nmax)                  !< Values of the vertically computational coordinate at the element boundaries
 
    !* Global
    logical :: SHALLOW_OCN_APPROX_FLAG = .true.  !< Flag whether the shallow ocean approximation is applied
    integer  :: NeGX               = 2           !< Number of finite element in the y-coordinate direction in each panel of the cubed-sphere mesh
    integer  :: NeGY               = 2           !< Number of finite element in the y-coordinate direction in each panel of the cubed-sphere mesh
    integer  :: NeZ                = 2           !< Number of finite element in the vertical direction in each MPI process
    integer  :: NLocalMeshPerPrc   = 6           !< Number of local mesh per MPI process 
    integer  :: Nprc               = 1           !< Total number of MPI process
    integer  :: PolyOrder_h        = 2           !< Polynomial order for the horizontal direction
    integer  :: PolyOrder_v        = 2           !< Polynomial order for the z-direction
    logical  :: LumpedMassMatFlag  = .false.     !< Flag whether a mass lumping is applied
 
    character(len=H_LONG) :: TOPO_IN_BASENAME    = ''                   !< Basename of the input file
    character(len=H_MID)  :: TOPO_IN_VARNAME     = 'topo'               !< Variable name of topography in the input file
    character(len=H_MID)  :: VERTICAL_COORD_NAME = "TERRAIN_FOLLOWING"  !< Type of the vertical coordinate
 
    logical :: COMM_USE_MPI_PC             = .false. !< Flag whether persistent communication is used in MPI
    logical :: COMM_USE_MPI_PC_FUJITSU_EXT
 
    character(len=H_SHORT) :: Element_operation_type = 'General' !< General or TensorProd3D
    character(len=H_SHORT) :: SpMV_storage_format    = 'ELL'     !< CSR or ELL
    
    namelist / param_ocean_mesh / &
      shallow_ocn_approx_flag,                     &
      dom_zmin, dom_zmax,                          &
      fz, isperiodicz,                             &
      negx, negy, nez, nlocalmeshperprc, nprc,     &
      polyorder_h, polyorder_v, lumpedmassmatflag, &
      element_operation_type,                      &
      spmv_storage_format,                         &
      vertical_coord_name,                         &
      topo_in_basename, topo_in_varname,           &
      comm_use_mpi_pc,                             &
      comm_use_mpi_pc_fujitsu_ext
    
    integer :: k
    logical :: is_spec_FZ
    
    integer :: ierr
    !-------------------------------------------
 
    log_newline
    log_info("OCN_MESH_setup",*) 'Setup'
 
    fz(:) = -1.0_rp
#ifdef __FUJITSU
    comm_use_mpi_pc_fujitsu_ext = .true.
#else
    comm_use_mpi_pc_fujitsu_ext = .false.
#endif
 
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_ocean_mesh,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("OCN_MESH_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("OCN_MESH_setup",*) 'Not appropriate names in namelist PARAM_OCEAN_MESH. Check!'
      call prc_abort
    endif
    log_nml(param_ocean_mesh)
 
    !----
 
    !- Setup the element
 
    call this%element%Init( polyorder_h, polyorder_v, lumpedmassmatflag )
    call this%element_v1D%Init( polyorder_v, lumpedmassmatflag )
 
    !- Setup the mesh
 
    is_spec_fz = .true.
    do k=1, nez+1
      if (fz(k) < 0.0_rp) then
        is_spec_fz = .false.
      end if
    end do
    if (is_spec_fz) then
      call this%mesh%Init( &
        negx, negy, nez, rplanet, dom_zmin, dom_zmax,          &
        this%element, nlocalmeshperprc, nproc=nprc,            &
        fz=fz(1:nez+1), shallow_approx=shallow_ocn_approx_flag )
    else
      call this%mesh%Init( &
        negx, negy, nez, rplanet, dom_zmin, dom_zmax, &
        this%element, nlocalmeshperprc, nproc=nprc,   &
        shallow_approx=shallow_ocn_approx_flag        )
    end if
    
    call this%mesh%Generate()
    
    !-
 
    call this%OceanMesh_Init( this%mesh )
    call this%PrepairElementOperation( element_operation_type, spmv_storage_format )
 
    !- Set topography & vertical coordinate
    
    call this%Read_topography_file( topo_in_basename, topo_in_varname, &
      this%mesh%mesh2D, dom_zmin )
 
    this%vcoord_type_id = meshutil_get_vcoord_typeid( vertical_coord_name )
    ! We tentatively disable the setup of vertical coordinate in ocean component. 
    ! For future, we will implement a more flexible setup of vertical coordinate.
    ! call this%Setup_vcoordinate()
 
    !-
    this%comm_use_mpi_pc             = comm_use_mpi_pc
    this%comm_use_mpi_pc_fujitsu_ext = comm_use_mpi_pc_fujitsu_ext
 
    return
  end subroutine oceanmeshgm_init
 
  !> Finalize an object to manage computational mesh
  !!
!OCL SERIAL
  subroutine oceanmeshgm_final(this)
    implicit none
 
    class(oceanmeshgm), intent(inout) :: this
    integer :: commid
    !-------------------------------------------
 
    do commid=1, this%communicator_num
      call this%comm_list(commid)%Final()
    end do
 
    call this%mesh%Final()
    call this%OceanMesh_Final()
    
    return
  end subroutine oceanmeshgm_final
 
  !> Create a communicator for data communication on global computational domain
  !!
!OCL SERIAL
  subroutine oceanmeshgm_create_communicator( this, sfield_num, hvfield_num, htensorfield_num, &
    var_manager, field_list, commid, &
    field_list_is, field_list_ie     )
    implicit none
    class(oceanmeshgm), target, intent(inout) :: this
    integer, intent(in) :: sfield_num        !< Number of scalar fields
    integer, intent(in) :: hvfield_num       !< Number of horizontal vector fields
    integer, intent(in) :: htensorfield_num  !< Number of horizontal tensor fields
    class(modelvarmanager), intent(inout) :: var_manager !< Object to manage variables
    class(meshfield3d), intent(in) :: field_list(:)      !< Array of 3D fields
    integer, intent(out) :: commid           !< Communicator ID
    integer, intent(in), optional :: field_list_is  !< Starting index of the field list for communication (optional)
    integer, intent(in), optional :: field_list_ie  !< Ending index of the field list for communication (optional)
    !-----------------------------------------------------
 
    commid = this%Get_communicatorID( ocn_mesh_max_commnuicator_num )
    call this%comm_list(commid)%Init( sfield_num, hvfield_num, htensorfield_num, this%mesh )
    if ( this%comm_use_mpi_pc ) then
      call this%comm_list(commid)%Prepare_PC( this%comm_use_mpi_pc_fujitsu_ext )
    end if
    call var_manager%MeshFieldComm_Prepare( this%comm_list(commid), field_list, field_list_is, field_list_ie )
 
    return
  end subroutine oceanmeshgm_create_communicator
 
!OCL SERIAL
  subroutine oceanmeshgm_setup_restartfile1( this, restart_file, var_num, dim_name_postfix )
    implicit none
    class(oceanmeshgm), target, intent(inout) :: this
    class(file_restart_meshfield_component), intent(inout) :: restart_file
    integer, intent(in) :: var_num
    character(*), intent(in) :: dim_name_postfix
    !------------------------------------------------
 
    call restart_file%Init( 'OCEAN', var_num, dim_name_postfix, meshcubedsphere3d=this%mesh )
    return
  end subroutine oceanmeshgm_setup_restartfile1
 
!OCL SERIAL
  subroutine oceanmeshgm_setup_restartfile2( this, restart_file, &
    in_basename, in_postfix_timelabel,                           &
    out_basename, out_postfix_timelabel,                         &
    out_dtype, out_title, var_num, dim_name_postfix              )
    implicit none
    class(oceanmeshgm), target, intent(inout) :: this
    class(file_restart_meshfield_component), intent(inout) :: restart_file
    character(*), intent(in) :: in_basename
    logical, intent(in) :: in_postfix_timelabel
    character(*), intent(in) :: out_basename
    logical, intent(in) :: out_postfix_timelabel
    character(*), intent(in) :: out_title
    character(*), intent(in) :: out_dtype  
    integer, intent(in) :: var_num  
    character(*), intent(in) :: dim_name_postfix
    !-----------------------------------------------------------
 
    call restart_file%Init( 'OCEAN', in_basename, in_postfix_timelabel, &
      out_basename, out_postfix_timelabel, out_dtype, out_title,        &
      var_num, dim_name_postfix,                                        &
      meshcubedsphere3d=this%mesh )
 
  end subroutine oceanmeshgm_setup_restartfile2
 
!> Calculate horizontal vector components in longitude-latitude coordinates
!OCL SERIAL
  subroutine oceanmeshgm_calc_uvmet( this, U, V, &
    Umet, Vmet )
 
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_cs2lonlatvec
    implicit none
    class(oceanmeshgm), target, intent(in) :: this
    type(meshfield3d), intent(in) :: U
    type(meshfield3d), intent(in) :: V
    type(meshfield3d), intent(inout) :: Umet
    type(meshfield3d), intent(inout) :: Vmet
 
    integer :: n
    integer :: ke, ke2D
    type(localmesh3d), pointer :: lcmesh
    class(elementbase3d), pointer :: elem
    !------------------------------------------
 
    do n=1, this%mesh%LOCAL_MESH_NUM
      lcmesh => this%mesh%lcmesh_list(n)
      elem => lcmesh%refElem3D
      call cubedspherecoordcnv_cs2lonlatvec( &
        lcmesh%panelID, lcmesh%pos_en(:,:,1), lcmesh%pos_en(:,:,2), & ! (in)
        lcmesh%gam, elem%Np * lcmesh%Ne,                            & ! (in)
        u%local(n)%val(:,lcmesh%NeS:lcmesh%NeE),                    & ! (in)
        v%local(n)%val(:,lcmesh%NeS:lcmesh%NeE),                    & ! (in)
        umet%local(n)%val(:,lcmesh%NeS:lcmesh%NeE),                 & ! (out)
        vmet%local(n)%val(:,lcmesh%NeS:lcmesh%NeE)                  ) ! (out)
    end do
 
    return
  end subroutine oceanmeshgm_calc_uvmet
 
!> Setup the vertical coordinate
!OCL SERIAL
  subroutine oceanmeshgm_setup_vcoordinate( this )
    use scale_meshfieldcomm_cubedspheredom2d, only: meshfieldcommcubedspheredom2d
    implicit none
    class(oceanmeshgm), TARGET, INTENT(INOUT) :: this
 
    type(meshfieldcommcubedspheredom3d) :: comm3D
    type(meshfieldcommcubedspheredom2d) :: comm2D
    !-------------------------------------------------
 
    call comm2d%Init( 1, 0, 0, this%mesh%mesh2D )
    call comm3d%Init( 2, 1, 0, this%mesh )
 
    call this%topography%SetVCoordinate( this%ptr_mesh,   &
      this%vcoord_type_id, this%mesh%zmax_gl, comm3d, comm2d )
 
    call comm2d%Final()
    call comm3d%Final()
 
    return
  end subroutine oceanmeshgm_setup_vcoordinate
 
end module mod_ocean_mesh_gm