From 958210aa99519c2533c5ba4b22bd2fcef78aef22 Mon Sep 17 00:00:00 2001 From: Jeffrey Duda Date: Mon, 14 Jul 2025 17:16:53 +0000 Subject: [PATCH 1/6] Add 5 diagnostics variables to output. All are time-maxima (since last output time): -Instantaneous precipitation rate -Composite reflectivity -1-minute average precipitation rate -5-minute average precipitation rate -10-minute average precipitation rate --- src/core_atmosphere/Registry.xml | 23 ++++ src/core_atmosphere/mpas_atm_core.F | 23 +++- .../mpas_atmphys_driver_microphysics.F | 104 ++++++++++++++++-- .../physics/mpas_atmphys_interface.F | 9 +- .../physics/mpas_atmphys_vars.F | 6 +- 5 files changed, 149 insertions(+), 16 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index d537b63f8a..56cc169e30 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -665,6 +665,7 @@ + @@ -1125,12 +1126,14 @@ + + @@ -2759,6 +2762,10 @@ description="10 cm maximum radar reflectivity" packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/> + + @@ -2790,6 +2797,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"/> + + + + + + @@ -2802,6 +2821,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"/> + + diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index bb04c4addb..ef4bb8e76f 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -955,11 +955,32 @@ 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 + 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 #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) + 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 diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index caf053e0ff..d9835b5856 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -136,6 +136,7 @@ subroutine allocate_microphysics(configs) !surface precipitation: if(.not.allocated(rainnc_p) ) allocate(rainnc_p(ims:ime,jms:jme) ) if(.not.allocated(rainncv_p)) allocate(rainncv_p(ims:ime,jms:jme)) + if(.not.allocated(prate_p) ) allocate(prate_p(ims:ime,jms:jme) ) microp_select: select case(trim(microp_scheme)) case ("mp_thompson","mp_thompson_aerosols","mp_wsm6") @@ -212,6 +213,7 @@ subroutine allocate_microphysics(configs) if(.not.allocated(nr_p) ) allocate(nr_p(ims:ime,kms:kme,jms:jme)) if(.not.allocated(refl10cm_p) ) allocate(refl10cm_p(ims:ime,kms:kme,jms:jme) ) if(.not.allocated(refl10cm_1km_p) ) allocate(refl10cm_1km_p(ims:ime,jms:jme) ) + if(.not.allocated(refl10cm_max_p)) allocate(refl10cm_max_p(ims:ime,jms:jme) ) if(.not.allocated(frainnc_p) ) allocate(frainnc_p(ims:ime,jms:jme) ) if(.not.allocated(max_hail_diameter_sfc_p) ) allocate(max_hail_diameter_sfc_p(ims:ime,jms:jme) ) if(.not.allocated(max_hail_diameter_column_p) ) allocate(max_hail_diameter_column_p(ims:ime,jms:jme) ) @@ -301,6 +303,7 @@ subroutine deallocate_microphysics(configs) !surface precipitation: if(allocated(rainnc_p) ) deallocate(rainnc_p ) if(allocated(rainncv_p)) deallocate(rainncv_p) + if(allocated(prate_p)) deallocate(prate_p) microp_select: select case(trim(microp_scheme)) case ("mp_thompson","mp_thompson_aerosols","mp_wsm6") @@ -377,6 +380,7 @@ subroutine deallocate_microphysics(configs) if(allocated(nr_p) ) deallocate(nr_p ) if(allocated(refl10cm_p) ) deallocate(refl10cm_p ) if(allocated(refl10cm_1km_p)) deallocate(refl10cm_1km_p) + if(allocated(refl10cm_max_p)) deallocate(refl10cm_max_p) if(allocated(frainnc_p) ) deallocate(frainnc_p) if(allocated(max_hail_diameter_sfc_p) ) deallocate(max_hail_diameter_sfc_p) if(allocated(max_hail_diameter_column_p) ) deallocate(max_hail_diameter_column_p) @@ -837,7 +841,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. & @@ -854,12 +860,12 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !over ice otherwise. call compute_relhum(diag,its,ite) - end if +! end if 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 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') @@ -925,6 +931,7 @@ subroutine precip_from_MPAS(configs,diag_physics,its,ite) do i = its, ite rainncv_p(i,j) = 0._RKIND rainnc_p(i,j) = 0._RKIND + prate_p(i,j) = 0._RKIND enddo enddo @@ -1009,12 +1016,20 @@ 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_1min,rolling_precip_5min,rolling_precip_10min !local variables and arrays: - integer:: i,j,k + character(len=3) :: desc + integer:: i,j,k,z + integer :: n_sub_windows + integer, dimension(5) :: sub_limit, sub_window_len + logical, dimension(5), save :: init_rolling_precip_array = .false. + real(kind=RKIND), dimension(5) :: n_dt_sub, sub_window_time real(kind=RKIND):: rho_a !----------------------------------------------------------------------------------------------------------------- @@ -1028,10 +1043,49 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) 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) + + n_sub_windows = 3 + sub_window_len(:n_sub_windows) = 60*(/1,5,10/) ! min -> sec + 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))) .le. 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,*) sub_window_len(z)/60. + call mpas_log_write('WARNING: model time step is coarser than '//desc//' minute(s) and thus the ave_prate_[]min array will be either all zeros or meaningless!') + init_rolling_precip_array(z) = .true. + end if + end do + ! if all is okay, allocate rolling precip arrays just once and never again + if (.not. init_rolling_precip_array(1)) then + allocate(rolling_precip_1min(sub_limit(1),its:ite)) + rolling_precip_1min = 0._RKIND + init_rolling_precip_array(1) = .true. + end if + if (.not. init_rolling_precip_array(2)) then + allocate(rolling_precip_5min(sub_limit(2),its:ite)) + rolling_precip_5min = 0._RKIND + init_rolling_precip_array(2) = .true. + end if + if (.not. init_rolling_precip_array(3)) then + allocate(rolling_precip_10min(sub_limit(3),its:ite)) + rolling_precip_10min = 0._RKIND + init_rolling_precip_array(3) = .true. + end if do i = its,ite precipw(i) = 0._RKIND @@ -1049,15 +1103,42 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) !time-step precipitation: rainncv(i) = rainnc_p(i,j) + 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 + ! update rolling precipitation totals for sub-history-interval average-precipitation-rate calculation + ! 1-minute prate + IF (allocated(rolling_precip_1min)) THEN + DO k = 1,sub_limit(1)-1 + rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) + END DO + rolling_precip_1min(1,i) = rainnc_p(i,j) + ave_prate_1min(i) = max(ave_prate_1min(i),(rolling_precip_1min(sub_limit(1),i) - rolling_precip_1min(1,i))/sub_window_time(1)) + END IF + ! 5-minute prate + IF (allocated(rolling_precip_5min)) THEN + DO k = 1,sub_limit(2)-1 + rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:) + END DO + rolling_precip_5min(1,i) = rainnc_p(i,j) + ave_prate_5min(i) = max(ave_prate_5min(i),(rolling_precip_5min(sub_limit(2),i) - rolling_precip_5min(1,i))/sub_window_time(2)) + END IF + ! 10-minute prate + IF (allocated(rolling_precip_10min)) THEN + DO k = 1,sub_limit(3)-1 + rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:) + END DO + rolling_precip_10min(1,i) = rainnc_p(i,j) + ave_prate_10min(i) = max(ave_prate_10min(i),(rolling_precip_10min(sub_limit(3),i) - rolling_precip_10min(1,i))/sub_window_time(3)) + END IF + + 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 enddo enddo @@ -1172,6 +1253,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite) enddo refl10cm_max(i) = maxval(dBZ1d(:)) + refl10cm_max_p(i,j) = refl10cm_max(i) w1 = (zp(kp+1)-1000.)/(zp(kp+1)-zp(kp)) w2 = 1.0 - w1 refl10cm_1km(i) = w1*dBZ1d(kp) + w2*dBZ1d(kp+1) @@ -1205,6 +1287,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite) enddo refl10cm_max(i) = maxval(dBZ1d(:)) + refl10cm_max_p(i,j) = refl10cm_max(i) w1 = (zp(kp+1)-1000.)/(zp(kp+1)-zp(kp)) w2 = 1.0 - w1 refl10cm_1km(i) = w1*dBZ1d(kp) + w2*dBZ1d(kp+1) @@ -1349,7 +1432,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) @@ -1395,7 +1478,6 @@ subroutine compute_hourly_max_radar_reflectivity(configs,diag_physics,its,ite) do i = its,ite refl10cm_1km_max(i) = max(refl10cm_1km_max(i),refl10cm_1km(i)) enddo - case default end select microp_select diff --git a/src/core_atmosphere/physics/mpas_atmphys_interface.F b/src/core_atmosphere/physics/mpas_atmphys_interface.F index 63ceeece46..314225575a 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_interface.F +++ b/src/core_atmosphere/physics/mpas_atmphys_interface.F @@ -680,7 +680,7 @@ subroutine microphysics_from_MPAS(configs,mesh,state,time_lev,diag,diag_physics, real(kind=RKIND),dimension(:,:),pointer :: volg,volh,zrw,zgw,zhw real(kind=RKIND),dimension(:,:),pointer :: rainprod,evapprod real(kind=RKIND),dimension(:,:),pointer :: refl10cm - 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 real(kind=RKIND),dimension(:,:),pointer :: re_cloud,re_ice,re_snow real(kind=RKIND),dimension(:,:),pointer :: qc_bl, cldfrac_bl real(kind=RKIND),dimension(:,:),pointer :: rthmpten,rqvmpten,rqcmpten,rqrmpten,rqimpten,rqsmpten,rqgmpten @@ -863,10 +863,12 @@ subroutine microphysics_from_MPAS(configs,mesh,state,time_lev,diag,diag_physics, call mpas_pool_get_dimension(state,'index_nr',index_nr) 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) ni => scalars(index_ni,:,:) nr => scalars(index_nr,:,:) do j = jts,jte do i = its,ite + refl10cm_max_p(i,j) = refl10cm_compref_max(i) max_hail_diameter_sfc_p(i,j) = max_hail_diameter_sfc(i) max_hail_diameter_column_p(i,j) = max_hail_diameter_column(i) enddo @@ -1173,7 +1175,7 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te real(kind=RKIND),dimension(:,:),pointer :: volg,volh,zrw,zgw,zhw real(kind=RKIND),dimension(:,:),pointer :: rainprod,evapprod real(kind=RKIND),dimension(:,:),pointer :: refl10cm - 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 :: re_cloud,re_ice,re_snow real(kind=RKIND),dimension(:,:),pointer :: rthmpten,rqvmpten,rqcmpten,rqrmpten,rqimpten,rqsmpten,rqgmpten real(kind=RKIND),dimension(:,:),pointer :: rncmpten,rnimpten,rnrmpten,rnifampten,rnwfampten @@ -1400,8 +1402,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 diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F index 7412b57438..b3de04ffa9 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F @@ -284,7 +284,8 @@ module mpas_atmphys_vars graupelnc_p, &! graupelncv_p, &! sr_p, &! - nwfa2d_p + nwfa2d_p, & + prate_p real(kind=RKIND),dimension(:,:),allocatable:: & frainnc_p @@ -303,7 +304,8 @@ module mpas_atmphys_vars refl10cm_p, &! refl10cm_cu_p ! real(kind=RKIND),dimension(:,:),allocatable:: & - refl10cm_1km_p + refl10cm_1km_p, & + refl10cm_max_p real(kind=RKIND),dimension(:,:),allocatable:: & max_hail_diameter_sfc_p, max_hail_diameter_column_p !================================================================================================================= From c8e3e3dac2589b8268be666006c36c50be50f616 Mon Sep 17 00:00:00 2001 From: Jeffrey Duda Date: Fri, 18 Jul 2025 17:43:31 +0000 Subject: [PATCH 2/6] Fix to correct 1-/5-/10-minute average prate calculation. --- .../mpas_atmphys_driver_microphysics.F | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index d9835b5856..f01d0b0af7 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -1057,7 +1057,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) 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))) .le. 1e-6) then + 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 @@ -1111,27 +1111,27 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) ! update rolling precipitation totals for sub-history-interval average-precipitation-rate calculation ! 1-minute prate IF (allocated(rolling_precip_1min)) THEN - DO k = 1,sub_limit(1)-1 + DO k = sub_limit(1)-1,1,-1 rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) END DO - rolling_precip_1min(1,i) = rainnc_p(i,j) - ave_prate_1min(i) = max(ave_prate_1min(i),(rolling_precip_1min(sub_limit(1),i) - rolling_precip_1min(1,i))/sub_window_time(1)) + rolling_precip_1min(1,i) = rainnc(i) + ave_prate_1min(i) = max((rolling_precip_1min(1,i) - rolling_precip_1min(sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) END IF ! 5-minute prate IF (allocated(rolling_precip_5min)) THEN - DO k = 1,sub_limit(2)-1 + DO k = sub_limit(2)-1,1,-1 rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:) END DO - rolling_precip_5min(1,i) = rainnc_p(i,j) - ave_prate_5min(i) = max(ave_prate_5min(i),(rolling_precip_5min(sub_limit(2),i) - rolling_precip_5min(1,i))/sub_window_time(2)) + rolling_precip_5min(1,i) = rainnc(i) + ave_prate_5min(i) = max((rolling_precip_5min(1,i) - rolling_precip_5min(sub_limit(2),i))/sub_window_time(2),ave_prate_5min(i)) END IF ! 10-minute prate IF (allocated(rolling_precip_10min)) THEN - DO k = 1,sub_limit(3)-1 + DO k = sub_limit(3)-1,1,-1 rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:) END DO - rolling_precip_10min(1,i) = rainnc_p(i,j) - ave_prate_10min(i) = max(ave_prate_10min(i),(rolling_precip_10min(sub_limit(3),i) - rolling_precip_10min(1,i))/sub_window_time(3)) + rolling_precip_10min(1,i) = rainnc(i) + ave_prate_10min(i) = max((rolling_precip_10min(1,i) - rolling_precip_10min(sub_limit(3),i))/sub_window_time(3),ave_prate_10min(i)) END IF if(l_acrain .and. config_bucket_rainnc.gt.0._RKIND .and. & @@ -1270,7 +1270,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite) ! Reflectivity is already computed in 3D in the refl10cm_p array, so just read out the values if(.not.allocated(dBz1d)) allocate(dBZ1d(kts:kte)) if(.not.allocated(zp) ) allocate(zp(kts:kte) ) - call mpas_log_write("Computing refl10cm_1km_max") + !call mpas_log_write("Computing refl10cm_1km_max") !call mpas_log_write("max before = $r dbz", realArgs=(/maxval(refl10cm_1km_max)/)) do j = jts,jte ! Assuming here that jts == jte (2D arrays) do i = its,ite @@ -1420,7 +1420,7 @@ subroutine compute_hourly_max_radar_reflectivity(configs,diag_physics,its,ite) !local pointers: character(len=StrKIND),pointer:: microp_scheme character(len=StrKIND),pointer:: nssl_moments - real(kind=RKIND),dimension(:),pointer:: refl10cm_1km_max, refl10cm_1km + real(kind=RKIND),dimension(:),pointer:: refl10cm_1km_max, refl10cm_1km, refl10cm_max !local variables and arrays: integer:: i,j,k,kp From 2d1b7d4a1c22d3b8b2894c64f0fd7acade37a08a Mon Sep 17 00:00:00 2001 From: Jeffrey Duda Date: Fri, 25 Jul 2025 22:18:18 +0000 Subject: [PATCH 3/6] Fix error in character string length in warning message. --- .../physics/mpas_atmphys_driver_microphysics.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index c233a73124..d5aa1d77a4 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -1024,7 +1024,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) real(kind=RKIND),dimension(:,:),allocatable,save:: rolling_precip_1min,rolling_precip_5min,rolling_precip_10min !local variables and arrays: - character(len=3) :: desc + character(len=16) :: desc integer:: i,j,k,z integer :: n_sub_windows integer, dimension(5) :: sub_limit, sub_window_len @@ -1065,8 +1065,8 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) sub_window_time(z) = sub_window_len(z) end if if (sub_limit(z) == 1) then - write(desc,*) sub_window_len(z)/60. - call mpas_log_write('WARNING: model time step is coarser than '//desc//' minute(s) and thus the ave_prate_[]min array will be either all zeros or meaningless!') + 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.') init_rolling_precip_array(z) = .true. end if end do From b2189ad2b9ad471a2d2c09b7b52c2e3fe9827ee6 Mon Sep 17 00:00:00 2001 From: Jeffrey Duda Date: Tue, 30 Sep 2025 20:41:02 +0000 Subject: [PATCH 4/6] Requested changes to code to contract loopaple code. --- .../mpas_atmphys_driver_microphysics.F | 81 ++++++++++--------- 1 file changed, 41 insertions(+), 40 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index d5aa1d77a4..fcbcdbf583 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -1021,14 +1021,16 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) 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_1min,rolling_precip_5min,rolling_precip_10min +! real(kind=RKIND),dimension(:,:),allocatable,save:: rolling_precip_1min,rolling_precip_5min,rolling_precip_10min + real(kind=RKIND),dimension(:,:,:),allocatable,save:: rolling_precip !local variables and arrays: character(len=16) :: desc integer:: i,j,k,z integer :: n_sub_windows integer, dimension(5) :: sub_limit, sub_window_len - logical, dimension(5), save :: init_rolling_precip_array = .false. +! logical, dimension(5), save :: init_rolling_precip_array = .false. + logical, save :: alloc_rolling_precip_array = .false. real(kind=RKIND), dimension(5) :: n_dt_sub, sub_window_time real(kind=RKIND):: rho_a @@ -1054,6 +1056,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) n_sub_windows = 3 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 @@ -1067,24 +1070,14 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) 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.') - init_rolling_precip_array(z) = .true. + alloc_rolling_precip_array = .true. end if end do - ! if all is okay, allocate rolling precip arrays just once and never again - if (.not. init_rolling_precip_array(1)) then - allocate(rolling_precip_1min(sub_limit(1),its:ite)) - rolling_precip_1min = 0._RKIND - init_rolling_precip_array(1) = .true. - end if - if (.not. init_rolling_precip_array(2)) then - allocate(rolling_precip_5min(sub_limit(2),its:ite)) - rolling_precip_5min = 0._RKIND - init_rolling_precip_array(2) = .true. - end if - if (.not. init_rolling_precip_array(3)) then - allocate(rolling_precip_10min(sub_limit(3),its:ite)) - rolling_precip_10min = 0._RKIND - init_rolling_precip_array(3) = .true. + ! if all is okay, allocate rolling precip array once and then never again + if (.not. alloc_rolling_precip_array) then + allocate(rolling_precip(maxval(/sub_limit(1),sub_limit(2),sub_limit(3)/),its:ite)) + rolling_precip = 0._RKIND + alloc_rolling_precip_array = .true. end if do i = its,ite @@ -1110,29 +1103,37 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) ! update rolling precipitation totals for sub-history-interval average-precipitation-rate calculation ! 1-minute prate - IF (allocated(rolling_precip_1min)) THEN - DO k = sub_limit(1)-1,1,-1 - rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) - END DO - rolling_precip_1min(1,i) = rainnc(i) - ave_prate_1min(i) = max((rolling_precip_1min(1,i) - rolling_precip_1min(sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) - END IF - ! 5-minute prate - IF (allocated(rolling_precip_5min)) THEN - DO k = sub_limit(2)-1,1,-1 - rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:) + IF (allocated(rolling_precip)) THEN +! IF (allocated(rolling_precip_1min)) THEN + DO z = 1,3 + DO k = sub_limit(z)-1,1,-1 +! rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) + rolling_precip(z,k+1,:) = rolling_precip(z,k,:) END DO - rolling_precip_5min(1,i) = rainnc(i) - ave_prate_5min(i) = max((rolling_precip_5min(1,i) - rolling_precip_5min(sub_limit(2),i))/sub_window_time(2),ave_prate_5min(i)) - END IF - ! 10-minute prate - IF (allocated(rolling_precip_10min)) THEN - DO k = sub_limit(3)-1,1,-1 - rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:) - END DO - rolling_precip_10min(1,i) = rainnc(i) - ave_prate_10min(i) = max((rolling_precip_10min(1,i) - rolling_precip_10min(sub_limit(3),i))/sub_window_time(3),ave_prate_10min(i)) +! rolling_precip_1min(1,i) = rainnc(i) + rolling_precip(z,1,i) = rainnc(i) + END DO + ave_prate_1min(i) = max((rolling_precip(1,1,i) - rolling_precip(1,sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) + ave_prate_5min(i) = max((rolling_precip(2,1,i) - rolling_precip(2,sub_limit(1),i))/sub_window_time(2),ave_prate_5min(i)) + ave_prate_10min(i) = max((rolling_precip(3,1,i) - rolling_precip(3,sub_limit(1),i))/sub_window_time(3),ave_prate_10min(i)) +! ave_prate_1min(i) = max((rolling_precip_1min(1,i) - rolling_precip_1min(sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) END IF +! ! 5-minute prate +! IF (allocated(rolling_precip_5min)) THEN +! DO k = sub_limit(2)-1,1,-1 +! rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:) +! END DO +! rolling_precip_5min(1,i) = rainnc(i) +! ave_prate_5min(i) = max((rolling_precip_5min(1,i) - rolling_precip_5min(sub_limit(2),i))/sub_window_time(2),ave_prate_5min(i)) +! END IF +! ! 10-minute prate +! IF (allocated(rolling_precip_10min)) THEN +! DO k = sub_limit(3)-1,1,-1 +! rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:) +! END DO +! rolling_precip_10min(1,i) = rainnc(i) +! ave_prate_10min(i) = max((rolling_precip_10min(1,i) - rolling_precip_10min(sub_limit(3),i))/sub_window_time(3),ave_prate_10min(i)) +! END IF if(l_acrain .and. config_bucket_rainnc.gt.0._RKIND .and. & rainnc(i).gt.config_bucket_rainnc) then @@ -1271,7 +1272,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite) ! Reflectivity is already computed in 3D in the refl10cm_p array, so just read out the values if(.not.allocated(dBz1d)) allocate(dBZ1d(kts:kte)) if(.not.allocated(zp) ) allocate(zp(kts:kte) ) - !call mpas_log_write("Computing refl10cm_1km_max") + call mpas_log_write("Computing refl10cm_1km_max") !call mpas_log_write("max before = $r dbz", realArgs=(/maxval(refl10cm_1km_max)/)) do j = jts,jte ! Assuming here that jts == jte (2D arrays) do i = its,ite From cb5a1a6d89bb172198352b987a88afb3e8c47c86 Mon Sep 17 00:00:00 2001 From: Jeffrey Duda Date: Thu, 6 Nov 2025 18:19:57 +0000 Subject: [PATCH 5/6] Add IF statements to control when maximum precipitation rate calculations should be computed based on choice of microphysics scheme. --- .../mpas_atmphys_driver_microphysics.F | 129 +++++++++--------- 1 file changed, 68 insertions(+), 61 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index a5fb0eaa68..6171f5ca19 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -1023,15 +1023,16 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) 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_1min,rolling_precip_5min,rolling_precip_10min - real(kind=RKIND),dimension(:,:,:),allocatable,save:: rolling_precip + real(kind=RKIND),dimension(:,:),allocatable,save:: rolling_precip !local variables and arrays: character(len=16) :: desc integer:: i,j,k,z - integer :: n_sub_windows + integer :: n_sub_windows, longest_window_n integer, dimension(5) :: sub_limit, sub_window_len ! logical, dimension(5), save :: init_rolling_precip_array = .false. logical, save :: alloc_rolling_precip_array = .false. + logical :: calc_prate_max real(kind=RKIND), dimension(5) :: n_dt_sub, sub_window_time real(kind=RKIND):: rho_a @@ -1055,31 +1056,37 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) 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) - n_sub_windows = 3 - 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 - if (.not. alloc_rolling_precip_array) then - allocate(rolling_precip(maxval(/sub_limit(1),sub_limit(2),sub_limit(3)/),its:ite)) - rolling_precip = 0._RKIND - alloc_rolling_precip_array = .true. - end if + calc_prate_max = .false. + if (microp_scheme == "mp_thompson" .or. microp_scheme == "mp_tempo" .or. microp_scheme == "mp_thompson_aerosols" .or. microp_scheme == "mp_nssl2m" .or. microp_scheme == "mp_wsm6") calc_prate_max = .true. + + if (calc_prate_max) THEN + n_sub_windows = 3 + 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 @@ -1097,44 +1104,44 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) !time-step precipitation: rainncv(i) = rainnc_p(i,j) - prate_max(i) = max(rainnc_p(i,j)/dt_microp,prate_max(i)) + 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) - ! update rolling precipitation totals for sub-history-interval average-precipitation-rate calculation - ! 1-minute prate - IF (allocated(rolling_precip)) THEN -! IF (allocated(rolling_precip_1min)) THEN - DO z = 1,3 - DO k = sub_limit(z)-1,1,-1 -! rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) - rolling_precip(z,k+1,:) = rolling_precip(z,k,:) - END DO -! rolling_precip_1min(1,i) = rainnc(i) - rolling_precip(z,1,i) = rainnc(i) - END DO - ave_prate_1min(i) = max((rolling_precip(1,1,i) - rolling_precip(1,sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) - ave_prate_5min(i) = max((rolling_precip(2,1,i) - rolling_precip(2,sub_limit(1),i))/sub_window_time(2),ave_prate_5min(i)) - ave_prate_10min(i) = max((rolling_precip(3,1,i) - rolling_precip(3,sub_limit(1),i))/sub_window_time(3),ave_prate_10min(i)) -! ave_prate_1min(i) = max((rolling_precip_1min(1,i) - rolling_precip_1min(sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) - END IF -! ! 5-minute prate -! IF (allocated(rolling_precip_5min)) THEN -! DO k = sub_limit(2)-1,1,-1 -! rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:) -! END DO -! rolling_precip_5min(1,i) = rainnc(i) -! ave_prate_5min(i) = max((rolling_precip_5min(1,i) - rolling_precip_5min(sub_limit(2),i))/sub_window_time(2),ave_prate_5min(i)) -! END IF -! ! 10-minute prate -! IF (allocated(rolling_precip_10min)) THEN -! DO k = sub_limit(3)-1,1,-1 -! rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:) -! END DO -! rolling_precip_10min(1,i) = rainnc(i) -! ave_prate_10min(i) = max((rolling_precip_10min(1,i) - rolling_precip_10min(sub_limit(3),i))/sub_window_time(3),ave_prate_10min(i)) -! END IF + 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 +! IF (allocated(rolling_precip_1min)) THEN + DO k = longest_window_n-1,1,-1 +! rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) + rolling_precip(k+1,:) = rolling_precip(k,:) + END DO +! rolling_precip_1min(1,i) = rainnc(i) + 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)) + 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)) +! ave_prate_1min(i) = max((rolling_precip_1min(1,i) - rolling_precip_1min(sub_limit(1),i))/sub_window_time(1),ave_prate_1min(i)) + END IF +! ! 5-minute prate +! IF (allocated(rolling_precip_5min)) THEN +! DO k = sub_limit(2)-1,1,-1 +! rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:) +! END DO +! rolling_precip_5min(1,i) = rainnc(i) +! ave_prate_5min(i) = max((rolling_precip_5min(1,i) - rolling_precip_5min(sub_limit(2),i))/sub_window_time(2),ave_prate_5min(i)) +! END IF +! ! 10-minute prate +! IF (allocated(rolling_precip_10min)) THEN +! DO k = sub_limit(3)-1,1,-1 +! rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:) +! END DO +! rolling_precip_10min(1,i) = rainnc(i) +! ave_prate_10min(i) = max((rolling_precip_10min(1,i) - rolling_precip_10min(sub_limit(3),i))/sub_window_time(3),ave_prate_10min(i)) +! END IF + END IF ! calc_prate_max if(l_acrain .and. config_bucket_rainnc.gt.0._RKIND .and. & rainnc(i).gt.config_bucket_rainnc) then From def3905bb17b7392f9e504dd2e223eb91de2accb Mon Sep 17 00:00:00 2001 From: Jeffrey Duda Date: Tue, 16 Dec 2025 20:42:19 +0000 Subject: [PATCH 6/6] Corrected bug for sub-hourly max prate calculation --- .../physics/mpas_atmphys_driver_microphysics.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index 6171f5ca19..a32c6f6271 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -1115,8 +1115,8 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) IF (allocated(rolling_precip)) THEN ! IF (allocated(rolling_precip_1min)) THEN DO k = longest_window_n-1,1,-1 -! rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:) - rolling_precip(k+1,:) = rolling_precip(k,:) +! rolling_precip_1min(k+1,i) = rolling_precip_1min(k,i) + rolling_precip(k+1,i) = rolling_precip(k,i) END DO ! rolling_precip_1min(1,i) = rainnc(i) rolling_precip(1,i) = rainnc(i)