Adding hourly- (time-)maximum composite reflectivity to output

[Jeff_Duda]

I want to add a time-maximum composite reflectivity to be output from MPAS...specifically for the TEMPO microphysics. I see that 3D reflectivity is already calculated in the microphysics driver, and there's even a subroutine for calculating the hourly max values of 1-km AGL reflectivity. However, when I modified the code to track a max composite reflectivity value, the new field in the diag files is identical to the composite reflectivity that is already calculated by the existing code.

As a quick aside, I recently modified the MPAS code to output hourly-max precipitation rate and that product looks like what I expected, so I think I was successful in making the code changes for that particular field. So I have been trying to model the changes I made for the maximum precipitation rate for the maximum of composite reflectivity. But the two fields have a different trajectory through the model code:

• Precipitation rate can be directly derived from the rainncv variable which is calculated directly in the microphysics scheme subroutines (a level below the driver). Once back in the microphysics_driver, rainncv is not modified.
• On the other hand, 3D reflectivity is calculated in the individual microphysics scheme subroutines, but the computation to calculate the vertical max is performed in a separate subroutine towards the end of the driver code.

I noted the use of a

[var_name]_p

nomenclature for a lot of microphysics scheme variables. These are declared in mpas_atmphys_vars.F and presumably allow for certain arrays to be tracked outside of the individual microphysics scheme subroutines. So I do have this declared for the maximum composite reflectivity by tracking just the composite reflectivity through time in mpas_atmphys_vars.F:

 real(kind=RKIND),dimension(:,:),allocatable:: &
    refl10cm_1km_p, &
    refl10cm_max_p

This section at the end of the microphysics driver is where I could use the most help in understanding:

 !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
    if(trim(microp_scheme) == "mp_wsm6"     .or. &
       trim(microp_scheme) == "mp_thompson" .or. &
       trim(microp_scheme) == "mp_thompson_aerosols" .or. &
       trim(microp_scheme) == "mp_tempo" .or. &
       trim(microp_scheme) == "mp_nssl2m"   ) then
!       call mpas_log_write('Computing 3d reflectivity')
       call compute_radar_reflectivity(configs,diag_physics,its,ite)
    else
       call mpas_log_write('*** NOTICE: NOT computing simulated radar reflectivity')
       call mpas_log_write('            since WSM6 or Thompson microphysics scheme was not selected')
    end if

    !calculate the relative humidity over water if the temperature is strictly greater than 0.C,
    !over ice otherwise.
    call compute_relhum(diag,its,ite)

 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 compute_hourly_max_radar_reflectivity(configs,diag_physics,its,ite)
    else
       call mpas_log_write('*** NOTICE: NOT computing hourly max simulated radar reflectivity')
       call mpas_log_write('            since WSM6 or Thompson microphysics scheme was not selected')
    end if
 endif

 if(trim(microp_scheme) == "mp_tempo") then
    call mpas_log_write('Computing accumulated max hail diameter at the surface')
    call compute_accumulated_max_hail(diag_physics,its,ite)
 endif

!... copy updated precipitation from the wrf-physics grid back to the geodesic-dynamics grid:
 call precip_to_MPAS(configs,diag_physics,its,ite)

!... copy updated cloud microphysics variables from the wrf-physics grid back to the geodesic-
!    dynamics grid:
 call microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,tend_physics,tend,its,ite)

From comments higher up in the driver source code file it seems like the "whatever_to_MPAS" subroutines covert the local variables to variables that the higher-level MPAS code needs to output them in the history files. There is one for precip and one for microphysics. The crux of my code additions is this:

                max_compref(i) = max(max_compref(i),refl10cm_max_p(i,j))

The uncertainty in my code then is where to place this command. For the hourly-max precipitation rate it is located in the precip_to_MPAS subroutine:

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

...and that seems to be working. But putting the analogous line in microphysics_to_MPAS (in a select-case statement for 'mp_tempo') is resulting in that field looking the same as composite reflectivity. So I don't know what I'm missing and could use some help. Thanks.