Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add cloud top temperature and pressure estimates from NCL #1069

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Add cloud top temperature and pressure estimates from NCL #1069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/core_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -951,6 +951,8 @@
<var name="q2"/>
<var name="t2m"/>
<var name="th2m"/>
<var name="ctt"/>
<var name="ctpres"/>
<var name="gsw"/>
<var name="glw"/>
<var name="acsnow"/>
Expand Down
3 changes: 2 additions & 1 deletion src/core_atmosphere/diagnostics/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,13 @@ DIAGNOSTIC_MODULES = \
mpas_convective_diagnostics.o \
mpas_pv_diagnostics.o \
mpas_soundings.o \
calc_ctt.o \

mpas_isobaric_diagnostics.o: mpas_atm_diagnostics_utils.o

mpas_cloud_diagnostics.o: mpas_atm_diagnostics_utils.o

mpas_convective_diagnostics.o: mpas_atm_diagnostics_utils.o
mpas_convective_diagnostics.o: mpas_atm_diagnostics_utils.o calc_ctt.o

mpas_pv_diagnostics.o: mpas_atm_diagnostics_utils.o

Expand Down
6 changes: 6 additions & 0 deletions src/core_atmosphere/diagnostics/Registry_convective.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,5 +55,11 @@
<var name="wind_speed_level1_max" type="real" dimensions="nCells Time" units="m s^{-1}"
description="Maximum wind speed in lowest model level since last output"/>

<var name="ctt" type="real" dimensions="nCells Time" units="deg C" missing_value="-9999."
description="Cloud top temperature derived from NCL routines"/>

<var name="ctpres" type="real" dimensions="nCells Time" units="hPa"
description="Cloud top pressure derived from NCL routines"/>

</var_struct>

101 changes: 101 additions & 0 deletions src/core_atmosphere/diagnostics/calc_ctt.F
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
module calculate_ctt

contains

SUBROUTINE get_ctt(prs, tk, qci, qcw, ctt, haveqci,&
fill_nocloud, missing, opt_thresh, nz, prsctt)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nz, haveqci, fill_nocloud
REAL(KIND=8), DIMENSION(nz), INTENT(IN) :: prs, tk, qci, qcw
REAL, INTENT(OUT) :: ctt, prsctt
REAL(KIND=8), INTENT(IN) :: missing
REAL(KIND=8), INTENT(IN) :: opt_thresh
REAL(KIND=8), DIMENSION(nz):: pf

! LOCAL VARIABLES
INTEGER i,j,k,ripk
REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, ratmix
REAL(KIND=8) :: arg2, agl_hgt, vt

REAL(KIND=8) :: p1, p2

REAL(KIND=8), PARAMETER :: ABSCOEFI = .17D0 ! cloud ice absorption coefficient in m^2/g , was 0.272 as default

! https://www.ecmwf.int/sites/default/files/elibrary/2019/19056-improved-ice-cloud-modeling-capabilities-community-radiative-transfer-model.pdf
REAL(KIND=8), PARAMETER :: ABSCOEF = .145D0 ! cloud water absorption coefficient in m^2/g
REAL(KIND=8), PARAMETER :: CELKEL = 273.15D0 ! cloud water absorption coefficient in m^2/g
REAL(KIND=8), parameter :: gravity = 9.80616D0 !< Constant: Acceleration due to gravity [m s-2]

DO k=1,nz-1
ripk = nz-k+1
pf(k) = .5D0*(prs(ripk) + prs(ripk-1))
END DO

opdepthd = 0.D0
k = 0
prsctt = -1

! Integrate downward from model top, calculating path at full
! model vertical levels.

DO k=2,nz
opdepthu = opdepthd
ripk = nz - k + 1

IF (k .NE. 1) THEN
dp = 100.D0*(pf(k) - pf(k-1)) ! should be in Pa
ELSE
dp = 200.D0*(pf(1) - prs(nz)) ! should be in Pa
END IF

IF (haveqci .EQ. 0) then
IF (tk(ripk) .LT. CELKEL) then
! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed
opdepthd = opdepthu + ABSCOEFI*qcw(ripk)*1000.D0 * dp/gravity
ELSE
opdepthd = opdepthu + ABSCOEF*qcw(ripk)*1000.D0 * dp/gravity
END IF
ELSE
opdepthd = opdepthd + (ABSCOEF*qcw(ripk)*1000.D0 + ABSCOEFI*qci(ripk)*1000.D0)*dp/gravity
END IF

IF (opdepthd .LT. opt_thresh .AND. k .LT. nz) THEN
CYCLE

ELSE IF (opdepthd .LT. opt_thresh .AND. k .EQ. nz) THEN
IF (fill_nocloud .EQ. 0) THEN
prsctt = prs(1)
ENDIF
EXIT
ELSE
fac = (1. - opdepthu)/(opdepthd - opdepthu)
prsctt = pf(k-1) + fac*(pf(k) - pf(k-1))
prsctt = REAL(MIN(prs(1), MAX(prs(nz), prsctt)))
EXIT
END IF
END DO

! prsctt should only be 0 if fill values are used
IF (prsctt .GT. -1) THEN
DO k=2,nz
ripk = nz - k + 1
p1 = prs(ripk+1)
p2 = prs(ripk)
IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN
fac = (prsctt - p1)/(p2 - p1)
arg1 = fac*(tk(ripk) - tk(ripk+1)) - CELKEL
ctt = REAL(tk(ripk+1) + arg1)
EXIT
END IF
END DO
ELSE
ctt = REAL(missing)
END IF

RETURN

END SUBROUTINE get_ctt

end module calculate_ctt
49 changes: 46 additions & 3 deletions src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ module mpas_convective_diagnostics
use mpas_derived_types, only : MPAS_pool_type, MPAS_clock_type, MPAS_LOG_ERR, MPAS_LOG_CRIT
use mpas_kind_types, only : RKIND, R8KIND
use mpas_log, only : mpas_log_write
use calculate_ctt

type (MPAS_pool_type), pointer :: mesh
type (MPAS_pool_type), pointer :: state
Expand Down Expand Up @@ -240,6 +241,8 @@ subroutine convective_diagnostics_compute()
real (kind=RKIND), dimension(:), pointer :: cin
real (kind=RKIND), dimension(:), pointer :: lcl
real (kind=RKIND), dimension(:), pointer :: lfc
real (kind=RKIND), dimension(:), pointer :: ctt
real (kind=RKIND), dimension(:), pointer :: ctpres
real (kind=RKIND), dimension(:), pointer :: srh_0_1km
real (kind=RKIND), dimension(:), pointer :: srh_0_3km
real (kind=RKIND), dimension(:), pointer :: uzonal_surface
Expand All @@ -262,7 +265,7 @@ subroutine convective_diagnostics_compute()
real (kind=RKIND), dimension(:,:), pointer :: pressure_p
real (kind=RKIND), dimension(:,:), pointer :: pressure_base
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
integer, pointer :: index_qv
integer, pointer :: index_qv, index_qc, index_qi

real (kind=RKIND), parameter :: dev_motion = 7.5, z_bunker_bot = 0., z_bunker_top = 6000.
real (kind=RKIND) :: u_storm, v_storm, u_srh_bot, v_srh_bot, u_srh_top, v_srh_top
Expand All @@ -271,13 +274,17 @@ subroutine convective_diagnostics_compute()
real (kind=RKIND), dimension(:), allocatable :: dudz, dvdz, zp
real (kind=RKIND), dimension(:), allocatable :: zrel, srh
real (kind=RKIND), dimension(:), allocatable :: p_in, t_in, td_in
real (kind=8), dimension(:), allocatable :: p_in_dble
real (kind=RKIND), dimension(:), allocatable :: qc_in, qvapor_in, qi_in
real (kind=RKIND) :: ctt_out, ctpres_out

real (kind=RKIND), dimension(:,:), allocatable :: temperature, dewpoint

real (kind=RKIND) :: evp

logical :: need_cape, need_cin, need_lcl, need_lfc, need_srh_01km, need_srh_03km, need_uzonal_sfc, need_uzonal_1km, &
need_uzonal_6km, need_umerid_sfc, need_umerid_1km, need_umerid_6km, need_tsfc, need_tdsfc
need_uzonal_6km, need_umerid_sfc, need_umerid_1km, need_umerid_6km, need_tsfc, need_tdsfc, need_ctt, &
need_ctpres
logical :: need_any_diags

need_any_diags = .false.
Expand Down Expand Up @@ -310,6 +317,10 @@ subroutine convective_diagnostics_compute()
need_any_diags = need_any_diags .or. need_tsfc
need_tdsfc = MPAS_field_will_be_written('dewpoint_surface')
need_any_diags = need_any_diags .or. need_tdsfc
need_ctt = MPAS_field_will_be_written('ctt')
need_any_diags = need_any_diags .or. need_ctt
need_ctpres = MPAS_field_will_be_written('ctpres')
need_any_diags = need_any_diags .or. need_ctpres

if (need_any_diags) then

Expand All @@ -322,6 +333,8 @@ subroutine convective_diagnostics_compute()
call mpas_pool_get_array(diag, 'cin', cin)
call mpas_pool_get_array(diag, 'lcl', lcl)
call mpas_pool_get_array(diag, 'lfc', lfc)
call mpas_pool_get_array(diag, 'ctt', ctt)
call mpas_pool_get_array(diag, 'ctpres', ctpres)
call mpas_pool_get_array(diag, 'srh_0_1km', srh_0_1km)
call mpas_pool_get_array(diag, 'srh_0_3km', srh_0_3km)
call mpas_pool_get_array(diag, 'uzonal_surface', uzonal_surface)
Expand All @@ -344,6 +357,8 @@ subroutine convective_diagnostics_compute()
call mpas_pool_get_array(diag, 'pressure_p', pressure_p)

call mpas_pool_get_dimension(state, 'index_qv', index_qv)
call mpas_pool_get_dimension(state, 'index_qc', index_qc)
call mpas_pool_get_dimension(state, 'index_qi', index_qi)

allocate(temperature(nVertLevels,nCells))
allocate(dewpoint(nVertLevels,nCells))
Expand All @@ -358,6 +373,11 @@ subroutine convective_diagnostics_compute()
allocate(t_in(nVertLevels))
allocate(td_in(nVertLevels))

allocate(qc_in(nVertLevels))
allocate(qi_in(nVertLevels))
allocate(qvapor_in(nVertLevels))
allocate(p_in_dble(nVertLevels))

do iCell = 1,nCells
do k = 1,nVertLevels
temperature(k,iCell) = (theta_m(k,iCell)/(1._RKIND+rvord*scalars(index_qv,k,iCell)))*exner(k,iCell)
Expand Down Expand Up @@ -397,6 +417,25 @@ subroutine convective_diagnostics_compute()
umeridional_6km(iCell) = column_height_value(umeridional(1:nVertLevels,iCell), zp, 6000.0_RKIND, nVertLevels)
end if

! calculate ctt or HCF
if (need_ctt .OR. need_ctpres) then
! from NCL https://github.com/NCAR/wrf-python/blob/develop/fortran/wrf_fctt.f90
do iCell=1, nCellsSolve
p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell))
t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell)
qvapor_in(1:nVertLevels) = scalars(index_qv,1:nVertLevels,iCell)
qc_in(1:nVertLevels) = scalars(index_qc,1:nVertLevels,iCell)
qi_in(1:nVertLevels) = scalars(index_qi,1:nVertLevels,iCell)
zp(1:nVertLevels) = 0.5*(height(1:nVertLevels,iCell)+height(2:nVertlevels+1,iCell)) - height(1,iCell)

p_in_dble = p_in
call get_ctt(p_in_dble,DBLE(temperature(1:nVertLevels,iCell)), DBLE(scalars(index_qi,1:nVertLevels,iCell)), &
DBLE(scalars(index_qi,1:nVertLevels,iCell)) , ctt_out, &
1, 0, -9999.D0, 1.D0, nVertLevels, ctpres_out)
ctt(iCell) = ctt_out
ctpres(iCell) = ctpres_out/100.
endif

! storm-relative helicity
! first, calculate storm motion, using Bunkers formula for right-moving storms

Expand Down Expand Up @@ -482,7 +521,11 @@ subroutine convective_diagnostics_compute()
deallocate(p_in)
deallocate(t_in)
deallocate(td_in)
end if
deallocate(qc_in)
deallocate(qi_in)
deallocate(qvapor_in)
deallocate(p_in_dble)
end if ! any diags

end subroutine convective_diagnostics_compute

Expand Down