Skip to content

Commit

Permalink
fix the moist adiabat calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@abtawfik
abtawfik committed Nov 16, 2016
1 parent cc04f37 commit 964101d
Showing 1 changed file with 83 additions and 13 deletions.
96 changes: 83 additions & 13 deletions convect_trigger_pot/ncl_portable/ctp_hilow.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,6 @@
! Converted to F90 code -- Joshua Roundy on Apr 2015
! Modified by -- A.B. Tawfik on June 2015
!
!---------------------------------------------------------------------------------



!---------------------------------------------------------------------------------
!
! subroutines: Calculates the CTP and HiLow
! Assumes left most dimension is level
Expand All @@ -53,7 +48,8 @@
! in mass.
!
!---------------------------------------------------------------------------------
subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, t2m_in , q2m_in , psfc_in, CTP, HILOW , missing )
subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, &
t2m_in , q2m_in , psfc_in, CTP, HILOW , missing )

implicit none
!
Expand Down Expand Up @@ -110,6 +106,7 @@ subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, t2m_in , q2m_in , ps
real(4) :: moist_lapse, dz
real(4) :: pmid, tmid, qmid
real(4) :: qseg_old
real(4) :: qsat
!-----------------------------------------------------------------------------


Expand Down Expand Up @@ -145,9 +142,9 @@ subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, t2m_in , q2m_in , ps
!-- indices below and above
!-- desired pressure level
!--------------------------------
do nn=1,nlev
level_index: do nn=1,nlev
ilev(nn) = real(nn)
end do
end do level_index

!-----------------------------------------------------------------------------------
!-- Make sure there are no missing critical inputs --> if so return CTP HILOW missing
Expand Down Expand Up @@ -299,7 +296,7 @@ subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, t2m_in , q2m_in , ps
tpar_old = temp100
qseg_old = qhum100
CTP = 0.
do nn=1,nsegments
ctp_depth: do nn=1,nsegments

!----------------------------------------------------
!-- Pressure increment inbetween defined levels (Pa)
Expand Down Expand Up @@ -331,10 +328,12 @@ subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, t2m_in , q2m_in , ps
tmid = ( (t_segment*log(p_segment) + tseg_old*log(p_old)) / log(p_segment*p_old) )
qmid = ( (q_segment*log(p_segment) + qseg_old*log(p_old)) / log(p_segment*p_old) )
dz = (p_old-p_segment) / (grav * pmid /(Rd*tmid*((1.+(qmid/ep)) / (1. + qmid))))

moist_lapse = (grav/cp) * ( (1. + (Lv * qmid)/( Rd*tmid )) / &
(1. + (Lv**2 * qmid)/(cp*Rv*tmid**2)) )
qsat = saturation_specific_humidity( pmid, tmid, missing)

moist_lapse = (grav/cp) * ( (1. + (Lv * qsat)/( Rd*tmid )) / &
(1. + (Lv**2 * qsat)/(cp*Rv*tmid**2)) )


!----------------------------------------------------
!--- Get the parcel temperature
!----------------------------------------------------
Expand All @@ -359,7 +358,7 @@ subroutine ctp_hi_low ( nlev_in, tlev_in, qlev_in, plev_in, t2m_in , q2m_in , ps
qseg_old = q_segment
p_old = p_segment

end do
end do ctp_depth


return
Expand Down Expand Up @@ -418,3 +417,74 @@ subroutine dew_point( qlev, plev, tdew, missing )

end subroutine dew_point
!---------------------------------------------------------------------------------









!---------------------------------------------------------------------------------
!
! subroutines: Calculates the saturation specific humidity [kg/kg]
! following the AMS glossary definition
!
!---------------------------------------------------------------------------------
real(4) function saturation_specific_humidity( p, t, missing )

implicit none

!
! Input/Output Variables
!
real(4), intent(in ) :: missing !** missing value
real(4), intent(in ) :: p !** pressure [Pa]
real(4), intent(in ) :: t !** dry bulb temp. [K]
!
! Local variables
!
real(4) :: press, esat
real(4) :: numerator, denomenator
real(4), parameter :: t0 = 273.15, ep=0.622, es0=6.11, a=17.269, b=35.86
real(4), parameter :: onemep= 1.0 - ep

!---------------------------------------------------------------------------------

!--------------------------------------------
!--- Initialization and preliminary calculations
!--------------------------------------------
saturation_specific_humidity = missing

!--------------------------------------------
!--- Perform a quick check for missing values
!--------------------------------------------
if( t.eq.missing .or. p.eq.missing ) return

!--------------------------------------------
!--- Convert pressure to hPa
!--------------------------------------------
press = p/1e2

!--------------------------------------------
!--- Split out the numerator and denomenator
!--------------------------------------------
numerator = ep* (es0*exp((a*( t-t0))/( t-b)))
denomenator = press-onemep*(es0*exp((a*( t-t0))/( t-b)))

!--------------------------------------------
!--- Intermediate calculation
!--------------------------------------------
esat = numerator/denomenator

!--------------------------------------------
!--- Vapor pressure and co
!--------------------------------------------
saturation_specific_humidity = esat / (1 + esat)

end function saturation_specific_humidity
!---------------------------------------------------------------------------------



0 comments on commit 964101d

Please sign in to comment.