Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
7 changes: 7 additions & 0 deletions src/core_atmosphere/Registry.xml
Original file line number Diff line number Diff line change
Expand Up @@ -1026,8 +1026,15 @@
#endif

<!-- Begin convective diagnostics defined in diagnostics/Registry_convective.xml -->
<var name="rli"/>
<var name="rki"/>
<var name="rwi"/>
<var name="sfccape"/>
<var name="sfccin"/>
<var name="cape"/>
<var name="cin"/>
<var name="mlcape"/>
<var name="mlcin"/>
<var name="lcl"/>
<var name="lfc"/>
<var name="srh_0_1km"/>
Expand Down
5 changes: 4 additions & 1 deletion src/core_atmosphere/diagnostics/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -11,19 +11,22 @@ DIAGNOSTIC_MODULES = \
mpas_pv_diagnostics.o \
mpas_soundings.o \
mpas_vertical_interpolation.o \
mpas_stability_indices.o

mpas_isobaric_diagnostics.o: mpas_atm_diagnostics_utils.o mpas_vertical_interpolation.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 mpas_stability_indices.o

mpas_pv_diagnostics.o: mpas_atm_diagnostics_utils.o

mpas_soundings.o:

mpas_vertical_interpolation.o:

mpas_stability_indices.o:

################### Generally no need to modify below here ###################


Expand Down
27 changes: 25 additions & 2 deletions src/core_atmosphere/diagnostics/Registry_convective.xml
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,34 @@

<var_struct name="diag" time_levs="1">

<var name="rli" type="real" dimensions="nCells Time" units="K"
description="Lifted index"/>

<var name="rki" type="real" dimensions="nCells Time" units="K"
description="K index"/>

<var name="rwi" type="real" dimensions="nCells Time" units="K"
description="Showalter index"/>

<var name="sfccape" type="real" dimensions="nCells Time" units="J kg^{-1}"
description="Surface Convective available potential energy"/>

<var name="sfccin" type="real" dimensions="nCells Time" units="J kg^{-1}"
description="Surface Convective inhibition"/>

<var name="cape" type="real" dimensions="nCells Time" units="J kg^{-1}"
description="Convective available potential energy"/>
description="Most Unstable Convective available potential energy"/>

<var name="cin" type="real" dimensions="nCells Time" units="J kg^{-1}"
description="Convective inhibition"/>
description="Most Unstable Convective inhibition"/>

<var name="mlcape" type="real" dimensions="nCells Time" units="J k
g^{-1}"
description="Mixed Layer Convective available potential energy"/>

<var name="mlcin" type="real" dimensions="nCells Time" units="J k
g^{-1}"
description="Mixed Layer Convective inhibition"/>

<var name="lcl" type="real" dimensions="nCells Time" units="m"
description="Lifted condensation level"/>
Expand Down
3 changes: 3 additions & 0 deletions src/core_atmosphere/diagnostics/Registry_isobaric.xml
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@
<var name="temperature_isobaric" type="real" dimensions="t_iso_levels nCells Time" units="K"
description="Temperature interpolated to isobaric surfaces defined in t_iso_levels"/>

<var name="theta_isobaric" type="real" dimensions="t_iso_levels nCells Time" units="K"
description="Theta interpolated to isobaric surfaces defined in t_iso_levels"/>

<var name="dewpoint_isobaric" type="real" dimensions="t_iso_levels nCells Time" units="K"
description="Dewpoint temperature interpolated to isobaric surfaces defined in t_iso_levels"/>

Expand Down
98 changes: 93 additions & 5 deletions src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F
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 stability_indices

type (MPAS_pool_type), pointer :: mesh
type (MPAS_pool_type), pointer :: state
Expand Down Expand Up @@ -236,8 +237,12 @@ subroutine convective_diagnostics_compute()
integer, pointer :: nCells, nCellsSolve, nVertLevels, nVertLevelsP1

! Fields that are computed in this routine
real (kind=RKIND), dimension(:), pointer :: sfccape
real (kind=RKIND), dimension(:), pointer :: sfccin
real (kind=RKIND), dimension(:), pointer :: cape
real (kind=RKIND), dimension(:), pointer :: cin
real (kind=RKIND), dimension(:), pointer :: mlcape
real (kind=RKIND), dimension(:), pointer :: mlcin
real (kind=RKIND), dimension(:), pointer :: lcl
real (kind=RKIND), dimension(:), pointer :: lfc
real (kind=RKIND), dimension(:), pointer :: srh_0_1km
Expand All @@ -251,6 +256,10 @@ subroutine convective_diagnostics_compute()
real (kind=RKIND), dimension(:), pointer :: temperature_surface
real (kind=RKIND), dimension(:), pointer :: dewpoint_surface

real (kind=RKIND), dimension(:), pointer :: rli
real (kind=RKIND), dimension(:), pointer :: rki
real (kind=RKIND), dimension(:), pointer :: rwi

! Other fields used in the computation of convective diagnostics
! defined above
real (kind=RKIND), dimension(:,:), pointer :: height
Expand All @@ -268,6 +277,9 @@ subroutine convective_diagnostics_compute()
real (kind=RKIND) :: u_storm, v_storm, u_srh_bot, v_srh_bot, u_srh_top, v_srh_top
real (kind=RKIND) :: u_mean, v_mean, u_shear, v_shear, shear_magnitude
real (kind=RKIND) :: b_term, cape_out, cin_out

real (kind=RKIND) :: rli_out, rki_out, rwi_out

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
Expand All @@ -276,16 +288,34 @@ subroutine convective_diagnostics_compute()

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, &
logical :: need_sfccape, need_sfccin
logical :: need_cape, need_cin, need_mlcape, need_mlcin, 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
logical :: need_any_diags
logical :: need_rli, need_rki, need_rwi

need_any_diags = .false.

need_rli = MPAS_field_will_be_written('rli')
need_any_diags = need_any_diags .or. need_rli
need_rki = MPAS_field_will_be_written('rki')
need_any_diags = need_any_diags .or. need_rki
need_rwi = MPAS_field_will_be_written('rwi')
need_any_diags = need_any_diags .or. need_rwi

need_sfccape = MPAS_field_will_be_written('sfccape')
need_any_diags = need_any_diags .or. need_sfccape
need_sfccin = MPAS_field_will_be_written('sfccin')
need_any_diags = need_any_diags .or. need_sfccin

need_cape = MPAS_field_will_be_written('cape')
need_any_diags = need_any_diags .or. need_cape
need_cin = MPAS_field_will_be_written('cin')
need_any_diags = need_any_diags .or. need_cin
need_mlcape = MPAS_field_will_be_written('mlcape')
need_any_diags = need_any_diags .or. need_mlcape
need_mlcin = MPAS_field_will_be_written('mlcin')
need_any_diags = need_any_diags .or. need_mlcin
need_lcl = MPAS_field_will_be_written('lcl')
need_any_diags = need_any_diags .or. need_lcl
need_lfc = MPAS_field_will_be_written('lfc')
Expand Down Expand Up @@ -318,8 +348,13 @@ subroutine convective_diagnostics_compute()
call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(mesh, 'nVertLevelsP1', nVertLevelsP1)

call mpas_pool_get_array(diag, 'sfccape', sfccape)
call mpas_pool_get_array(diag, 'sfccin', sfccin)

call mpas_pool_get_array(diag, 'cape', cape)
call mpas_pool_get_array(diag, 'cin', cin)
call mpas_pool_get_array(diag, 'mlcape', mlcape)
call mpas_pool_get_array(diag, 'mlcin', mlcin)
call mpas_pool_get_array(diag, 'lcl', lcl)
call mpas_pool_get_array(diag, 'lfc', lfc)
call mpas_pool_get_array(diag, 'srh_0_1km', srh_0_1km)
Expand All @@ -333,6 +368,10 @@ subroutine convective_diagnostics_compute()
call mpas_pool_get_array(diag, 'temperature_surface', temperature_surface)
call mpas_pool_get_array(diag, 'dewpoint_surface', dewpoint_surface)

call mpas_pool_get_array(diag, 'rli', rli)
call mpas_pool_get_array(diag, 'rki', rki)
call mpas_pool_get_array(diag, 'rwi', rwi)

call mpas_pool_get_array(mesh, 'zgrid', height)
call mpas_pool_get_array(diag, 'uReconstructMeridional', umeridional)
call mpas_pool_get_array(diag, 'uReconstructZonal', uzonal)
Expand Down Expand Up @@ -451,6 +490,21 @@ subroutine convective_diagnostics_compute()
end do

! calculate cape and cin

if (need_sfccape .or. need_sfccin) then
do iCell=1, nCellsSolve
p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND
t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) - 273.15_RKIND
td_in(1:nVertLevels) = dewpoint(1:nVertLevels,iCell) - 273.15_RKIND

call getcape( 1, nVertLevels, p_in, t_in, td_in, cape_out, cin_out )

sfccape(iCell) = cape_out
sfccin(iCell) = cin_out

end do
end if

if (need_cape .or. need_cin) then
do iCell=1, nCellsSolve
p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND
Expand All @@ -463,14 +517,48 @@ subroutine convective_diagnostics_compute()
! /(17.625-log(relhum(k,iCell))-((17.625*t_in(k))/(243.04+t_in(k))))
! end do

call getcape( nVertLevels, p_in, t_in, td_in, cape_out, cin_out )
call getcape( 2, nVertLevels, p_in, t_in, td_in, cape_out, cin_out )

cape(iCell) = cape_out
cin(iCell) = cin_out

end do
end if

if (need_mlcape .or. need_mlcin) then
do iCell=1, nCellsSolve
p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND
t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) - 273.15_RKIND
td_in(1:nVertLevels) = dewpoint(1:nVertLevels,iCell) - 273.15_RKIND

call getcape( 3, nVertLevels, p_in, t_in, td_in, cape_out, cin_out )

mlcape(iCell) = cape_out
mlcin(iCell) = cin_out

end do
end if

! Stability indices

if (need_rli .or. need_rki .or. need_rwi ) then
do iCell=1, nCellsSolve
p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND
t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) - 273.15_RKIND
td_in(1:nVertLevels) = dewpoint(1:nVertLevels,iCell) - 273.15_RKIND

! print*,'Pressure (1:nVertLevels) ',p_in(1),p_in(nVertLevels)
! Enter variables from surface to top

call INDX1 (nVertLevels, p_in, t_in, td_in, rli_out, rki_out, rwi_out )

rli(iCell) = rli_out
rki(iCell) = rki_out
rwi(iCell) = rwi_out

end do
end if

deallocate(temperature)
deallocate(dewpoint)

Expand Down Expand Up @@ -592,10 +680,10 @@ end function integral_zpoint
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------

subroutine getcape( nk , p_in , t_in , td_in, cape , cin )
subroutine getcape( source, nk , p_in , t_in , td_in, cape , cin )
implicit none

integer, intent(in) :: nk
integer, intent(in) :: source, nk
real (kind=RKIND), dimension(nk), intent(in) :: p_in,t_in,td_in
real (kind=RKIND), intent(out) :: cape,cin

Expand Down Expand Up @@ -639,7 +727,7 @@ subroutine getcape( nk , p_in , t_in , td_in, cape , cin )
! results,larger number makes code
! go faster)

integer, parameter :: source = 2 ! Source parcel:
! integer, parameter :: source = 2 ! Source parcel:
! 1 = surface
! 2 = most unstable (max theta-e)
! 3 = mixed-layer (specify ml_depth)
Expand Down
22 changes: 20 additions & 2 deletions src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ module mpas_isobaric_diagnostics
need_kzm_isobaric, need_ni_isobaric, need_nr_isobaric, &
need_qr_isobaric, need_qs_isobaric, need_qv_isobaric, &
need_qc_isobaric, need_qg_isobaric, need_qi_isobaric, &
need_zgeo_isobaric
need_zgeo_isobaric, need_theta_isobaric


contains
Expand Down Expand Up @@ -152,6 +152,7 @@ subroutine isobaric_diagnostics_compute()
need_zgeo_isobaric = .false.
need_cldfrac_isobaric = .false.
need_temperature_isobaric = .false.
need_theta_isobaric = .false.
need_dewpoint_isobaric = .false.
need_relhum_isobaric = .false.
need_w_isobaric = .false.
Expand Down Expand Up @@ -198,11 +199,12 @@ subroutine isobaric_diagnostics_compute()
need_zgeo_isobaric = MPAS_field_will_be_written('zgeo_isobaric')
need_any_diags = need_any_diags .or. need_zgeo_isobaric


need_cldfrac_isobaric = MPAS_field_will_be_written('cldfrac_isobaric')
need_any_diags = need_any_diags .or. need_cldfrac_isobaric
need_temperature_isobaric = MPAS_field_will_be_written('temperature_isobaric')
need_any_diags = need_any_diags .or. need_temperature_isobaric
need_theta_isobaric = MPAS_field_will_be_written('theta_isobaric')
need_any_diags = need_any_diags .or. need_theta_isobaric
need_dewpoint_isobaric = MPAS_field_will_be_written('dewpoint_isobaric')
need_any_diags = need_any_diags .or. need_dewpoint_isobaric
need_relhum_isobaric = MPAS_field_will_be_written('relhum_isobaric')
Expand Down Expand Up @@ -325,6 +327,7 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag, diag_physics, tend_ph

real (kind=RKIND), dimension(:,:), pointer :: cldfrac_isobaric
real (kind=RKIND), dimension(:,:), pointer :: temperature_isobaric
real (kind=RKIND), dimension(:,:), pointer :: theta_isobaric
real (kind=RKIND), dimension(:,:), pointer :: dewpoint_isobaric
real (kind=RKIND), dimension(:,:), pointer :: relhum_isobaric
real (kind=RKIND), dimension(:,:), pointer :: w_isobaric
Expand Down Expand Up @@ -450,6 +453,7 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag, diag_physics, tend_ph

call mpas_pool_get_array(diag, 'cldfrac_isobaric', cldfrac_isobaric)
call mpas_pool_get_array(diag, 'temperature_isobaric', temperature_isobaric)
call mpas_pool_get_array(diag, 'theta_isobaric', theta_isobaric)
call mpas_pool_get_array(diag, 'dewpoint_isobaric', dewpoint_isobaric)
call mpas_pool_get_array(diag, 'relhum_isobaric', relhum_isobaric)
call mpas_pool_get_array(diag, 'w_isobaric', w_isobaric)
Expand Down Expand Up @@ -674,6 +678,20 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag, diag_physics, tend_ph
! call mpas_log_write('--- end interpolate temperature:')
end if

if ( need_theta_isobaric ) then

! THETA = T * (1000./P1)**Rd_Cp

do k = 1, nIsoLevelsT
do icell = 1, nCells
theta_isobaric(k,icell) = &
temperature_isobaric(k,icell)*(100000.0/t_iso_levels(k))**0.2857142
enddo
enddo

! call mpas_log_write('--- end interpolate theta:')
end if


if (need_mp_thompson .eqv. .true. .or. need_bl_mynn .eqv. .true.) then
if (NEED_NI_ISOBARIC) then
Expand Down
Loading