scale_atmos_phy_cp_kf_jmapplib.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_cp_kf_jmapplib
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
 
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_cp_kf_jmapplib_setup
  public :: atmos_phy_cp_kf_jmapplib_finalize
  public :: atmos_phy_cp_kf_jmapplib_tendency
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
#ifdef JMAPPLIB
  character(len=6), private :: ATMOS_PHY_CP_KF_JMAPPLIB_type = "KF1701"
  real(RP), allocatable, private  :: ishall_counter(:,:)
  real(RP), allocatable, private  :: ideep_counter (:,:)
#endif
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_cp_kf_jmapplib_setup( &
       KA, KS, KE, IA, JA, &
       dx, dt )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       pre00 => const_pre00
#ifdef JMAPPLIB
    use conv_grid, only: &
       conv_grid_ini
    use conv_kf_const, only: &
       conv_kf_const_ini
    use conv_kf_parm, only: &
       conv_kf_parm_ini, &
       conv_kf_parm_get_stepcu
    use conv_kf_parm1701, only: &
       conv_kf_parm_ini1701, &
       conv_kf_parm_get_stepcu1701
    use conv_kf_lut, only: &
       conv_kf_lut_ini
#endif
    implicit none
 
    integer,  intent(in) :: ka, ks, ke
    integer,  intent(in) :: ia, ja
 
    real(rp), intent(in) :: dx
    real(dp), intent(in) :: dt
 
#ifdef JMAPPLIB
 
    integer  :: atmos_phy_cp_kf_jmapplib_radius_type   = 3
    real(rp) :: atmos_phy_cp_kf_jmapplib_dtlcl_fact    = 0.5_rp
    real(rp) :: atmos_phy_cp_kf_jmapplib_dpmin         = 5000.0_rp
    real(rp) :: atmos_phy_cp_kf_jmapplib_dlifetime     = 600.0_rp
    real(rp) :: atmos_phy_cp_kf_jmapplib_slifetime     = 600.0_rp
    real(rp) :: atmos_phy_cp_kf_jmapplib_detlq_qr_rate = 0.0_rp
    real(rp) :: atmos_phy_cp_kf_jmapplib_detic_qs_rate = 1.0_rp
 
    namelist / param_atmos_phy_cp_kf_jmapplib / &
         atmos_phy_cp_kf_jmapplib_type, &
         atmos_phy_cp_kf_jmapplib_radius_type, &
         atmos_phy_cp_kf_jmapplib_dtlcl_fact, &
         atmos_phy_cp_kf_jmapplib_dpmin, &
         atmos_phy_cp_kf_jmapplib_dlifetime, &
         atmos_phy_cp_kf_jmapplib_slifetime, &
         atmos_phy_cp_kf_jmapplib_detlq_qr_rate, &
         atmos_phy_cp_kf_jmapplib_detic_qs_rate
 
    integer :: kmax
 
    integer :: ierr
    integer :: k
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_CP_KF_JMAPPLIB_setup",*) 'Setup'
    log_info("ATMOS_PHY_CP_KF_JMAPPLIB_setup",*) 'KF scheme implemented in the JMA Physics Process Library'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_cp_kf_jmapplib,iostat=ierr)
    if( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_CP_KF_JMAPPLIB_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_CP_KF_JMAPPLIB. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_cp_kf_jmapplib)
 
 
    kmax = ke - ks + 1
    call conv_grid_ini( kmax )
    call conv_kf_const_ini( dx, pre00, dt )
 
    select case ( atmos_phy_cp_kf_jmapplib_type )
    case ( "KF" )
 
       call conv_kf_parm_ini( &
            dt,                                                           &
            kfrad_var_in        = atmos_phy_cp_kf_jmapplib_radius_type,   &
            cudt_in             = -1.0_rp,                                &
            dtlcl_fct_in        = atmos_phy_cp_kf_jmapplib_dtlcl_fact,    &
            dpmin_in            = atmos_phy_cp_kf_jmapplib_dpmin,         &
            cu_lifetime_min_in  = atmos_phy_cp_kf_jmapplib_dlifetime,     &
            shallow_lifetime_in = atmos_phy_cp_kf_jmapplib_slifetime,     &
            detlq_qr_r_in       = atmos_phy_cp_kf_jmapplib_detlq_qr_rate, &
            detic_qs_r_in       = atmos_phy_cp_kf_jmapplib_detic_qs_rate  )
 
!       call conv_kf_parm_get_stepcu(stepcu) ! (out)
 
    case ( "KF1701" )
 
       call conv_kf_parm_ini1701( &
            dt,                                                           &
            kfrad_var_in        = atmos_phy_cp_kf_jmapplib_radius_type,   &
            kfrad_in            = 1000.0_rp,                              &
            cudt_in             = -1.0_rp,                                &
            dpmin_in            = atmos_phy_cp_kf_jmapplib_dpmin,         &
            cu_lifetime_min_in  = atmos_phy_cp_kf_jmapplib_dlifetime,     &
            shallow_lifetime_in = atmos_phy_cp_kf_jmapplib_slifetime,     &
            detlq_qr_r_in       = atmos_phy_cp_kf_jmapplib_detlq_qr_rate, &
            detic_qs_r_in       = atmos_phy_cp_kf_jmapplib_detic_qs_rate, &
            tkemax_in           = 1.0_rp                                  )
            
!       call conv_kf_parm_get_stepcu1701(stepcu) ! (out)
 
    case default
       log_error("ATMOS_PHY_CP_KF_JMAPPLIB_setup",*) 'ATMOS_PHY_CP_KF_TYPE must be either "KF" or "KF1701"'
       call prc_abort
    end select
 
    call conv_kf_lut_ini
 
    allocate( ishall_counter(ia,ja), ideep_counter(ia,ja) )
    ishall_counter(:,:) = 0.0_rp
    ideep_counter(:,:) = 0.0_rp
 
#else
 
    log_error("ATMOS_PHY_CP_KF_JMAPPLIB_setup",*) 'To use "KF-JMAPPLIB", compile SCALE with "SCALE_ENABLE_JMAPPLIB=T" option.'
    call prc_abort
 
#endif
 
    return
  end subroutine atmos_phy_cp_kf_jmapplib_setup
 
  !-----------------------------------------------------------------------------
  !! Finalize
  
  subroutine atmos_phy_cp_kf_jmapplib_finalize
    implicit none
 
#ifdef JMAPPLIB
    deallocate( ishall_counter, ideep_counter )
#endif
 
    return
  end subroutine atmos_phy_cp_kf_jmapplib_finalize
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_cp_kf_jmapplib_tendency( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       DENS, U, V, W, TEMP, POTT,       &
       PRES, EXNER,                     &
       QDRY, QV, QC, QI,                &
       us, PBLH, SFLX_BUOY,             &
       CZ, FZ,                          &
       DENS_t, RHOT_t,                  &
       RHOQV_t, RHOQ_t,                 &
       w0avg, nca,                      &
       SFLX_rain, SFLX_snow, SFLX_prec, &
       cloudtop, cloudbase              )
    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc, &
       i_hr, &
       i_hi, &
       i_hs, &
       i_hg
#ifdef JMAPPLIB
    use conv_kf_main, only: &
       conv_kf_run
    use conv_kf_main1701, only: &
       conv_kf_run1701
    use pbl_diag, only: &
       pbl_diag_buoyancy_excess
#endif
    implicit none
 
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp), intent(in) :: dens (ka,ia,ja)
    real(rp), intent(in) :: u    (ka,ia,ja)
    real(rp), intent(in) :: v    (ka,ia,ja)
    real(rp), intent(in) :: w    (ka,ia,ja)
    real(rp), intent(in) :: temp (ka,ia,ja)
    real(rp), intent(in) :: pott (ka,ia,ja)
    real(rp), intent(in) :: pres (ka,ia,ja)
    real(rp), intent(in) :: exner(ka,ia,ja)
    real(rp), intent(in) :: qdry (ka,ia,ja)
    real(rp), intent(in) :: qv   (ka,ia,ja)
    real(rp), intent(in) :: qc   (ka,ia,ja)
    real(rp), intent(in) :: qi   (ka,ia,ja)
    real(rp), intent(in) :: us       (ia,ja)
    real(rp), intent(in) :: pblh     (ia,ja)
    real(rp), intent(in) :: sflx_buoy(ia,ja)
    real(rp), intent(in) :: cz   (ka,ia,ja)
    real(rp), intent(in) :: fz (0:ka,ia,ja)
 
    real(rp), intent(out)   :: dens_t (ka,ia,ja)
    real(rp), intent(inout) :: rhot_t (ka,ia,ja)
    real(rp), intent(inout) :: rhoqv_t(ka,ia,ja)
    real(rp), intent(inout) :: rhoq_t (ka,ia,ja,n_hyd)
    real(rp), intent(inout) :: w0avg  (ka,ia,ja)
    real(rp), intent(inout) :: nca      (ia,ja)
    real(rp), intent(inout) :: sflx_rain(ia,ja)
    real(rp), intent(inout) :: sflx_snow(ia,ja)
    real(rp), intent(inout) :: sflx_prec(ia,ja)
    real(rp), intent(inout) :: cloudtop (ia,ja)
    real(rp), intent(inout) :: cloudbase(ia,ja)
 
#ifdef JMAPPLIB
 
    integer, parameter :: iconv_run = 1
 
    real(rp) :: tend_pt(ks:ke)
    real(rp) :: tend_qv(ks:ke)
    real(rp) :: tend_qc(ks:ke)
    real(rp) :: tend_qr(ks:ke)
    real(rp) :: tend_qi(ks:ke)
    real(rp) :: tend_qs(ks:ke)
    real(rp) :: tend_qg(ks:ke)
    real(rp) :: dens_lc(ks:ke)
    real(rp) :: qv_lc  (ks:ke)
    real(rp) :: qc_lc
    real(rp) :: qi_lc
    real(rp) :: z_f    (ks:ke)
    real(rp) :: dz_f   (ks:ke)
 
    real(rp) :: lcl_temp(ia,ja)
    real(rp) :: abe_out
    real(rp) :: tv_ex
    real(rp) :: ptv_ex
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    log_progress(*) "atmosphere / physics / cumulus / KF-JMAPPLIB"
 
    !$omp parallel do &
    !$omp private(tend_pt,tend_qv,tend_qc,tend_qr,tend_qi,tend_qs,tend_qg, &
    !$omp         dens_lc,qv_lc,qc_lc,qi_lc,z_f,dz_f, &
    !$omp         ptv_ex,tv_ex,abe_out)
    do j = js, je
    do i = is, ie
 
       do k = ks, ke
          dens_lc(k) = dens(k,i,j) * qdry(k,i,j)
          qv_lc(k) = qv(k,i,j) / qdry(k,i,j)
          tend_pt(k) = rhot_t(k,i,j) / dens(k,i,j)
          tend_qv(k) = rhoqv_t(k,i,j) / dens_lc(k)
          tend_qc(k) = rhoq_t(k,i,j,i_hc) / dens_lc(k)
          tend_qr(k) = rhoq_t(k,i,j,i_hr) / dens_lc(k)
          tend_qi(k) = rhoq_t(k,i,j,i_hi) / dens_lc(k)
          tend_qs(k) = rhoq_t(k,i,j,i_hs) / dens_lc(k)
          tend_qg(k) = rhoq_t(k,i,j,i_hg) / dens_lc(k)
 
          z_f(k) = cz(k,i,j) - fz(ks-1,i,j)
          dz_f(k) = fz(k,i,j) - fz(k-1,i,j)
       end do
 
       select case ( atmos_phy_cp_kf_jmapplib_type )
       case ( "KF" )
 
          call conv_kf_run( &
               iconv_run,                                                           & ! (in)
               u(ks:ke,i,j), v(ks:ke,i,j), w(ks:ke,i,j), temp(ks:ke,i,j), qv_lc(:), & ! (in)
               pres(ks:ke,i,j), exner(ks:ke,i,j), dens_lc(:), z_f(:), dz_f(:),      & ! (in)
               nca(i,j),                                                            & ! (inout)
               sflx_prec(i,j), sflx_rain(i,j), sflx_snow(i,j),                      & ! (inout)
               cloudtop(i,j), cloudbase(i,j), abe_out,                              & ! (inout)
               ishall_counter(i,j), ideep_counter(i,j),                             & ! (inout)
               lcl_temp(i,j),                                                       & ! (inout)
               w0avg(ks:ke,i,j),                                                    & ! (inout)
               tend_pt(:), tend_qv(:), tend_qc(:),                                  & ! (inout)
               tend_qi(:), tend_qr(:), tend_qs(:)                                   ) ! (inout)
 
          tend_qg(ks:ke) = 0.0_rp
 
       case ( "KF1701" )
 
          qc_lc = qc(ks,i,j) / qdry(ks,i,j)
          qi_lc = qi(ks,i,j) / qdry(ks,i,j)
          call pbl_diag_buoyancy_excess( &
               pott(ks,i,j), qv_lc(ks),qc_lc, qi_lc,              & ! (in)
               pblh(i,j), us(i,j), sflx_buoy(i,j), exner(ks,i,j), & ! (in)
               ptv_ex, tv_ex                                      ) ! (out)
 
          call conv_kf_run1701( &
               iconv_run,                                                           & ! (in)
               u(ks:ke,i,j), v(ks:ke,i,j), w(ks:ke,i,j), temp(ks:ke,i,j), qv_lc(:), & ! (in)
               pres(ks:ke,i,j), exner(ks:ke,i,j), dens_lc(:), z_f(:), dz_f(:),      & ! (in)
               pblh(i,j), tv_ex,                                                    & ! (in)
               nca(i,j),                                                            & ! (inout)
               sflx_prec(i,j), sflx_rain(i,j), sflx_snow(i,j),                      & ! (inout)
               cloudtop(i,j), cloudbase(i,j), abe_out,                              & ! (inout)
               ishall_counter(i,j), ideep_counter(i,j),                             & ! (inout)
               lcl_temp(i,j),                                                       & ! (inout)
               w0avg(ks:ke,i,j),                                                    & ! (inout)
               tend_pt(:), tend_qv(:), tend_qc(:),                                  & ! (inout)
               tend_qi(:), tend_qr(:),                                              & ! (inout)
               tend_qs(:), tend_qg(:)                                               ) ! (inout)
 
       end select
 
       if ( nca(i,j) < -1 ) nca(i,j) = -1
 
       do k = ks, ke
          rhot_t(k,i,j)      = tend_pt(k) * dens(k,i,j)
          rhoqv_t(k,i,j)      = tend_qv(k) * dens_lc(k)
          rhoq_t(k,i,j,i_hc) = tend_qc(k) * dens_lc(k)
          rhoq_t(k,i,j,i_hr) = tend_qr(k) * dens_lc(k)
          rhoq_t(k,i,j,i_hi) = tend_qi(k) * dens_lc(k)
          rhoq_t(k,i,j,i_hs) = tend_qs(k) * dens_lc(k)
          rhoq_t(k,i,j,i_hg) = tend_qg(k) * dens_lc(k)
          dens_t(k,i,j)      = rhoqv_t(k,i,j) &
                              + rhoq_t(k,i,j,i_hc) + rhoq_t(k,i,j,i_hr) &
                              + rhoq_t(k,i,j,i_hi) + rhoq_t(k,i,j,i_hs) + rhoq_t(k,i,j,i_hg)
       end do
 
    end do
    end do
 
#endif
 
    return
  end subroutine atmos_phy_cp_kf_jmapplib_tendency
 
end module scale_atmos_phy_cp_kf_jmapplib