mod_atmos_phy_rd.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Physics / Radiation
!!
!! @par Description
!!          Module for radiation component in atmospheric model
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_atmos_phy_rd
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !  
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prof
  use scale_const, only: &
    undef8 => const_undef8
 
  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, &
    elementbase1d, elementbase2d, elementbase3d
 
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield2d, meshfield3d
  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_atmos_phy_rd_vars, only: atmosphyrdvars
  use mod_atmos_vars_container, only: &
    atmosvarscontainer
 
  use scale_atm_phy_rd_dgm_simple, only: atmphyradsimple
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  !> Derived type to manage a component of radiation in atmospheric model
  !!
  type, extends(modelcomponentproc), public :: atmosphyrd
    integer :: rd_typeid          !< Type id of radiation scheme
    integer :: rd_solarins_typeid !< Type id of solar insolation scheme
 
    type(atmosphyrdvars) :: vars               !< Object to manage variables with radiation
    type(meshfield2d), pointer :: sfc_temp_ptr !< Pointer to an object with surface temperature field
    type(meshfield2d), pointer :: sfc_alb_ptr  !< Pointer to an object with surface albedo field
 
    integer :: atm_var_container_typeid        !< Type ID of variable container for radiation
 
    !-
    type(atmphyradsimple) :: simple_rad !< Object to manage a simplified-radiation scheme
 
    ! type(MeshField3D) :: radflux_LW_dn_save
    ! type(MeshField3D) :: radflux_LW_up_save
  contains
    procedure, public :: setup => atmosphyrd_setup
    procedure, public :: calc_tendency => atmosphyrd_calc_tendency
    procedure, public :: update => atmosphyrd_update
    procedure, public :: finalize => atmosphyrd_finalize
    procedure, public :: setsfcvars => atmosphyrd_set_sfcvars
    procedure, private :: calc_tendency_core => atmosphyrd_calc_tendency_core
  end type atmosphyrd
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, parameter :: rd_typeid_simple  = 1           !< Type ID of a simple radiation scheme
 
  integer, parameter :: rd_insolation_typeid_none   = 0 !< Type ID of no solar insolation
  integer, parameter :: rd_insolation_typeid_simple = 1 !< Type ID of a simple solar insolation scheme
contains
 
  !> Setup a component of radiation in atmospheric model
  !!
  !! @param model_mesh Object to manage computational mesh of atmospheric model 
  !! @param tm_parent_comp Object to mange a temporal scheme in a parent component
  !!
  subroutine atmosphyrd_setup( this, model_mesh, tm_parent_comp )
    use mod_atmos_mesh, only: atmosmesh
    use scale_time_manager, only: time_manager_component
    use mod_atmos_vars, only: atm_vars_container_primary_id
    use scale_atm_phy_rd_solarins_simple, only: &
      atm_phy_rd_solarins_simple_setup
    implicit none
    class(atmosphyrd), intent(inout) :: this
    class(modelmeshbase), target, intent(in) :: model_mesh
    class(time_manager_component), intent(inout) :: tm_parent_comp
 
    real(dp) :: time_dt                             = undef8 !< Timestep for radiation
    character(len=H_SHORT) :: time_dt_unit          = 'SEC'  !< Unit of timestep
 
    character(len=H_MID) :: rd_type          = 'NONE'        !< Type of a radiation scheme [NONE, GRAYRAD]
    character(len=H_MID) :: rd_solarins_type = 'SIMPLE'      !< Type of solar insolation scheme [NONE, SIMPLE, REAL]
    integer :: atm_var_container_typeid
 
    namelist /param_atmos_phy_rd/ &
      time_dt,             &
      time_dt_unit,        &
      rd_type,             &
      rd_solarins_type,    &
      atm_var_container_typeid
 
    class(atmosmesh), pointer     :: atm_mesh
    class(meshbase), pointer      :: ptr_mesh
    class(localmesh3d), pointer :: lcmesh3d
    class(elementbase3d), pointer :: elem3d
 
    class(meshbase3d), pointer :: mesh3d
 
    integer :: ierr
    !-----------------------------------------------------
 
    if (.not. this%IsActivated()) return
 
    log_newline
    log_info("ATMOS_PHY_RD_setup",*) 'Setup'
 
    atm_var_container_typeid = atm_vars_container_primary_id
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_rd,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_PHY_RD_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_PHY_RD_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_RD. Check!'
      call prc_abort
    endif
    log_nml(param_atmos_phy_rd)
 
    this%atm_var_container_typeid = atm_var_container_typeid
 
    !--- Register this component in the time manager
    
    call tm_parent_comp%Regist_process( 'ATMOS_PHY_RD', time_dt, time_dt_unit, & ! (in)
      this%tm_process_id )                                                       ! (out) 
 
    !--- Set the type of radiation scheme
    
    select case ( rd_type )
    case ("SIMPLE")
      this%RD_TYPEID = rd_typeid_simple
      call this%simple_rad%Init()
    case default
      log_error("ATMOS_PHY_RD_setup",*) 'Not appropriate RD_TYPE. Check!'
      call prc_abort
    end select
 
    select case ( rd_solarins_type )
    case ("SIMPLE")
      this%RD_SOLARINS_TYPEID = rd_insolation_typeid_simple
      call atm_phy_rd_solarins_simple_setup()
    case ("NONE")
      this%RD_SOLARINS_TYPEID = rd_insolation_typeid_none
    case default
      log_error("ATMOS_PHY_RD_setup",*) 'Not appropriate RD_SOLARINS_TYPE. Check!'
      call prc_abort
    end select
 
 
    !- Initialize the variables 
    call this%vars%Init( model_mesh )
 
    select type(model_mesh)
    class is (atmosmesh)
      atm_mesh => model_mesh
    end select
    mesh3d => atm_mesh%ptr_mesh 
    ! call this%radflux_LW_dn_save%Init( "RADFLUX_LW_dn", "W/m2", mesh3D )
    ! call this%radflux_LW_up_save%Init( "RADFLUX_LW_up", "W/m2", mesh3D )
 
    return
  end subroutine atmosphyrd_setup
 
!OCL SERIAL
  subroutine atmosphyrd_set_sfcvars( this, sfc_temp, sfc_alb )
    implicit none
    class(atmosphyrd), intent(inout) :: this
    type(meshfield2d), target, intent(in) :: sfc_temp
    type(meshfield2d), target, intent(in) :: sfc_alb
    !-----------------------------------------------------
    this%SFC_TEMP_ptr => sfc_temp
    this%SFC_ALB_ptr => sfc_alb
    return
  end subroutine atmosphyrd_set_sfcvars
 
  !> Calculate tendencies associated with radiation in atmospheric model
  !!
  !!
  !! @param model_mesh Object to manage computational mesh of atmospheric model 
  !! @param prgvars_list Object to manage prognostic variables with atmospheric dynamical core
  !! @param trcvars_list Object to manage auxiliary variables 
  !! @param forcing_list Object to manage forcing terms
  !! @param is_update Flag to speicfy whether the tendencies are updated in this call
  !!
!OCL SERIAL
  subroutine atmosphyrd_calc_tendency( &
    this, model_mesh, prgvars_list, trcvars_list, &
    auxvars_list, forcing_list, is_update         )
    use scale_atm_phy_rd_dgm_common, only: atm_phy_rd_dgm_calc_heating
    use scale_atm_phy_rd_solarins_simple, only: atm_phy_rd_solarins_simple_get
    use mod_atmos_vars, only: &
      atmosvars_getlocalmeshprgvars,     &
      atmosvars_getlocalmeshphyauxvars,  &
      atmosvars_getlocalmeshqtrc_qv,     &
      atmosvars_getlocalmeshphytends
    use mod_atmos_phy_rd_vars, only: &
      rd_rhoh_id => atmos_phy_rd_rhoh_id, &
      solins_id => atmos_phy_rd_aux2d_solins_id,         &
      cos_sza_id => atmos_phy_rd_aux2d_cossza_id,        &
      sflx_sw_up_id => atmos_phy_rd_aux2d_sflx_sw_up_id, &
      sflx_sw_dn_id => atmos_phy_rd_aux2d_sflx_sw_dn_id, &
      sflx_lw_up_id => atmos_phy_rd_aux2d_sflx_lw_up_id, &
      sflx_lw_dn_id => atmos_phy_rd_aux2d_sflx_lw_dn_id, &
      tomflx_lw_up_id => atmos_phy_rd_aux2d_tomflx_lw_up_id
    implicit none
    class(atmosphyrd), 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(localmesh3d), pointer :: lcmesh
    class(localmesh2d), pointer :: lcmesh2d
 
    class(localmeshfieldbase), pointer :: ddens, momx, momy, momz, drhot
    class(localmeshfieldbase), pointer :: rtot, cvtot, cptot
    class(localmeshfieldbase), pointer :: dens_hyd, pres_hyd
    class(localmeshfieldbase), pointer :: pres, pt
    class(localmeshfieldbase), pointer :: qv, qv_tp
    class(localmeshfieldbase), pointer :: rd_rhoh
    class(localmeshfieldbase), pointer :: dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp, rhoh_p
 
    integer :: n
    integer :: ke
    !------------------------------------------------------------------------
 
    if (.not. this%IsActivated()) return
 
    log_progress(*) 'atmosphere / physics / radiation'
 
    call model_mesh%GetModelMesh( mesh )
 
    if (is_update) then
    
      do n=1, mesh%LOCAL_MESH_NUM
        call prof_rapstart('ATM_PHY_RD_get_localmesh_ptr', 2)    
 
        !- Get pointers to the fields in the variable container for cloud microphysics
        call atmosvars_getlocalmeshprgvars( n,  &
          mesh, prgvars_list, auxvars_list,       &
          ddens, momx, momy, momz, drhot,         &
          dens_hyd, pres_hyd, rtot, cvtot, cptot, &
          lcmesh                                  )
 
        call atmosvars_getlocalmeshphyauxvars( n,  &
          mesh, auxvars_list,                      &
          pres, pt )
 
        call atmosvars_getlocalmeshqtrc_qv( n, &
          mesh, trcvars_list, forcing_list,    &
          qv, qv_tp )
        call prof_rapend('ATM_PHY_RD_get_localmesh_ptr', 2)
 
        !- Calculate tendencies associated with radiation
 
        call prof_rapstart('ATM_PHY_RD_cal_tend', 2)
        
        lcmesh2d => lcmesh%lcmesh2D
 
        select case ( this%RD_SOLARINS_TYPEID )
        case ( rd_insolation_typeid_simple )
          call atm_phy_rd_solarins_simple_get( &
            this%vars%auxvars2D(solins_id)%local(n)%val(:,lcmesh2d%NeS:lcmesh2d%NeE),  & ! (out)
            this%vars%auxvars2D(cos_sza_id)%local(n)%val(:,lcmesh2d%NeS:lcmesh2d%NeE), & ! (out)
            lcmesh2d%lat, lcmesh2d%Ne * lcmesh2d%refElem2D%Np )
        end select
 
        call this%vars%tends(rd_rhoh_id)%GetLocalMeshField( n, rd_rhoh )
        call this%calc_tendency_core( rd_rhoh%val,                & ! (out)
          this%vars%auxvars2D(sflx_sw_up_id)%local(n)%val,        & ! (out)
          this%vars%auxvars2D(sflx_sw_dn_id)%local(n)%val,        & ! (out)
          this%vars%auxvars2D(sflx_lw_up_id)%local(n)%val,        & ! (out)
          this%vars%auxvars2D(sflx_lw_dn_id)%local(n)%val,        & ! (out)
          this%vars%auxvars2D(tomflx_lw_up_id)%local(n)%val,      & ! (out)
          this%vars%auxvars2D(solins_id)%local(n)%val,            & ! (in)
          this%vars%auxvars2D(cos_sza_id)%local(n)%val,           & ! (in)
          ddens%val, pres%val, qv%val,                            & ! (in)
          this%SFC_TEMP_ptr%local(n)%val,                         & ! (in)
          this%SFC_ALB_ptr%local(n)%val,                          & ! (in)
          dens_hyd%val, rtot%val, cvtot%val,                      & ! (in)
          lcmesh, lcmesh%refElem3D, lcmesh2d, lcmesh2d%refElem2D, & ! (in)
          model_mesh%element3D_operation )                          ! (in)
        call prof_rapend('ATM_PHY_RD_cal_tend', 2)
      end do
    
    end if
 
    !- Add tendencies calculated in this component to the total tendencies
 
    do n=1, mesh%LOCAL_MESH_NUM
      call atmosvars_getlocalmeshphytends( n,        &
        mesh, forcing_list,                          &
        dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp, &
        rhoh_p, lcmesh3d=lcmesh )
 
      call this%vars%tends(rd_rhoh_id)%GetLocalMeshField( n, rd_rhoh )
      !$omp parallel private(ke)
      !$omp do
      do ke=lcmesh%NeS, lcmesh%NeE
        rhoh_p %val(:,ke) = rhoh_p %val(:,ke) + rd_rhoh%val(:,ke)
      end do
      !$omp end parallel
    end do
 
    return
  end subroutine atmosphyrd_calc_tendency
 
  !> Update variables in a component of radiation in atmospheric model
  !!
  !! @param model_mesh Object to manage computational mesh of atmospheric model 
  !! @param prgvars_list Object to manage prognostic variables with atmospheric dynamical core
  !! @param trcvars_list Object to manage auxiliary variables 
  !! @param forcing_list Object to manage forcing terms
  !! @param is_update Flag to speicfy whether the tendencies are updated in this call
  !!
!OCL SERIAL  
  subroutine atmosphyrd_update( this, model_mesh, &
    prgvars_list, trcvars_list,                   &
    auxvars_list, forcing_list, is_update         )  
    
    implicit none
    class(atmosphyrd), 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
    !--------------------------------------------------
    return
  end subroutine atmosphyrd_update
 
  !> Finalize a component of radiation in atmospheric model
  !!
!OCL SERIAL  
  subroutine atmosphyrd_finalize( this )
    implicit none
    class(atmosphyrd), intent(inout) :: this
 
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
 
    select case ( this%RD_TYPEID )
    case( rd_typeid_simple )
      call this%simple_rad%Final()
    end select
 
    call this%vars%Final()
    return
  end subroutine atmosphyrd_finalize
 
!- private ------------------------------------------------
!OCL SERIAL
  subroutine atmosphyrd_calc_tendency_core( this, &
    RHOH,                                           & ! (out)
    sflx_sw_up, sflx_sw_dn, sflx_lw_up, sflx_lw_dn, & ! (out)
    tomflx_lw_up,                                   & ! (out)
    solins, cos_sza,                                & ! (in)
    ddens, pres, qv, sfc_temp, sfc_alb,             & ! (in)
    dens_hyd, rtot, cvtot,                          & ! (in)
    lcmesh, elem3d, lcmesh2d, elem2d,               & ! (in)
    elem3d_operation )                                ! (in)
    use scale_atmos_phy_rd_common, only: &
      i_up, i_dn, i_lw, i_sw
    use scale_element_operation_base, only: elementoperationbase3d
    use scale_atm_phy_rd_dgm_common, only: atm_phy_rd_dgm_calc_heating
 
    use scale_meshfield_base, only: meshfield3d
    use scale_file_history_meshfield, only: file_history_meshfield_in
    implicit none
    class(atmosphyrd), intent(inout) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    class(localmesh2d), intent(in) :: lcmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(out) :: rhoh(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: sflx_sw_up(elem2d%np,lcmesh2d%nea)
    real(rp), intent(out) :: sflx_sw_dn(elem2d%np,lcmesh2d%nea)
    real(rp), intent(out) :: sflx_lw_up(elem2d%np,lcmesh2d%nea)
    real(rp), intent(out) :: sflx_lw_dn(elem2d%np,lcmesh2d%nea)
    real(rp), intent(out) :: tomflx_lw_up(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: solins(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: cos_sza(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: ddens(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: pres(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: qv(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: sfc_temp(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: sfc_alb(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: dens_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rtot(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: cvtot(elem3d%np,lcmesh%nea)
    class(elementoperationbase3d), intent(in) :: elem3d_operation
 
    real(rp) :: flux_rad(elem3d%np,lcmesh%ne,2,2,2)
    real(rp) :: flux_rad_top(elem3d%nnode_h1d**2,lcmesh%ne2d,2,2,2)
    real(rp) :: sflux_rad_up(elem3d%nnode_h1d**2,lcmesh%ne2d,2,2)
    real(rp) :: sflux_rad_dn(elem3d%nnode_h1d**2,lcmesh%ne2d,2,2)
    real(rp) :: temp_(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: dens_(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: pres_(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: qv_(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    
    ! real(RP) :: flux_up(elem3D%Nnode_v,lcmesh%NeZ,elem3D%Nnode_h1D**2,lcmesh%Ne2D)
    ! real(RP) :: flux_dn(elem3D%Nnode_v,lcmesh%NeZ,elem3D%Nnode_h1D**2,lcmesh%Ne2D)
    ! real(RP) :: flux_net(elem3D%Nnode_v,lcmesh%NeZ,elem3D%Nnode_h1D**2,lcmesh%Ne2D)
    ! real(RP) :: flux_net_sfc(elem3D%Nnode_h1D**2,lcmesh%Ne2D)
    ! real(RP) :: flux_net_toa(elem3D%Nnode_h1D**2,lcmesh%Ne2D)
    ! real(RP) :: flux_net_tom(elem3D%Nnode_h1D**2,lcmesh%Ne2D)
 
    integer :: ke, ke_h, ke_z
    integer :: p, ph, pz
 
    real(rp) :: dens_tmp
    !--------------------------------------------------
 
    !$omp parallel do private(ke,p, dens_tmp) collapse(2)
    do ke_z=1, lcmesh%NeZ
    do ke_h=1, lcmesh%Ne2D
      ke = ke_h + (ke_z-1)*lcmesh%Ne2D
      do pz=1, elem3d%Nnode_v
      do ph=1, elem3d%Nnode_h1D**2
        p = ph + (pz-1)*elem3d%Nnode_h1D**2
        dens_tmp = dens_hyd(p,ke) + ddens(p,ke)
        
        pres_(pz,ke_z,ph,ke_h) = pres(p,ke)
        temp_(pz,ke_z,ph,ke_h) = pres(p,ke) / ( rtot(p,ke) * dens_tmp )
        dens_(pz,ke_z,ph,ke_h) = dens_tmp
        qv_(pz,ke_z,ph,ke_h)   = qv(p,ke)
      end do
      end do
    end do
    end do
 
    select case( this%RD_TYPEID )
    case ( rd_typeid_simple )
      call this%simple_rad%calculate_rad_flux( flux_rad(:,:,:,:,1),            & ! (out)
        flux_rad_top(:,:,:,:,1), sflux_rad_up(:,:,:,1), sflux_rad_dn(:,:,:,1), & ! (out)
        solins, pres_, temp_, dens_, qv_, sfc_temp, sfc_alb,                   & ! (in)
        lcmesh, elem3d, lcmesh2d, elem2d )                                       ! (in)
      
      !$omp parallel do
      do ke=lcmesh2d%NeS, lcmesh2d%NeE
        sflx_sw_dn(:,ke) = sflux_rad_dn(:,ke,i_sw,1)
        sflx_lw_dn(:,ke) = sflux_rad_dn(:,ke,i_lw,1)
        sflx_sw_up(:,ke) = sflux_rad_up(:,ke,i_sw,1)
        sflx_lw_up(:,ke) = sflux_rad_up(:,ke,i_lw,1)
 
        tomflx_lw_up(:,ke) = flux_rad_top(:,ke,i_lw,1,1)
      end do
 
!     !$omp parallel do
!     do ke=lcmesh%NeS, lcmesh%NeE
!       this%radflux_LW_dn_save%local(1)%val(:,ke) = flux_rad(:,ke,I_LW,I_dn,1)
!       this%radflux_LW_up_save%local(1)%val(:,ke) = flux_rad(:,ke,I_LW,I_up,1)
!     end do
    end select
 
    call atm_phy_rd_dgm_calc_heating( rhoh, &
      flux_rad(:,:,:,:,1), ddens, dens_hyd, cvtot, &
      lcmesh, elem3d, elem2d, &
      elem3d_operation )
 
!      call FILE_HISTORY_meshfield_in( this%radflux_LW_dn_save, 'RADFLUX_LW_dn' )
!      call FILE_HISTORY_meshfield_in( this%radflux_LW_up_save, 'RADFLUX_LW_up' )
      return
  end subroutine atmosphyrd_calc_tendency_core
end module mod_atmos_phy_rd