From 2c128ea04c6e346c2e0d5fb2a26a573540155d5d Mon Sep 17 00:00:00 2001 From: thomasuhoh Date: Fri, 9 Jun 2023 08:53:18 +0200 Subject: [PATCH] Add cloud top temperature and pressure estimates from NCL routines to diag output --- src/core_atmosphere/Registry.xml | 2 + src/core_atmosphere/diagnostics/Makefile | 3 +- .../diagnostics/Registry_convective.xml | 6 ++ src/core_atmosphere/diagnostics/calc_ctt.F | 101 ++++++++++++++++++ .../diagnostics/mpas_convective_diagnostics.F | 49 ++++++++- 5 files changed, 157 insertions(+), 4 deletions(-) create mode 100644 src/core_atmosphere/diagnostics/calc_ctt.F diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index b3127cce3b..7f01e5e6fa 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -951,6 +951,8 @@ + + diff --git a/src/core_atmosphere/diagnostics/Makefile b/src/core_atmosphere/diagnostics/Makefile index 614bc1c137..4dc6febccf 100644 --- a/src/core_atmosphere/diagnostics/Makefile +++ b/src/core_atmosphere/diagnostics/Makefile @@ -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 diff --git a/src/core_atmosphere/diagnostics/Registry_convective.xml b/src/core_atmosphere/diagnostics/Registry_convective.xml index 59ad1c2f58..c94ce7a6fe 100644 --- a/src/core_atmosphere/diagnostics/Registry_convective.xml +++ b/src/core_atmosphere/diagnostics/Registry_convective.xml @@ -55,5 +55,11 @@ + + + + diff --git a/src/core_atmosphere/diagnostics/calc_ctt.F b/src/core_atmosphere/diagnostics/calc_ctt.F new file mode 100644 index 0000000000..18a41c4b06 --- /dev/null +++ b/src/core_atmosphere/diagnostics/calc_ctt.F @@ -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 diff --git a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F index 163ee3774f..63e017e5a3 100644 --- a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F @@ -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 @@ -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 @@ -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 @@ -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. @@ -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 @@ -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) @@ -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)) @@ -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) @@ -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 @@ -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