mod_ocean_dyn.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Ocean Dynamics
!!
!! @par Description
!!          Module for oceanic dynamical process
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_ocean_dyn
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !  
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prof
  use scale_const, only: &
    undef8 => const_undef8
 
  use scale_sparsemat, only: sparsemat
  use scale_timeint_rk, only: timeint_rk
  
  use scale_mesh_base, only: meshbase
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
 
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_element_base, only: &
    elementbase, elementbase2d, elementbase3d
    
  use scale_meshfield_base, only: meshfieldbase
  use scale_localmeshfield_base, only: &
    localmeshfieldbase, localmeshfieldbaselist
 
  use scale_model_mesh_manager, only: modelmeshbase
  use scale_model_var_manager, only: modelvarmanager
  use scale_model_component_proc, only:  modelcomponentproc
 
!  use mod_ocean_dyn_dgm_driver, only: OceanDynDGMDriver_hydro3d
  
  use mod_ocean_mesh, only: oceanmesh
  use mod_ocean_vars, only: &
    oceanvars_getlocalmeshprgvars
  use mod_ocean_dyn_vars, only: &
    oceandynvars,                      &
    oceandynauxvars_getlocalmeshfields
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  !> Derived type to manage a component of oceanic dynamics
  !!
  type, extends(modelcomponentproc), public :: oceandyn
!    type(OceanDynDGMDriver_hydro3d) :: dyncore_driver !< A driver object to manage a oceanic dynamical core
 
    type(oceandynvars) :: dyn_vars                     !< An object to manage variables in a component of oceanic dynamics
    integer :: eqs_type                                !< Type of governing equations: 0=NONE, 1=SLAB
 
    class(meshbase3d), pointer :: mesh3d
    real(rp) :: dtsec                                  !< Timestep for a oceanic dynamical core
  contains
    procedure, public :: setup => oceandyn_setup 
    procedure, public :: calc_tendency => oceandyn_calc_tendency
    procedure, public :: update => oceandyn_update
    procedure, public :: finalize => oceandyn_finalize
  end type oceandyn
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !-----------------------------------------------------------------------------
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: setup_coriolis_parameter
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
 
  integer :: ocean_dyn_eqs_fixed_init_temp = 0
  integer :: ocean_dyn_eqs_slab            = 1
 
contains
 
!> Setup an object to manage a component of oceanic dynamics
!!
!! @param model_mesh Object to manage computational mesh of oceanic model
!! @param tm_parent_comp Object to mange a temporal scheme in a parent component
!!
!OCL SERIAL
  subroutine oceandyn_setup( this, model_mesh, tm_parent_comp )
    use mod_ocean_mesh, only: oceanmesh
    use scale_time_manager, only: time_manager_component
 
    implicit none
 
    class(oceandyn), intent(inout) :: this
    class(modelmeshbase), target, intent(in) :: model_mesh
    class(time_manager_component), intent(inout) :: tm_parent_comp
 
    character(len=H_MID) :: eqs_type             = "SLAB"            !< Type of governing equations: FIXED_TEMP or SLAB
    character(len=H_SHORT) :: tinteg_type        = 'ERK_SSP_3s3o'    !< Type of temporal scheme for a dry dynamical core
    character(len=H_SHORT) :: tinteg_type_tracer = 'ERK_SSP_3s3o'    !< Type of temporal scheme for tracer advection equations
    real(dp) :: time_dt                          = undef8            !< Timestep for a atmospheric dynamical core
    character(len=H_SHORT) :: time_dt_unit       = 'SEC'             !< Unit of timestep
    
!    logical :: MODALFILTER_FLAG           = .false. !< Flag to set whether a modal filtering is used
    
    namelist / param_ocean_dyn /       &
      eqs_type,                        &
      tinteg_type,                     &
      tinteg_type_tracer,              &      
      time_dt,                         &
      time_dt_unit                    
!      MODALFILTER_FLAG,                
      
    class(oceanmesh), pointer     :: ocean_mesh
    class(meshbase), pointer      :: ptr_mesh
    class(localmeshbase), pointer :: ptr_lcmesh
    class(elementbase3d), pointer :: elem3d
    integer :: n
    real(dp) :: dtsec
 
    integer :: ierr
    !--------------------------------------------------
 
    if (.not. this%IsActivated()) return
    log_info('OceanDyn_setup',*)
 
    eqs_type = 'SLAB'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_ocean_dyn,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("OceanDyn_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("OceanDyn_setup",*) 'Not appropriate names in namelist PARAM_OCEAN_DYN. Check!'
      call prc_abort
    endif
    log_nml(param_ocean_dyn)
    
    !-
    select case (trim(eqs_type))
    case ('FIXED_INIT_TEMP')
      this%eqs_type = ocean_dyn_eqs_fixed_init_temp
    case ('SLAB')
      this%eqs_type = ocean_dyn_eqs_slab
    case default
      log_info('OceanDyn_setup',*) 'Not appropriate EQS_TYPE. Check!'
      call prc_abort
    end select
 
    !- get mesh --------------------------------------------------
 
    call model_mesh%GetModelMesh( ptr_mesh )
    select type(model_mesh)
    class is (oceanmesh)
      ocean_mesh => model_mesh
    end select
    this%mesh3D => ocean_mesh%ptr_mesh
 
    !- Setup the temporal integrator
 
    call tm_parent_comp%Regist_process( 'OCEAN_DYN', time_dt, time_dt_unit, & ! (in)
      this%tm_process_id )                                                    ! (out)
 
    this%dtsec = tm_parent_comp%process_list(this%tm_process_id)%dtsec
 
    !- initialize the variables 
    call this%dyn_vars%Init( model_mesh )
 
    call setup_coriolis_parameter( this%dyn_vars, ocean_mesh )
 
    !- Initialize a module for 3D dynamical core 
    ! call this%dyncore_driver%Init( EQS_TYPE, &
    !   TINTEG_TYPE, dtsec,                    &
    !   MODALFILTER_FLAG,                      &
    !   atm_mesh )
 
    return
  end subroutine oceandyn_setup
 
 
!> Calculate tendencies associated with oceanic dynamics
!!
!! Because the tendencies with oceanic dynamical cores are treated in OceanDyn_update,
!! no calculation is performed in this subroutine.
!!
!! @param model_mesh Object to manage computational mesh of oceanic model 
!! @param prgvars_list Object to mange prognostic variables with oceanic dynamical core
!! @param trcvars_list Object to mange auxiliary variables 
!! @param forcing_list Object to mange forcing terms
!! @param is_update Flag to specify whether the tendencies are updated in this call
!!
!OCL SERIAL  
  subroutine oceandyn_calc_tendency( this, model_mesh, prgvars_list, trcvars_list, auxvars_list, forcing_list, is_update )
    implicit none
    
    class(oceandyn), intent(inout) :: this
    class(modelmeshbase), intent(in) :: model_mesh
    class(modelvarmanager), intent(inout) :: prgvars_list
    class(modelvarmanager), intent(inout) :: trcvars_list    
    class(modelvarmanager), intent(inout) :: auxvars_list
    class(modelvarmanager), intent(inout) :: forcing_list
    logical, intent(in) :: is_update
 
    !--------------------------------------------------
 
    if (.not. this%IsActivated()) return
    !LOG_INFO('OceanDyn_tendency',*)
    return  
  end subroutine oceandyn_calc_tendency
 
 
!> Update variables with a component of oceanic dynamics
!! The tendencies with oceanic dynamical cores are evaluated and the prognostic variables is updated.
!!
!! @param model_mesh Object to manage computational mesh of oceanic model 
!! @param prgvars_list Object to mange prognostic variables with oceanic dynamical core
!! @param trcvars_list Object to mange auxiliary variables 
!! @param forcing_list Object to mange forcing terms
!! @param is_update Flag to speicfy whether the tendencies are updated in this call
!!
!OCL SERIAL
  subroutine oceandyn_update( this, model_mesh, prgvars_list, trcvars_list, auxvars_list, forcing_list, is_update )
    use scale_meshfield_base, only: &
      meshfield2d, meshfield3d
    use mod_ocean_vars, only: &
      prgvar_therm_id,     &
      phytend_rhoh_id
    implicit none
 
    class(oceandyn), intent(inout) :: this
    class(modelmeshbase), intent(in) :: model_mesh
    class(modelvarmanager), intent(inout) :: prgvars_list
    class(modelvarmanager), intent(inout) :: trcvars_list    
    class(modelvarmanager), intent(inout) :: auxvars_list
    class(modelvarmanager), intent(inout) :: forcing_list
    logical, intent(in) :: is_update
 
    class(meshbase), pointer :: mesh
    class(meshbase3d), pointer :: mesh3d
 
    integer :: idom
    class(localmesh3d), pointer :: lmesh
    class(localmesh2d), pointer :: lmesh2d
 
    class(meshfield3d), pointer :: therm, rhoh
    !--------------------------------------------------
  
    if ( this%eqs_type == ocean_dyn_eqs_fixed_init_temp ) then
      log_info('OceanDyn_update',*) 'No update for fixed initial temperature'
      return
    end if
 
    call prof_rapstart( 'OCEAN_DYN_update', 1)
 
    call model_mesh%GetModelMesh( mesh )
    select type(mesh)
    class is (meshbase3d)
      mesh3d => mesh
    end select
 
 
    call prof_rapstart( 'OCEAN_DYN_core', 2)
    call prgvars_list%Get3D( prgvar_therm_id, therm )
    call forcing_list%Get3D( phytend_rhoh_id, rhoh )
    do idom=1, mesh3d%LOCAL_MESH_NUM
      lmesh => mesh3d%lcmesh_list(idom)
      lmesh2d => lmesh%lcmesh2D
      call ocn_dyn_update_lc( therm%local(idom)%val, &
        rhoh%local(idom)%val,                                          &
        this%dtsec, lmesh, lmesh%refElem3D, lmesh2d, lmesh2d%refElem2D )
 
    end do
 
    call prof_rapend( 'OCEAN_DYN_core', 2)
 
    !---------------------------
    call prof_rapend( 'OCEAN_DYN_update', 1)
 
    return  
  end subroutine oceandyn_update
 
!> Finalize an object to manage a component of oceanic dynamics
!!
!OCL SERIAL
  subroutine oceandyn_finalize( this )
    implicit none
    class(oceandyn), intent(inout) :: this
 
    integer :: n
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
    log_info('OceanDyn_finalize',*)
 
    ! call this%dyncore_driver%Final()
    call this%dyn_vars%Final()
 
    return  
  end subroutine oceandyn_finalize
 
  !--- private ---------------
 
  !> Setup Coriolis parameter
!OCL SERIAL
  subroutine setup_coriolis_parameter( this, ocn_mesh )
 
    use scale_coriolis_param, only: get_coriolis_parameter
    use scale_mesh_cubedom3d, only: meshcubedom3d
    implicit none
 
    class(oceandynvars), target, intent(inout) :: this
    class(oceanmesh), target, intent(in) :: ocn_mesh
 
    class(localmeshfieldbase), pointer :: coriolis
    class(localmesh3d), pointer :: lcmesh3d
    class(localmesh2d), pointer :: lcmesh2d
    integer :: n
 
    character(len=H_SHORT) :: coriolis_type !< Type of coriolis force: 'PLANE', 'SPHERE'
    real(rp) :: coriolis_f0   = 0.0_rp
    real(rp) :: coriolis_beta = 0.0_rp
    real(rp) :: coriolis_y0
 
    namelist /param_ocean_dyn_coriolis/ &
      coriolis_type,                         &                
      coriolis_f0, coriolis_beta, coriolis_y0
        
    class(meshbase3d), pointer :: mesh3d
    class(meshcubedom3d), pointer :: meshcube
    integer :: ierr
    !---------------------------------------------------------------
 
    mesh3d => ocn_mesh%ptr_mesh
 
    coriolis_type = 'NONE'
 
    select type(mesh3d)
    type is (meshcubedom3d)
      coriolis_y0 = 0.5_rp*(mesh3d%ymax_gl +  mesh3d%ymin_gl)
    end select
 
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_ocean_dyn_coriolis,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("OCEAN_DYN_setup_coriolis",*) 'Not found namelist. Default used.'
    else if( ierr > 0 ) then !--- fatal error
      log_error("OCEAN_DYN_setup_coriolis",*) 'Not appropriate names in namelist PARAM_OCEAN_DYN_CORIOLIS. Check!'
      call prc_abort
    end if
    log_nml(param_ocean_dyn_coriolis)
 
    do n = 1, mesh3d%LOCAL_MESH_NUM
      call oceandynauxvars_getlocalmeshfields( n, mesh3d, this%AUXVARS2D_manager, &
        coriolis, lcmesh3d )
      lcmesh2d => lcmesh3d%lcmesh2D
 
      call get_coriolis_parameter( &
        coriolis%val(:,lcmesh2d%NeS:lcmesh2d%NeE),                       & ! (out)
        coriolis_type, lcmesh2d%refElem2D%Np * lcmesh2d%Ne,              & ! (in)
        lcmesh2d%pos_en(:,:,2), coriolis_f0, coriolis_beta, coriolis_y0, & ! (in)
        lcmesh3d%lat2D                                                   ) ! (in)
      
      !$acc update device( coriolis%val )
    end do
 
    return
  end subroutine setup_coriolis_parameter
 
!-- private
!OCL SERIAL
  subroutine ocn_dyn_update_lc( THERM, &
    RHOH, dt, lmesh, elem, lmesh2D, elem2D )
    use scale_const, only: &
      dwatr => const_dwatr
    use scale_atmos_hydrometeor, only: &
       cv_water      
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(inout) :: therm(elem%np,lmesh%nea)
    real(rp), intent(in) :: rhoh(elem%np,lmesh%nea)
    real(rp), intent(in) :: dt
 
    integer :: ke
    real(rp) :: coef
    !----------------------------------------------------------------------
 
    coef = dt / ( cv_water * dwatr )
    !$omp parallel do
    do ke=lmesh%NeS, lmesh%NeE
      therm(:,ke) = therm(:,ke) + rhoh(:,ke) * coef
    end do
    return
  end subroutine ocn_dyn_update_lc
end module mod_ocean_dyn