scale_atmos_sub_thermodyn.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "macro_thermodyn.h"
#include "inc_openmp.h"
module scale_atmos_thermodyn
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer

  use scale_const, only: &
     eps   => const_eps,   &
     rdry  => const_rdry,  &
     cpdry => const_cpdry, &
     cpvap => const_cpvap, &
     cvdry => const_cvdry, &
     cl    => const_cl,    &
     ci    => const_ci,    &
     rvap  => const_rvap,  &
     lhv0  => const_lhv0,  &
     lhs0  => const_lhs0,  &
     lhf0  => const_lhf0,  &
     psat0 => const_psat0, &
     pre00 => const_pre00, &
     tem00 => const_tem00
  !-----------------------------------------------------------------------------
  implicit none
  private

  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_thermodyn_setup

  public :: atmos_thermodyn_qd
  public :: atmos_thermodyn_cv
  public :: atmos_thermodyn_cp
  public :: atmos_thermodyn_pott
  public :: atmos_thermodyn_rhoe
  public :: atmos_thermodyn_rhot
  public :: atmos_thermodyn_temp_pres
  public :: atmos_thermodyn_temp_pres_e
  public :: atmos_thermodyn_templhv
  public :: atmos_thermodyn_templhs
  public :: atmos_thermodyn_templhf
  public :: atmos_thermodyn_entr

  interface atmos_thermodyn_qd
     module procedure atmos_thermodyn_qd_0d
     module procedure atmos_thermodyn_qd_3d
  end interface atmos_thermodyn_qd

  interface atmos_thermodyn_cp
     module procedure atmos_thermodyn_cp_0d
     module procedure atmos_thermodyn_cp_3d
  end interface atmos_thermodyn_cp

  interface atmos_thermodyn_cv
     module procedure atmos_thermodyn_cv_0d
     module procedure atmos_thermodyn_cv_3d
  end interface atmos_thermodyn_cv

  interface atmos_thermodyn_rhoe
     module procedure atmos_thermodyn_rhoe_0d
     module procedure atmos_thermodyn_rhoe_3d
  end interface atmos_thermodyn_rhoe

  interface atmos_thermodyn_rhot
     module procedure atmos_thermodyn_rhot_0d
     module procedure atmos_thermodyn_rhot_3d
  end interface atmos_thermodyn_rhot

  interface atmos_thermodyn_pott
     module procedure atmos_thermodyn_pott_0d
     module procedure atmos_thermodyn_pott_3d
  end interface atmos_thermodyn_pott

  interface atmos_thermodyn_temp_pres
     module procedure atmos_thermodyn_temp_pres_0d
     module procedure atmos_thermodyn_temp_pres_3d
  end interface atmos_thermodyn_temp_pres

  interface atmos_thermodyn_temp_pres_e
     module procedure atmos_thermodyn_temp_pres_e_0d
     module procedure atmos_thermodyn_temp_pres_e_3d
  end interface atmos_thermodyn_temp_pres_e

  interface atmos_thermodyn_templhv
     module procedure atmos_thermodyn_templhv_0d
     module procedure atmos_thermodyn_templhv_2d
     module procedure atmos_thermodyn_templhv_3d
  end interface atmos_thermodyn_templhv

  interface atmos_thermodyn_templhs
     module procedure atmos_thermodyn_templhs_0d
     module procedure atmos_thermodyn_templhs_2d
     module procedure atmos_thermodyn_templhs_3d
  end interface atmos_thermodyn_templhs

  interface atmos_thermodyn_templhf
     module procedure atmos_thermodyn_templhf_0d
     module procedure atmos_thermodyn_templhf_2d
     module procedure atmos_thermodyn_templhf_3d
  end interface atmos_thermodyn_templhf

  interface atmos_thermodyn_entr
     module procedure atmos_thermodyn_entr_0d
     module procedure atmos_thermodyn_entr_2d
     module procedure atmos_thermodyn_entr_3d
  end interface atmos_thermodyn_entr

  public :: atmos_thermodyn_tempre
  public :: atmos_thermodyn_tempre2

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public              :: thermodyn_emask = 0.0_rp 

  real(RP), public, allocatable :: aq_cp(:)
  real(RP), public, allocatable :: aq_cv(:)

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_setup
    use scale_process, only: &
       prc_mpistop
    use scale_const, only: &
       cpvap          => const_cpvap,          &
       cvvap          => const_cvvap,          &
       cl             => const_cl,             &
       ci             => const_ci,             &
       thermodyn_type => const_thermodyn_type
    implicit none

    integer :: n
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[THERMODYN] / Categ[ATMOS SHARE] / Origin[SCALElib]'

    allocate( aq_cp(qqa) )
    allocate( aq_cv(qqa) )

    if ( thermodyn_type == 'EXACT' ) then
       thermodyn_emask = 1.0_rp

       aq_cp(i_qv) = cpvap
       aq_cv(i_qv) = cvvap

       if ( qws /= 0 ) then
          do n = qws, qwe
             aq_cp(n) = cl
             aq_cv(n) = cl
          enddo
       endif

       if ( qis /= 0 ) then
          do n = qis, qie
             aq_cp(n) = ci
             aq_cv(n) = ci
          enddo
       endif
    elseif( thermodyn_type == 'SIMPLE' ) then
       thermodyn_emask = 0.0_rp

       aq_cp(i_qv) = cpdry
       aq_cv(i_qv) = cvdry

       if ( qws /= 0 ) then
          do n = qws, qwe
             aq_cp(n) = cvdry
             aq_cv(n) = cvdry
          enddo
       endif

       if ( qis /= 0 ) then
          do n = qis, qie
             aq_cp(n) = cvdry
             aq_cv(n) = cvdry
          enddo
       endif
    endif

    return
  end subroutine atmos_thermodyn_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_qd_0d( &
       qdry, &
       q     )
    implicit none

    real(RP), intent(out) :: qdry
    real(RP), intent(in)  :: q(qa)

    integer :: iqw
    !-----------------------------------------------------------------------------

    qdry = 1.0_rp
#ifndef DRY
    do iqw = qqs, qqe
       qdry = qdry - q(iqw)
    enddo
#endif

    return
  end subroutine atmos_thermodyn_qd_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_qd_3d( &
       qdry, &
       q     )
    implicit none

    real(RP), intent(out) :: qdry(ka,ia,ja)
    real(RP), intent(in)  :: q   (ka,ia,ja,qa)

    integer :: k, i, j, iqw
    !-----------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

!       qdry(k,i,j) = 1.0_RP
!       do iqw = QQS, QQE
!          qdry(k,i,j) = qdry(k,i,j) - q(k,i,j,iqw)
!       enddo
       calc_qdry( qdry(k,i,j), q, k, i, j, iqw )

    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_qd_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_cp_0d( &
       CPtot, &
       q,     &
       qdry   )
    implicit none

    real(RP), intent(out) :: CPtot
    real(RP), intent(in)  :: q(qa)
    real(RP), intent(in)  :: qdry

    integer :: iqw
    !---------------------------------------------------------------------------

    cptot = qdry * cpdry
#ifndef DRY
    do iqw = qqs, qqe
       cptot = cptot + q(iqw) * aq_cp(iqw)
    enddo
#endif

    return
  end subroutine atmos_thermodyn_cp_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_cp_3d( &
       CPtot, &
       q,     &
       qdry   )
    implicit none

    real(RP), intent(out) :: CPtot(ka,ia,ja)
    real(RP), intent(in)  :: q    (ka,ia,ja,qa)
    real(RP), intent(in)  :: qdry (ka,ia,ja)

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

!       CPtot(k,i,j) = qdry(k,i,j) * CPdry
!       do iqw = QQS, QQE
!          CPtot(k,i,j) = CPtot(k,i,j) + q(k,i,j,iqw) * AQ_CP(iqw)
!       enddo

       calc_cp(cptot(k,i,j), qdry(k,i,j), q, k, i, j, iqw, cpdry, aq_cp)

    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_cp_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_cv_0d( &
       CVtot, &
       q,     &
       qdry   )
    implicit none

    real(RP), intent(out) :: CVtot
    real(RP), intent(in)  :: q(qa)
    real(RP), intent(in)  :: qdry

    integer :: iqw
    !---------------------------------------------------------------------------

    cvtot = qdry * cvdry
#ifndef DRY
    do iqw = qqs, qqe
       cvtot = cvtot + q(iqw) * aq_cv(iqw)
    enddo
#endif

    return
  end subroutine atmos_thermodyn_cv_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_cv_3d( &
       CVtot, &
       q,     &
       qdry   )
    implicit none

    real(RP), intent(out) :: CVtot(ka,ia,ja)
    real(RP), intent(in)  :: q    (ka,ia,ja,qa)
    real(RP), intent(in)  :: qdry (ka,ia,ja)

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

!       CVtot(k,i,j) = qdry(k,i,j) * CVdry
!       do iqw = QQS, QQE
!          CVtot(k,i,j) = CVtot(k,i,j) + q(k,i,j,iqw) * AQ_CV(iqw)
!       enddo

       calc_cv(cvtot(k,i,j), qdry(k,i,j), q, k, i, j, iqw, cvdry, aq_cv)

    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_cv_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_pott_0d( &
       pott, &
       temp, &
       pres, &
       q     )
    implicit none

    real(RP), intent(out) :: pott  ! potential temperature [K]
    real(RP), intent(in)  :: temp  ! temperature           [K]
    real(RP), intent(in)  :: pres  ! pressure              [Pa]
    real(RP), intent(in)  :: q(qa) ! mass concentration    [kg/kg]

    real(RP) :: qdry
    real(RP) :: Rtot, CVtot, RovCP

    integer  :: iqw
    !---------------------------------------------------------------------------

#ifdef DRY
    cvtot = cvdry
    rtot  = rdry
#else
    qdry  = 1.0_rp
    cvtot = 0.0_rp
    do iqw = qqs, qqe
       qdry  = qdry  - q(iqw)
       cvtot = cvtot + q(iqw) * aq_cv(iqw)
    enddo
    cvtot = cvdry * qdry + cvtot
    rtot  = rdry  * qdry + rvap * q(i_qv)
#endif

    rovcp = rtot / ( cvtot + rtot )

    pott  = temp * ( pre00 / pres )**rovcp

    return
  end subroutine atmos_thermodyn_pott_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_pott_3d( &
       pott, &
       temp, &
       pres, &
       q     )
    implicit none

    real(RP), intent(out) :: pott(ka,ia,ja)    ! potential temperature [K]
    real(RP), intent(in)  :: temp(ka,ia,ja)    ! temperature           [K]
    real(RP), intent(in)  :: pres(ka,ia,ja)    ! pressure              [Pa]
    real(RP), intent(in)  :: q   (ka,ia,ja,qa) ! mass concentration    [kg/kg]

    real(RP) :: qdry
    real(RP) :: Rtot, CVtot, RovCP

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
#ifdef DRY
       cvtot = cvdry
       rtot  = rdry
#else
       qdry  = 1.0_rp
       cvtot = 0.0_rp
       do iqw = qqs, qqe
          qdry  = qdry  - q(k,i,j,iqw)
          cvtot = cvtot + q(k,i,j,iqw) * aq_cv(iqw)
       enddo
       cvtot = cvdry * qdry + cvtot
       rtot  = rdry  * qdry + rvap * q(k,i,j,i_qv)
#endif

       rovcp = rtot / ( cvtot + rtot )

       pott(k,i,j) = temp(k,i,j) * ( pre00 / pres(k,i,j) )**rovcp
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_pott_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_rhoe_0d( &
       rhoe, &
       rhot, &
       q     )
    implicit none

    real(RP), intent(out) :: rhoe
    real(RP), intent(in)  :: rhot
    real(RP), intent(in)  :: q(qa)

    real(RP) :: qdry
    real(RP) :: pres
    real(RP) :: Rtot, CVtot, CPovCV

    integer :: iqw
    !---------------------------------------------------------------------------

#ifdef DRY
    cvtot = cvdry
    rtot  = rdry
#else
    qdry  = 1.0_rp
    cvtot = 0.0_rp
    do iqw = qqs, qqe
       qdry  = qdry  - q(iqw)
       cvtot = cvtot + q(iqw) * aq_cv(iqw)
    enddo
    cvtot = cvdry * qdry + cvtot
    rtot  = rdry  * qdry + rvap * q(i_qv)
#endif

    cpovcv = ( cvtot + rtot ) / cvtot

    pres = pre00 * ( rhot * rtot / pre00 )**cpovcv

    rhoe = pres / rtot * cvtot

    return
  end subroutine atmos_thermodyn_rhoe_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_rhoe_3d( &
       rhoe, &
       rhot, &
       q     )
    implicit none

    real(RP), intent(out) :: rhoe(ka,ia,ja)
    real(RP), intent(in)  :: rhot(ka,ia,ja)
    real(RP), intent(in)  :: q   (ka,ia,ja,qa)

    real(RP) :: qdry
    real(RP) :: pres
    real(RP) :: Rtot, CVtot, CPovCV

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw,qdry,pres,Rtot,CVtot,CPovCV) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
#ifdef DRY
       cvtot = cvdry
       rtot  = rdry
#else
       qdry  = 1.0_rp
       cvtot = 0.0_rp
       do iqw = qqs, qqe
          qdry  = qdry  - q(k,i,j,iqw)
          cvtot = cvtot + q(k,i,j,iqw) * aq_cv(iqw)
       enddo
       cvtot = cvdry * qdry + cvtot
       rtot  = rdry  * qdry + rvap * q(k,i,j,i_qv)
#endif


       cpovcv = ( cvtot + rtot ) / cvtot

       pres = pre00 * ( rhot(k,i,j) * rtot / pre00 )**cpovcv

       rhoe(k,i,j) = pres / rtot * cvtot
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_rhoe_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_rhot_0d( &
       rhot, &
       rhoe, &
       q     )
    implicit none

    real(RP), intent(out) :: rhot
    real(RP), intent(in)  :: rhoe
    real(RP), intent(in)  :: q(qa)

    real(RP) :: qdry
    real(RP) :: pres
    real(RP) :: Rtot, CVtot, RovCP

    integer :: iqw
    !---------------------------------------------------------------------------

#ifdef DRY
    cvtot = cvdry
    rtot  = rdry
#else
    qdry  = 1.0_rp
    cvtot = 0.0_rp
    do iqw = qqs, qqe
       qdry  = qdry  - q(iqw)
       cvtot = cvtot + q(iqw) * aq_cv(iqw)
    enddo
    cvtot = cvdry * qdry + cvtot
    rtot  = rdry  * qdry + rvap * q(i_qv)
#endif

    rovcp = rtot / ( cvtot + rtot )

    pres = rhoe * rtot / cvtot

    rhot = rhoe / cvtot * ( pre00 / pres )**rovcp

    return
  end subroutine atmos_thermodyn_rhot_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_rhot_3d( &
       rhot, &
       rhoe, &
       q     )
    implicit none

    real(RP), intent(out) :: rhot(ka,ia,ja)
    real(RP), intent(in)  :: rhoe(ka,ia,ja)
    real(RP), intent(in)  :: q   (ka,ia,ja,qa)

    real(RP) :: qdry
    real(RP) :: pres
    real(RP) :: Rtot, CVtot, RovCP

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw,qdry,pres,Rtot,CVtot,RovCP) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

#ifdef DRY
       cvtot = cvdry
       rtot  = rdry
#else
       qdry  = 1.0_rp
       cvtot = 0.0_rp
       do iqw = qqs, qqe
          qdry  = qdry  - q(k,i,j,iqw)
          cvtot = cvtot + q(k,i,j,iqw) * aq_cv(iqw)
       enddo
       cvtot = cvdry * qdry + cvtot
       rtot  = rdry  * qdry + rvap * q(k,i,j,i_qv)
#endif

       rovcp = rtot / ( cvtot + rtot )

       pres = rhoe(k,i,j) * rtot / cvtot

       rhot(k,i,j) = rhoe(k,i,j) / cvtot * ( pre00 / pres )**rovcp
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_rhot_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_temp_pres_0d( &
       temp, &
       pres, &
       dens, &
       rhot, &
       q     )
    implicit none

    real(RP), intent(out) :: temp
    real(RP), intent(out) :: pres
    real(RP), intent(in)  :: dens
    real(RP), intent(in)  :: rhot
    real(RP), intent(in)  :: q(qa)

    real(RP) :: qdry
    real(RP) :: Rtot, CVtot, CPovCV

    integer :: iqw
    !---------------------------------------------------------------------------

#ifdef DRY
    cvtot = cvdry
    rtot  = rdry
#else
    qdry  = 1.0_rp
    cvtot = 0.0_rp
    do iqw = qqs, qqe
       qdry  = qdry  - q(iqw)
       cvtot = cvtot + q(iqw) * aq_cv(iqw)
    enddo
    cvtot = cvdry * qdry + cvtot
    rtot  = rdry  * qdry + rvap * q(i_qv)
#endif

    cpovcv = ( cvtot + rtot ) / cvtot

    pres = pre00 * ( rhot * rtot / pre00 )**cpovcv
    temp = pres / ( dens * rtot )

    return
  end subroutine atmos_thermodyn_temp_pres_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_temp_pres_3d( &
       temp, &
       pres, &
       dens, &
       rhot, &
       q     )
    implicit none

    real(RP), intent(out) :: temp(ka,ia,ja)
    real(RP), intent(out) :: pres(ka,ia,ja)
    real(RP), intent(in)  :: dens(ka,ia,ja)
    real(RP), intent(in)  :: rhot(ka,ia,ja)
    real(RP), intent(in)  :: q   (ka,ia,ja,qa)

    real(RP) :: qdry
    real(RP) :: Rtot, CVtot, CPovCV

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw,qdry,Rtot,CVtot,CPovCV) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
#ifdef DRY
       cvtot = cvdry
       rtot  = rdry
#else
       qdry  = 1.0_rp
       cvtot = 0.0_rp
       do iqw = qqs, qqe
          qdry  = qdry  - q(k,i,j,iqw)
          cvtot = cvtot + q(k,i,j,iqw) * aq_cv(iqw)
       enddo
       cvtot = cvdry * qdry + cvtot
       rtot  = rdry  * qdry + rvap * q(k,i,j,i_qv)
#endif

       cpovcv = ( cvtot + rtot ) / cvtot

       pres(k,i,j) = pre00 * ( rhot(k,i,j) * rtot / pre00 )**cpovcv
       temp(k,i,j) = pres(k,i,j) / ( dens(k,i,j) * rtot )
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_temp_pres_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_temp_pres_e_0d( &
       temp, &
       pres, &
       dens, &
       rhoe, &
       q     )
    implicit none

    real(RP), intent(out) :: temp
    real(RP), intent(out) :: pres
    real(RP), intent(in)  :: dens
    real(RP), intent(in)  :: rhoe
    real(RP), intent(in)  :: q(qa)

    real(RP) :: qdry
    real(RP) :: Rtot, CVtot

    integer :: iqw
    !---------------------------------------------------------------------------

#ifdef DRY
    cvtot = cvdry
    rtot  = rdry
#else
    qdry  = 1.0_rp
    cvtot = 0.0_rp
    do iqw = qqs, qqe
       qdry  = qdry  - q(iqw)
       cvtot = cvtot + q(iqw) * aq_cv(iqw)
    enddo
    cvtot = cvdry * qdry + cvtot
    rtot  = rdry  * qdry + rvap * q(i_qv)
#endif

    temp = rhoe / ( dens * cvtot )
    pres = dens * rtot * temp

    return
  end subroutine atmos_thermodyn_temp_pres_e_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_temp_pres_e_3d( &
       temp, &
       pres, &
       dens, &
       rhoe, &
       q     )
    implicit none

    real(RP), intent(out) :: temp(ka,ia,ja)
    real(RP), intent(out) :: pres(ka,ia,ja)
    real(RP), intent(in)  :: dens(ka,ia,ja)
    real(RP), intent(in)  :: rhoe(ka,ia,ja)
    real(RP), intent(in)  :: q   (ka,ia,ja,qa)

    real(RP) :: qdry
    real(RP) :: Rtot, CVtot

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw,qdry,Rtot,CVtot) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
#ifdef DRY
       cvtot = cvdry
       rtot  = rdry
#else
       qdry  = 1.0_rp
       cvtot = 0.0_rp
       do iqw = qqs, qqe
          qdry  = qdry  - q(k,i,j,iqw)
          cvtot = cvtot + q(k,i,j,iqw) * aq_cv(iqw)
       enddo
       cvtot = cvdry * qdry + cvtot
       rtot  = rdry  * qdry + rvap * q(k,i,j,i_qv)
#endif

       temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cvtot )
       pres(k,i,j) = dens(k,i,j) * rtot * temp(k,i,j)
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_temp_pres_e_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhv_0d( &
       lhv, &
       temp )
    implicit none

    real(RP), intent(out) :: lhv  ! potential temperature [K]
    real(RP), intent(in)  :: temp ! temperature           [K]
    !---------------------------------------------------------------------------

    lhv = lhv0 + ( cpvap - cl ) * ( temp - tem00 ) * thermodyn_emask

    return
  end subroutine atmos_thermodyn_templhv_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhv_2d( &
       lhv, &
       temp )
    implicit none

    real(RP), intent(out) :: lhv (ia,ja) ! potential temperature [K]
    real(RP), intent(in)  :: temp(ia,ja) ! temperature           [K]

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
       lhv(i,j) = lhv0 + ( cpvap - cl ) * ( temp(i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_templhv_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhv_3d( &
       lhv, &
       temp )
    implicit none

    real(RP), intent(out) :: lhv (ka,ia,ja) ! potential temperature [K]
    real(RP), intent(in)  :: temp(ka,ia,ja) ! temperature           [K]

    integer :: k, i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhv(k,i,j) = lhv0 + ( cpvap - cl ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_templhv_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhs_0d( &
       lhs, &
       temp )
    implicit none

    real(RP), intent(out) :: lhs   ! potential temperature [K]
    real(RP), intent(in)  :: temp  ! temperature           [K]
    !---------------------------------------------------------------------------

    lhs = lhs0 + ( cpvap - ci ) * ( temp - tem00 ) * thermodyn_emask

    return
  end subroutine atmos_thermodyn_templhs_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhs_2d( &
       lhs, &
       temp )
    implicit none

    real(RP), intent(out) :: lhs (ia,ja) ! potential temperature [K]
    real(RP), intent(in)  :: temp(ia,ja) ! temperature           [K]

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
       lhs(i,j) = lhs0 + ( cpvap - ci ) * ( temp(i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_templhs_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhs_3d( &
       lhs, &
       temp )
    implicit none

    real(RP), intent(out) :: lhs (ka,ia,ja) ! potential temperature [K]
    real(RP), intent(in)  :: temp(ka,ia,ja) ! temperature           [K]

    integer :: k, i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhs(k,i,j) = lhs0 + ( cpvap - ci ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_templhs_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhf_0d( &
       lhf, &
       temp )
    implicit none

    real(RP), intent(out) :: lhf  ! potential temperature [K]
    real(RP), intent(in)  :: temp ! temperature           [K]
    !---------------------------------------------------------------------------

    lhf = lhf0 + ( cl - ci ) * ( temp - tem00 ) * thermodyn_emask

    return
  end subroutine atmos_thermodyn_templhf_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhf_2d( &
       lhf, &
       temp )
    implicit none

    real(RP), intent(out) :: lhf (ia,ja) ! potential temperature [K]
    real(RP), intent(in)  :: temp(ia,ja) ! temperature           [K]

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
       lhf(i,j) = lhf0 + ( cl - ci ) * ( temp(i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_templhf_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_templhf_3d( &
       lhf, &
       temp )
    implicit none

    real(RP), intent(out) :: lhf (ka,ia,ja) ! potential temperature [K]
    real(RP), intent(in)  :: temp(ka,ia,ja) ! temperature           [K]

    integer :: k, i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhf(k,i,j) = lhf0 + ( cl - ci ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_templhf_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_entr_0d( &
       entr, &
       temp, &
       pres, &
       q     )
    implicit none

    real(RP), intent(out) :: entr
    real(RP), intent(in)  :: temp
    real(RP), intent(in)  :: pres
    real(RP), intent(in)  :: q(qa)

    real(RP) :: lhv, lhf
    real(RP) :: qdry, Rtot
    real(RP) :: logT_T0, pres_dry, pres_vap
    real(RP) :: sq

    integer :: iqw
    !---------------------------------------------------------------------------

    lhv = lhv0 + ( cpvap - cl ) * ( temp - tem00 ) * thermodyn_emask
    lhf = lhf0 + ( cl    - ci ) * ( temp - tem00 ) * thermodyn_emask

    qdry = 1.0_rp
    do iqw = qqs, qqe
       qdry = qdry - q(iqw)
    enddo
    rtot = rdry * qdry + rvap * q(i_qv)

    logt_t0 = log( temp / tem00 )

    ! dry air + vapor
    pres_dry = max( pres * qdry    * rdry / rtot, eps )
    pres_vap = max( pres * q(i_qv) * rvap / rtot, eps )

    entr = qdry    * cpdry * logt_t0                 &
         - qdry    * rdry  * log( pres_dry / pre00 ) &
         + q(i_qv) * cpvap * logt_t0                 &
         - q(i_qv) * rvap  * log( pres_vap / psat0 ) &
         + q(i_qv) * lhv / tem00

    ! hydrometeors
    do iqw = qqs, qqe
       sq = 0.0_rp

       if( iqw == i_qc ) sq = q(i_qc) * cl * logt_t0

       if( iqw == i_qr ) sq = q(i_qr) * cl * logt_t0

       if( iqw == i_qi ) sq = q(i_qi) * ci * logt_t0 &
                            - q(i_qi) * lhf / tem00

       if( iqw == i_qs ) sq = q(i_qs) * ci * logt_t0 &
                            - q(i_qs) * lhf / tem00

       if( iqw == i_qg ) sq = q(i_qg) * ci * logt_t0 &
                            - q(i_qg) * lhf / tem00

       entr = entr + sq
    enddo

    return
  end subroutine atmos_thermodyn_entr_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_entr_2d( &
       entr, &
       temp, &
       pres, &
       q     )
    implicit none

    real(RP), intent(out) :: entr(ia,ja)
    real(RP), intent(in)  :: temp(ia,ja)
    real(RP), intent(in)  :: pres(ia,ja)
    real(RP), intent(in)  :: q   (ia,ja,qa)

    real(RP) :: lhv, lhf
    real(RP) :: qdry, Rtot
    real(RP) :: logT_T0, pres_dry, pres_vap
    real(RP) :: sq

    integer :: i, j, iqw
    !---------------------------------------------------------------------------

    ! dry air + vapor
    do j = jsb, jeb
    do i = isb, ieb
       lhv = lhv0 + ( cpvap - cl ) * ( temp(i,j) - tem00 ) * thermodyn_emask

       qdry = 1.0_rp
       do iqw = qqs, qqe
          qdry = qdry - q(i,j,iqw)
       enddo
       rtot = rdry * qdry + rvap * q(i,j,i_qv)

       logt_t0 = log( temp(i,j) / tem00 )

       pres_dry = max( pres(i,j) * qdry        * rdry / rtot, eps )
       pres_vap = max( pres(i,j) * q(i,j,i_qv) * rvap / rtot, eps )

       entr(i,j) = qdry        * cpdry * logt_t0 &
                 - qdry        * rdry  * log( pres_dry / pre00 ) &
                 + q(i,j,i_qv) * cpvap * logt_t0 &
                 - q(i,j,i_qv) * rvap  * log( pres_vap / psat0 ) &
                 + q(i,j,i_qv) * lhv / tem00
    enddo
    enddo

    ! hydrometeors
    do j = jsb, jeb
    do i = isb, ieb
       lhf = lhf0 + ( cl - ci ) * ( temp(i,j) - tem00 ) * thermodyn_emask

       logt_t0 = log( temp(i,j) / tem00 )

       do iqw = qqs, qqe
          sq = 0.0_rp

          if( iqw == i_qc ) sq = q(i,j,i_qc) *   cl * logt_t0
          if( iqw == i_qr ) sq = q(i,j,i_qr) *   cl * logt_t0
          if( iqw == i_qi ) sq = q(i,j,i_qi) * ( ci * logt_t0 - lhf / tem00 )
          if( iqw == i_qs ) sq = q(i,j,i_qs) * ( ci * logt_t0 - lhf / tem00 )
          if( iqw == i_qg ) sq = q(i,j,i_qg) * ( ci * logt_t0 - lhf / tem00 )

          entr(i,j) = entr(i,j) + sq
       enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_entr_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_entr_3d( &
       entr, &
       temp, &
       pres, &
       q     )
    implicit none

    real(RP), intent(out) :: entr(ka,ia,ja)
    real(RP), intent(in)  :: temp(ka,ia,ja)
    real(RP), intent(in)  :: pres(ka,ia,ja)
    real(RP), intent(in)  :: q   (ka,ia,ja,qa)

    real(RP) :: lhv, lhf
    real(RP) :: qdry, Rtot
    real(RP) :: logT_T0, pres_dry, pres_vap
    real(RP) :: sq

    integer :: k, i, j, iqw
    !---------------------------------------------------------------------------

    ! dry air + vapor
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhv = lhv0 + ( cpvap - cl ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask

       qdry = 1.0_rp
       do iqw = qqs, qqe
          qdry = qdry - q(k,i,j,iqw)
       enddo
       rtot = rdry * qdry + rvap * q(k,i,j,i_qv)

       logt_t0 = log( temp(k,i,j) / tem00 )

       pres_dry = max( pres(k,i,j) * qdry          * rdry / rtot, eps )
       pres_vap = max( pres(k,i,j) * q(k,i,j,i_qv) * rvap / rtot, eps )

       entr(k,i,j) = qdry          * cpdry * logt_t0 &
                   - qdry          * rdry  * log( pres_dry / pre00 ) &
                   + q(k,i,j,i_qv) * cpvap * logt_t0 &
                   - q(k,i,j,i_qv) * rvap  * log( pres_vap / psat0 ) &
                   + q(k,i,j,i_qv) * lhv / tem00
    enddo
    enddo
    enddo

    ! hydrometeors
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhf = lhf0 + ( cl - ci ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask

       logt_t0 = log( temp(k,i,j) / tem00 )

       do iqw = qqs, qqe
          sq = 0.0_rp

          if( iqw == i_qc ) sq = q(k,i,j,i_qc) *   cl * logt_t0
          if( iqw == i_qr ) sq = q(k,i,j,i_qr) *   cl * logt_t0
          if( iqw == i_qi ) sq = q(k,i,j,i_qi) * ( ci * logt_t0 - lhf / tem00 )
          if( iqw == i_qs ) sq = q(k,i,j,i_qs) * ( ci * logt_t0 - lhf / tem00 )
          if( iqw == i_qg ) sq = q(k,i,j,i_qg) * ( ci * logt_t0 - lhf / tem00 )

          entr(k,i,j) = entr(k,i,j) + sq
       enddo
    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_entr_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_tempre( &
      temp, pres,         &
      Ein,  dens, qdry, q )
    implicit none

    real(RP), intent(out) :: temp(ka,ia,ja)    ! temperature
    real(RP), intent(out) :: pres(ka,ia,ja)    ! pressure
    real(RP), intent(in)  :: Ein (ka,ia,ja)    ! internal energy
    real(RP), intent(in)  :: dens(ka,ia,ja)    ! density
    real(RP), intent(in)  :: qdry(ka,ia,ja)    ! dry concentration
    real(RP), intent(in)  :: q   (ka,ia,ja,qa) ! water concentration

    real(RP) :: cv, Rmoist

    integer :: i, j, k, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw,cv,Rmoist) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

       calc_cv(cv, qdry(k,i,j), q, k, i, j, iqw, cvdry, aq_cv)
       calc_r(rmoist, q(k,i,j,i_qv), qdry(k,i,j), rdry, rvap)

       temp(k,i,j) = ein(k,i,j) / cv

       pres(k,i,j) = dens(k,i,j) * rmoist * temp(k,i,j)

    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_tempre

  !-----------------------------------------------------------------------------
  subroutine atmos_thermodyn_tempre2( &
      temp, pres,         &
      dens, pott, qdry, q )
    implicit none

    real(RP), intent(out) :: temp(ka,ia,ja)    ! temperature
    real(RP), intent(out) :: pres(ka,ia,ja)    ! pressure
    real(RP), intent(in)  :: dens(ka,ia,ja)    ! density
    real(RP), intent(in)  :: pott(ka,ia,ja)    ! potential temperature
    real(RP), intent(in)  :: qdry(ka,ia,ja)    ! dry concentration
    real(RP), intent(in)  :: q   (ka,ia,ja,qa) ! water concentration

    real(RP) :: Rmoist, cp

    integer :: i, j, k, iqw
    !---------------------------------------------------------------------------

    !$omp parallel do private(i,j,k,iqw,cp,Rmoist) OMP_SCHEDULE_ collapse(2)
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

       calc_cp(cp, qdry(k,i,j), q, k, i, j, iqw, cpdry, aq_cp)
       calc_r(rmoist, q(k,i,j,i_qv), qdry(k,i,j), rdry, rvap)
       calc_pre(pres(k,i,j), dens(k,i,j), pott(k,i,j), rmoist, cp, pre00)

       temp(k,i,j) = pres(k,i,j) / ( dens(k,i,j) * rmoist )

    enddo
    enddo
    enddo

    return
  end subroutine atmos_thermodyn_tempre2

end module scale_atmos_thermodyn