mod_ocean_mesh_rm.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Ocean / Mesh
!!
!! @par Description
!!          Module for mesh with ocean regional model
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_ocean_mesh_rm
  !-----------------------------------------------------------------------------
  !
  !++ 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_cubedom3d, only: meshcubedom3d
  use scale_element_base, only: elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_3d, only: localmesh3d
  use scale_meshfieldcomm_cubedom3d, only: meshfieldcommcubedom3d
  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 regional ocean model
  !!
  type, extends(oceanmesh), public :: oceanmeshrm
    type(meshcubedom3d) :: mesh                                               !< Object to manage a cubic computational mesh
    type(meshfieldcommcubedom3d) :: comm_list(ocn_mesh_max_commnuicator_num)  !< Object to manage data communication on a cubic computational mesh
  contains
    procedure :: init => oceanmeshrm_init
    procedure :: final => oceanmeshrm_final
    procedure :: create_communicator => oceanmeshrm_create_communicator
    procedure :: setup_restartfile1 => oceanmeshrm_setup_restartfile1
    procedure :: setup_restartfile2 => oceanmeshrm_setup_restartfile2
    procedure :: calc_uvmet => oceanmeshrm_calc_uvmet
    procedure :: setup_vcoordinate => oceanmeshrm_setup_vcoordinate
  end type oceanmeshrm
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-------------------
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
 
  !- Ocean Mesh RM -----------------------------------------
 
  !> Initialize an object to manage computational mesh
  !!
  subroutine oceanmeshrm_init( this )
    use scale_meshutil_vcoord, only: &
      meshutil_get_vcoord_typeid
 
    implicit none
    class(oceanmeshrm), target, intent(inout) :: this
 
    real(RP) :: dom_xmin         = 0.0_rp      !< Minimum x-coordinate value of the computational domain
    real(RP) :: dom_xmax         = 100.0e3_rp  !< Maximum x-coordinate value of the computational domain
    real(RP) :: dom_ymin         = 0.0_rp      !< Minimum y-coordinate value of the computational domain
    real(RP) :: dom_ymax         = 100.0e3_rp  !< Maximum y-coordinate value of the computational domain
    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  :: isPeriodicX       = .true.     !< Flag whether a periodic boundary condition is applied in the x-direction
    logical  :: isPeriodicY       = .true.     !< Flag whether a periodic boundary condition is applied in the y-direction
    logical  :: isPeriodicZ       = .false.    !< Flag whether a periodic boundary condition is applied in the vertical direction
 
    integer  :: NeX               = 2        !< Number of finite element in the x-direction in each MPI process
    integer  :: NeY               = 2        !< Number of finite element in the y-direction in each MPI process
    integer  :: NeZ               = 2        !< Number of finite element in the vertical direction in each MPI process
    integer  :: NprcX             = 1        !< Number of MPI process in the x-direction
    integer  :: NprcY             = 1        !< Number of MPI process in the y-direction
    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'     !< Storage format of sparse matrix (CSR or ELL)
 
    integer, parameter :: ATMOS_MESH_NLocalMeshPerPrc = 1
 
    integer, parameter :: FZ_nmax = 1000
    real(RP) :: FZ(FZ_nmax)                  !< Values of the vertically computational coordinate at the element boundaries
 
    namelist / param_ocean_mesh / &
      dom_xmin, dom_xmax,                          &
      dom_ymin, dom_ymax,                          &
      dom_zmin, dom_zmax,                          &
      fz,                                          &
      isperiodicx, isperiodicy, isperiodicz,       &
      nex, ney, nez,                               &
      polyorder_h, polyorder_v, lumpedmassmatflag, &
      element_operation_type,                      &
      spmv_storage_format,                         &
      nprcx, nprcy,                                &
      topo_in_basename, topo_in_varname,           &
      vertical_coord_name,                         &
      comm_use_mpi_pc,                             &
      comm_use_mpi_pc_fujitsu_ext
    
    integer :: k
    logical :: is_spec_FZ
 
    integer :: ierr
    !-------------------------------------------
 
    log_newline
    log_info("OCEAN_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("OCEAN_MESH_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
        log_error("OCEAN_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( &
        nprcx*nex, nprcy*ney, nez,                                 &
        dom_xmin, dom_xmax,dom_ymin, dom_ymax, dom_zmin, dom_zmax, &
        isperiodicx, isperiodicy, isperiodicz,                     &
        this%element, atmos_mesh_nlocalmeshperprc, nprcx, nprcy,   &
        fz=fz(1:nez+1)    )
    else
      call this%mesh%Init( &
        nprcx*nex, nprcy*ney, nez,                                 &
        dom_xmin, dom_xmax,dom_ymin, dom_ymax, dom_zmin, dom_zmax, &
        isperiodicx, isperiodicy, isperiodicz,                     &
        this%element, atmos_mesh_nlocalmeshperprc, nprcx, nprcy    )
    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 oceanmeshrm_init
 
  !> Finalize an object to manage computational mesh
  !!
  subroutine oceanmeshrm_final(this)
    implicit none
 
    class(oceanmeshrm), 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 oceanmeshrm_final
 
  !> Create a communicator for data communication on regional computational domain
  !!
  subroutine oceanmeshrm_create_communicator( this, sfield_num, hvfield_num, htensorfield_num, &
    var_manager, field_list, commid, &
    field_list_is, field_list_ie     )
    implicit none
    class(oceanmeshrm), 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 oceanmeshrm_create_communicator 
  
  subroutine oceanmeshrm_setup_restartfile1( this, restart_file, var_num, dim_name_postfix )
    implicit none
    class(oceanmeshrm), 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, mesh3d=this%mesh )
    return
  end subroutine oceanmeshrm_setup_restartfile1
 
  subroutine oceanmeshrm_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(oceanmeshrm), 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,                                       &
      mesh3d=this%mesh )
 
    return
  end subroutine oceanmeshrm_setup_restartfile2
 
  subroutine oceanmeshrm_calc_uvmet( this, U, V, &
    Umet, Vmet )
    implicit none
    class(oceanmeshrm), 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
    type(localmesh3d), pointer :: lcmesh
    !------------------------------------------
 
    do n=1, this%mesh%LOCAL_MESH_NUM
      lcmesh => this%mesh%lcmesh_list(n)
      !$omp parallel do
      do ke=lcmesh%NeS, lcmesh%NeE
        umet%local(n)%val(:,ke) = u%local(n)%val(:,ke)
        vmet%local(n)%val(:,ke) = v%local(n)%val(:,ke)
      end do
    end do
 
    return
  end subroutine oceanmeshrm_calc_uvmet
 
!> Setup the vertical coordinate
  subroutine oceanmeshrm_setup_vcoordinate( this )
    use scale_meshfieldcomm_rectdom2d, only: meshfieldcommrectdom2d
    implicit none
    class(oceanmeshrm), TARGET, INTENT(INOUT) :: this
 
    type(meshfieldcommcubedom3d) :: comm3D
    type(meshfieldcommrectdom2d) :: comm2D
    !-------------------------------------------------
 
    call comm2d%Init( 1, 0, 0, this%mesh%mesh2D )
    call comm3d%Init( 2, 1, 0, this%mesh )
 
    call this%topography%SetVCoordinate( this%mesh,   &
      this%vcoord_type_id, this%mesh%zmax_gl, comm3d, comm2d )
 
    call comm2d%Final()
    call comm3d%Final()
 
    return
  end subroutine oceanmeshrm_setup_vcoordinate
 
end module mod_ocean_mesh_rm