Skip to content

Add prate_max and other related precip variables to MPAS output #151

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

Open
Jeff-Duda wants to merge 13 commits into
base: gsl/develop
Choose a base branch
from
+121 −12
Open

Add prate_max and other related precip variables to MPAS output #151

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
958210a
Add 5 diagnostics variables to output. All are time-maxima (since las…
Jeff-Duda Jul 14, 2025
c8e3e3d
Fix to correct 1-/5-/10-minute average prate calculation.
Jeff-Duda Jul 18, 2025
6a09de4
Merge branch 'gsl/develop' into prate_max
Jeff-Duda Jul 18, 2025
2d1b7d4
Fix error in character string length in warning message.
Jeff-Duda Jul 25, 2025
fefdcba
Merge branch 'ufs-community:gsl/develop' into prate_max
Jeff-Duda Jul 30, 2025
976a44a
Merge branch 'prate_max' of github.com:Jeff-Duda/MPAS-Model into prat…
Jeff-Duda Jul 30, 2025
da0243e
Merge branch 'gsl/develop' into prate_max
clark-evans Sep 30, 2025
b2189ad
Requested changes to code to contract loopaple code.
Jeff-Duda Sep 30, 2025
1b132f3
Merge branch 'prate_max' of github.com:Jeff-Duda/MPAS-Model into prat…
Jeff-Duda Sep 30, 2025
cb5a1a6
Add IF statements to control when maximum precipitation rate calculat…
Jeff-Duda Nov 6, 2025
def3905
Corrected bug for sub-hourly max prate calculation
Jeff-Duda Dec 16, 2025
354d8c7
max_prate code cleanup. Includes removal of pseudo-SPPT temporarily a…
Jeff-Duda Feb 24, 2026
e1a7494
Update to prepare for PR of prate_max
Jeff-Duda Mar 13, 2026
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
27 changes: 27 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 @@ -682,6 +682,7 @@
<var name="max_hail_diameter_column"/>
<var name="max_hail_diameter_column_acc"/>
<var name="refl10cm_max"/>
<var name="refl10cm_compref_max"/>
Copy link
Collaborator

@AndersJensen-NOAA AndersJensen-NOAA Feb 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fix spacing

<var name="i_rainnc"/>
<var name="rainncv"/>
<var name="rainnc"/>
Expand Down Expand Up @@ -1142,13 +1143,15 @@
<var name="swupt"/>
<var name="rainc"/>
<var name="rainnc"/>
<var name="prate_max"/>
Copy link
Collaborator

@AndersJensen-NOAA AndersJensen-NOAA Feb 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fix spacing

<var name="vcpool"/>
<var name="refl10cm"/>
<var name="max_hail_diameter_sfc"/>
<var name="max_hail_diameter_sfc_acc"/>
<var name="max_hail_diameter_column"/>
<var name="max_hail_diameter_column_acc"/>
<var name="refl10cm_max"/>
<var name="refl10cm_compref_max"/>
Copy link
Collaborator

@AndersJensen-NOAA AndersJensen-NOAA Feb 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fix spacing

<var name="refl10cm_1km"/>
<var name="refl10cm_1km_max"/>
<var name="precipw"/>
Expand Down Expand Up @@ -2981,6 +2984,10 @@
units="s"
description="Relaxation time to initial values in 1-d OML"
possible_values="Non-negative real values"/>

<nml_option name="calc_prate_max" type="logical" default_value="false" in_defaults="false" units="-"
description="switch for calculating time-maximum precipitation rate for instantaneous, 1-minute, 5-minute, and 10-minute periods"
possible_values=".true. or .false."/>
</nml_record>

<var_struct name="diag_physics" time_levs="1">
Expand Down Expand Up @@ -3036,6 +3043,10 @@
description="10 cm maximum radar reflectivity"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="refl10cm_compref_max" type="real" dimensions="nCells Time" units="dBZ"
description="maximum composite 10 cm radar reflectivity since last output time"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="refl10cm_1km" type="real" dimensions="nCells Time" units="dBZ"
description="diagnosed 10 cm radar reflectivity at 1 km AGL"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>
Expand Down Expand Up @@ -3068,6 +3079,18 @@
description="accumulated total grid-scale precipitation"
packages="mp_kessler_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="ave_prate_1min" type="real" dimensions="nCells Time" units="mm s^{-1}"
description="maximum 1-minute precipitation rate since last output time"
packages="mp_kessler_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="ave_prate_5min" type="real" dimensions="nCells Time" units="mm s^{-1}"
description="maximum 5-minute precipitation rate since last output time"
packages="mp_kessler_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="ave_prate_10min" type="real" dimensions="nCells Time" units="mm s^{-1}"
description="maximum 10-minute precipitation rate since last output time"
packages="mp_kessler_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="frainnc" type="real" dimensions="nCells Time" units="mm"
description="accumulated total grid-scale freezing rain"
packages="mp_tempo_in"/>
Expand All @@ -3080,6 +3103,10 @@
description="accumulated grid-scale precipitation of graupel"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="prate_max" type="real" dimensions="nCells Time" units="mm s^{-1}"
description="maximum instantaneous precipitation rate since last output time"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>

<var name="rainprod" type="real" dimensions="nVertLevels nCells Time" units="s^{-1}"
description="rain production rate"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/>
Expand Down
24 changes: 22 additions & 2 deletions src/core_atmosphere/mpas_atm_core.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -981,13 +981,33 @@ subroutine atm_reset_diagnostics(diag_physics)

type (mpas_pool_type), pointer :: diag_physics

real (kind=RKIND), dimension(:), pointer :: refl10cm_1km_max, max_hail_diameter_sfc_acc, &
max_hail_diameter_column_acc
real (kind=RKIND), dimension(:), pointer :: refl10cm_1km_max, max_hail_diameter_sfc_acc, prate_max, &
refl10cm_compref_max, ave_prate_1min, ave_prate_5min, ave_prate_10min, max_hail_diameter_column_acc

#ifdef DO_PHYSICS
call mpas_pool_get_array(diag_physics, 'prate_max', prate_max)
call mpas_pool_get_array(diag_physics, 'refl10cm_1km_max', refl10cm_1km_max)
call mpas_pool_get_array(diag_physics, 'refl10cm_compref_max', refl10cm_compref_max)
call mpas_pool_get_array(diag_physics, 'ave_prate_1min', ave_prate_1min)
call mpas_pool_get_array(diag_physics, 'ave_prate_5min', ave_prate_5min)
call mpas_pool_get_array(diag_physics, 'ave_prate_10min', ave_prate_10min)
call mpas_pool_get_array(diag_physics, 'max_hail_diameter_sfc_acc', max_hail_diameter_sfc_acc)
call mpas_pool_get_array(diag_physics, 'max_hail_diameter_column_acc', max_hail_diameter_column_acc)
if(associated(prate_max)) then
prate_max(:) = 0.
endif
if(associated(refl10cm_compref_max)) then
refl10cm_compref_max(:) = -35.
endif
if(associated(ave_prate_1min)) then
ave_prate_1min(:) = 0.
endif
if(associated(ave_prate_5min)) then
ave_prate_5min(:) = 0.
endif
if(associated(ave_prate_10min)) then
ave_prate_10min(:) = 0.
endif
if(associated(refl10cm_1km_max)) then
refl10cm_1km_max(:) = 0.
endif
Expand Down
77 changes: 68 additions & 9 deletions src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -779,7 +779,9 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten

!ensure that we only call compute_radar_reflectivity() if we are using an MPS that supports
!the computation of simulated radar reflectivity:
if (l_diags) then
!Compute reflectivity at all time steps for supported MP schemes rather than just at history intervals
!Only computing reflectivity occasionally leads to bad looking time-maximum reflectivity values
if (l_diags) then
if(trim(microp_scheme) == "mp_wsm6" .or. &
trim(microp_scheme) == "mp_thompson" .or. &
trim(microp_scheme) == "mp_thompson_aerosols" .or. &
Expand All @@ -801,7 +803,6 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
if (do_diag_dbz_flag) then
if(trim(microp_scheme) == "mp_tempo" .or. &
trim(microp_scheme) == "mp_nssl2m" ) then
call mpas_log_write('Computing hourly max reflectivity')
call compute_hourly_max_radar_reflectivity(configs,diag_physics,its,ite)
else
call mpas_log_write('*** NOTICE: NOT computing hourly max simulated radar reflectivity')
Expand Down Expand Up @@ -951,29 +952,72 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite)

real(kind=RKIND),pointer:: config_bucket_rainnc
real(kind=RKIND),dimension(:),pointer:: precipw
real(kind=RKIND),dimension(:),pointer:: prate_max
real(kind=RKIND),dimension(:),pointer:: graupelnc,rainnc,snownc
real(kind=RKIND),dimension(:),pointer:: graupelncv,rainncv,snowncv,sr
real(kind=RKIND),dimension(:),pointer:: frainnc
real(kind=RKIND),dimension(:),pointer:: ave_prate_1min,ave_prate_5min,ave_prate_10min
real(kind=RKIND),dimension(:,:),allocatable,save:: rolling_precip

!local variables and arrays:
integer:: i,j,k
character(len=16) :: desc
integer:: i,j,k,z
integer, parameter :: n_sub_windows = 3
integer :: longest_window_n
integer, dimension(n_sub_windows) :: sub_limit, sub_window_len
logical, save :: alloc_rolling_precip_array = .false.
logical, pointer :: calc_prate_max
real(kind=RKIND), dimension(n_sub_windows) :: n_dt_sub, sub_window_time
real(kind=RKIND):: rho_a

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

call mpas_pool_get_config(configs,'config_microp_scheme',microp_scheme )
call mpas_pool_get_config(configs,'config_nssl_moments',nssl_moments)
call mpas_pool_get_config(configs,'config_bucket_rainnc',config_bucket_rainnc)
call mpas_pool_get_config(configs,'calc_prate_max',calc_prate_max)

call mpas_pool_get_array(diag_physics,'i_rainnc' ,i_rainnc )
call mpas_pool_get_array(diag_physics,'precipw' ,precipw )
call mpas_pool_get_array(diag_physics,'graupelnc' ,graupelnc )
call mpas_pool_get_array(diag_physics,'graupelncv',graupelncv)
call mpas_pool_get_array(diag_physics,'rainnc' ,rainnc )
call mpas_pool_get_array(diag_physics,'prate_max' ,prate_max )
call mpas_pool_get_array(diag_physics,'rainncv' ,rainncv )
call mpas_pool_get_array(diag_physics,'snownc' ,snownc )
call mpas_pool_get_array(diag_physics,'snowncv' ,snowncv )
call mpas_pool_get_array(diag_physics,'sr' ,sr )
call mpas_pool_get_array(diag_physics,'ave_prate_1min',ave_prate_1min)
call mpas_pool_get_array(diag_physics,'ave_prate_5min',ave_prate_5min)
call mpas_pool_get_array(diag_physics,'ave_prate_10min',ave_prate_10min)

if (calc_prate_max) THEN
sub_window_len(:n_sub_windows) = 60*(/1,5,10/) ! min -> sec
sub_limit = -1
do z = 1,n_sub_windows
n_dt_sub(z) = sub_window_len(z)/dt_microp
! Account for time steps that do not divide evenly into 1, 5, or 10 minutes
if (abs(fraction(n_dt_sub(z))) .gt. 1e-6) then
sub_limit(z) = nint(n_dt_sub(z)) + 1
sub_window_time(z) = dt_microp*(sub_limit(z)-1)
else
sub_limit(z) = sub_window_len(z)/dt_microp + 1
sub_window_time(z) = sub_window_len(z)
end if
if (sub_limit(z) == 1) then
write(desc,'(I3)') int(sub_window_len(z)/60.)
call mpas_log_write('WARNING: model time step is coarser than '//trim(adjustl(desc))//' minute(s) and thus the ave_prate_[1-/5-/10-]min array will be either all zeros or meaningless.')
alloc_rolling_precip_array = .true.
end if
end do
! if all is okay, allocate rolling precip array once and then never again
longest_window_n = maxval(sub_limit)
if (.not. alloc_rolling_precip_array) then
allocate(rolling_precip(longest_window_n,its:ite))
rolling_precip = 0._RKIND
alloc_rolling_precip_array = .true.
end if
END IF ! calc_prate_max

do i = its,ite
precipw(i) = 0._RKIND
Expand All @@ -991,15 +1035,30 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite)

!time-step precipitation:
rainncv(i) = rainnc_p(i,j)
IF (calc_prate_max) prate_max(i) = max(rainnc_p(i,j)/dt_microp,prate_max(i))

!accumulated precipitation:
rainnc(i) = rainnc(i) + rainncv(i)

if(l_acrain .and. config_bucket_rainnc.gt.0._RKIND .and. &
rainnc(i).gt.config_bucket_rainnc) then
i_rainnc(i) = i_rainnc(i) + 1
rainnc(i) = rainnc(i) - config_bucket_rainnc
endif
if(l_acrain .and. config_bucket_rainnc.gt.0._RKIND .and. &
rainnc(i).gt.config_bucket_rainnc) then
i_rainnc(i) = i_rainnc(i) + 1
rainnc(i) = rainnc(i) - config_bucket_rainnc
endif

IF (calc_prate_max) THEN
! update rolling precipitation totals for sub-history-interval average-precipitation-rate calculation
! 1-minute prate
IF (allocated(rolling_precip)) THEN
DO k = longest_window_n-1,1,-1
rolling_precip(k+1,i) = rolling_precip(k,i)
END DO
rolling_precip(1,i) = rainnc(i)
ave_prate_1min(i) = max((rolling_precip(1,i) - rolling_precip(sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i))
Copy link
Collaborator

@dustinswales dustinswales Feb 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This approach taken here always assumes 3 averaging intervals, which are always set to 1, 5 and 10 minutes.
What if someone wants to have more, or different, averaging intervals?

ave_prate_5min(i) = max((rolling_precip(1,i) - rolling_precip(sub_limit(2),i))/sub_window_time(2),ave_prate_5min(i))
ave_prate_10min(i) = max((rolling_precip(1,i) - rolling_precip(sub_limit(3),i))/sub_window_time(3),ave_prate_10min(i))
END IF
END IF ! calc_prate_max

enddo
enddo
Expand Down Expand Up @@ -1293,7 +1352,7 @@ subroutine compute_hourly_max_radar_reflectivity(configs,diag_physics,its,ite)
call mpas_pool_get_config(configs,'config_microp_scheme',microp_scheme)
call mpas_pool_get_config(configs,'config_nssl_moments',nssl_moments) ! not needed?

! call mpas_pool_get_array(diag_physics,'refl10cm_max',refl10cm_max)
!call mpas_pool_get_array(diag_physics,'refl10cm_max',refl10cm_max)
call mpas_pool_get_array(diag_physics,'refl10cm_1km',refl10cm_1km)
call mpas_pool_get_array(diag_physics,'refl10cm_1km_max',refl10cm_1km_max)

Expand Down
5 changes: 4 additions & 1 deletion src/core_atmosphere/physics/mpas_atmphys_interface.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1171,7 +1171,7 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
real(kind=RKIND),dimension(:,:),pointer :: nc,ni,nr,nifa,nwfa
real(kind=RKIND),dimension(:,:),pointer :: ns,ng,nh,nccn
real(kind=RKIND),dimension(:,:),pointer :: volg,volh,zrw,zgw,zhw
real(kind=RKIND),dimension(:) ,pointer :: max_hail_diameter_sfc, max_hail_diameter_column
real(kind=RKIND),dimension(:) ,pointer :: max_hail_diameter_sfc, max_hail_diameter_column, refl10cm_compref_max, compref
real(kind=RKIND),dimension(:,:),pointer :: rainprod,evapprod,refl10cm
real(kind=RKIND),dimension(:,:),pointer :: re_cloud,re_ice,re_snow
real(kind=RKIND),dimension(:,:),pointer :: rthmpten,rqvmpten,rqcmpten,rqrmpten,rqimpten,rqsmpten,rqgmpten
Expand Down Expand Up @@ -1400,8 +1400,11 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te

call mpas_pool_get_array(diag_physics,'max_hail_diameter_sfc' ,max_hail_diameter_sfc)
call mpas_pool_get_array(diag_physics,'max_hail_diameter_column' ,max_hail_diameter_column)
call mpas_pool_get_array(diag_physics,'refl10cm_compref_max',refl10cm_compref_max )
call mpas_pool_get_array(diag_physics,'refl10cm_max',compref )
do j = jts,jte
do i = its,ite
refl10cm_compref_max(i) = max(refl10cm_compref_max(i),compref(i))
max_hail_diameter_sfc(i) = max_hail_diameter_sfc_p(i,j)
max_hail_diameter_column(i) = max_hail_diameter_column_p(i,j)
enddo
Expand Down
Loading