Add prate_max and other related precip variables to MPAS output#151
Add prate_max and other related precip variables to MPAS output#151Jeff-Duda wants to merge 12 commits intoufs-community:gsl/developfrom
Conversation
…t output time): -Instantaneous precipitation rate -Composite reflectivity -1-minute average precipitation rate -5-minute average precipitation rate -10-minute average precipitation rate
jderrico-noaa
requested review from
AndersJensen-NOAA,
XiaSun-Atmos,
barlage,
joeolson42 and
mdtoyNOAA
|
@Jeff-Duda thanks for the PR -- I'm glad the painful authentication issues have been squared away! I'm seeing some conflicts with mpas_atm_core.F and mpas_atmphys_interface.F, all of which appear to tie to Anders' hail diagnostics added with v8.3.0-1.4. Could you please update your branch to v8.3.0-1.8 at your earliest convenience? That should resolve these conflicts. |
|
@clark-evans, so, I went through this late on Friday with @jderrico-noaa . She showed me my merge conflicts and I did correct those in my repo. I committed them and ran the push step, but it seems something hasn't gone right yet. How do I know what I still need to do to correct this? |
I can see it here on the web -- at the bottom of this discussion, it shows "This branch has conflicts that must be resolved." It gives me the option of looking at them in the web editor, which is how I was able to find where they exist. That doesn't get everything quite resolved, though. |
Oddly enough, I just got an email with the merge of gsl/develop into your branch...from five days ago. Not sure what's going on with that, but glad it resolved itself! There will be one last update to do before merging, and either Janet or I will give you a heads-up once your PR reaches to the top of the queue for review/merging. In the meantime, could you please run the prePR set of regression tests from Mike Barlage's mpas_testcase repo? Fuller documentation is available in this Google doc. |
I am running into severe errors trying to do this using that document. I do not know how to fix them. |
|
@Jeff-Duda Thanks for getting this PR together. Is your new refl10cm max variable different from ours, which is named refl10cm_max? I'll take a look at the code and provide feedback soon. |
Yes, @AndersJensen-NOAA , the new variable is a time-maximum of the composite (vertical-max) reflectivity, whereas refl10cm_max in the pre-existing code is instantaneous composite reflectivity. |
…e_max Correcting character string length error in microphysics driver file.
|
Is this PR waiting on any action from me? |
@Jeff-Duda The PR is in the queue. No actions from you are needed at this time. There may be some changes requested during review though. For example, we'll likely request the addition of config flags so that these diagnostics can be turned on or off. If you are not familiar with how to do that in MPAS, we can help. |
|
@Jeff-Duda Can you pull in the latest changes from gsl-develop into your PR branch? Our latest code has some automatic tests that will run. |
| init_rolling_precip_array(z) = .true. | ||
| end if | ||
| end do | ||
| ! if all is okay, allocate rolling precip arrays just once and never again |
There was a problem hiding this comment.
For lines 1074-1088, instead of repeating code, can this be in a loop?
There was a problem hiding this comment.
@AndersJensen-NOAA I can do this, but it will require a different array to hold all precip intervals in one array, which will render the check on line ~1068 moot; the code will not be able to distinguish whether the model time step is too coarse for just the 1-minute precip (but okay for 5- and 10-minute rolling precip) or too coarse for 1- and 5-minute precip but okay for 10-minute precip. I don't know how critical this warning actually is, though, as it is mostly self-explanatory and easy to account for.
| END IF | ||
| ! 10-minute prate | ||
| IF (allocated(rolling_precip_10min)) THEN | ||
| DO k = sub_limit(3)-1,1,-1 |
There was a problem hiding this comment.
Is it possible to check all allocations at once and then combine these since sub_limit(1) < sub_limit(2) < sub_limit(3)?
DO k = sub_limit(3)-1,1-1
rolling_precip_10min(k+1,:) = rolling_precip_10min(k,:)
IF k <= sub_limit(2) rolling_precip_5min(k+1,:) = rolling_precip_5min(k,:)
IF k <= sub_limit(1) rolling_precip_1min(k+1,:) = rolling_precip_1min(k,:)
@AndersJensen-NOAA How do I do that? |
I went ahead and took care of that for you through the web interface. You'll want to make sure to sync your local copy, though. |
What does it mean to "sync my local copy"? |
The copy at https://github.com/Jeff-Duda/MPAS-Model/tree/prate_max has been updated with the latest gsl/develop, so you'll want to go to where you were developing it (presumably on disk somewhere) and sync with the remote. See also: https://stackoverflow.com/questions/6373277/synchronizing-a-local-git-repository-with-a-remote-one. |
|
My apologies to everyone who received a request to review just now. That was an indavertent click. |
|
@clark-evans and @AndersJensen-NOAA I have made some code changes, committed them, git-pulled to sync my local fork, then git-pushed the changes back. Did I miss anything? What is currently needed from me? |
Hmmm, it doesn't look like the remote repo has yet updated with your changes -- the last commit is my syncing of the repo with gsl/develop ~5 h ago. |
|
Here's a record of what commands I issued: |
|
Sorry, I see now I forgot to re-enter 'git push'. I have just done that. |
…ions should be computed based on choice of microphysics scheme.
|
It has been awhile since this PR was addressed, unfortunately due to the shutdown. But I did make code changes and committed them. But then I got some errors for whatever tests are being automatically conducted. There seems to remain only one error, and I doubt that it is within my control to fix. See the attached image. 
|
Hi @Jeff-Duda -- we see that error crop up every so often here. Could you try to manually re-run the test? You should be able to do so from Github on the web - go to the PR, then choose "Checks," then choose the failing job, then at upper right you should have a button related to re-running tests. I select the "re-run failed checks" option, and that usually does the trick. |
|
Thanks @clark-evans . The "checks" tab has 0 issues for me (I don't even see an option to re-run anything). Does this mean there actually are no problems now? |
You may need to do it in your repo (e.g., from https://github.com/Jeff-Duda/MPAS-Model/actions/runs/19145980522/job/54723565349), as I don't think the testing workflows in this repo have yet been triggered for your PR. I think your codebase acquired the workflows after your last merge of gsl/develop into it, whereas on our end I have a "4 workflows awaiting approval" message on your PR. |
|
Thanks, @clark-evans! Looks like the tests are now successful! @AndersJensen-NOAA, perhaps that means you can now re-evaluate? |
|
What needs to happen here to get this PR moving again? It seems to have stalled. |
|
@clark-evans @Jeff-Duda |
| !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) ) |
There was a problem hiding this comment.
This isn't used -- delete
| !surface precipitation: | ||
| if(allocated(rainnc_p) ) deallocate(rainnc_p ) | ||
| if(allocated(rainncv_p)) deallocate(rainncv_p) | ||
| if(allocated(prate_p)) deallocate(prate_p) |
There was a problem hiding this comment.
This isn't used -- delete
| 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 |
There was a problem hiding this comment.
undo this change. l_diags is likely used by others.
There was a problem hiding this comment.
@AndersJensen-NOAA I must protest this required change.
While you may be right about l_diags, the problem is that with this proposed code, all fields that are advertised as "hourly maximum reflectivity" are in fact, not that. They are the maximum over just two time steps, one at each end of an hour. This is improper and unhelpful. In order to get a proper analysis of time-maximum reflectivity, the reflectivity must be computed at each time step.
This change is attempting to fix an error rather than make a preferential change.
| call compute_relhum(diag,its,ite) | ||
|
|
||
| end if | ||
| ! end if |
| 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') |
| sr_p, &! | ||
| nwfa2d_p | ||
| nwfa2d_p, & | ||
| prate_p |
| real(kind=RKIND),dimension(:,:),allocatable:: & | ||
| refl10cm_1km_p | ||
| refl10cm_1km_p, & | ||
| refl10cm_max_p |
| <var name="max_hail_diameter_column"/> | ||
| <var name="max_hail_diameter_column_acc"/> | ||
| <var name="refl10cm_max"/> | ||
| <var name="refl10cm_compref_max"/> |
| <var name="swupt"/> | ||
| <var name="rainc"/> | ||
| <var name="rainnc"/> | ||
| <var name="prate_max"/> |
| <var name="max_hail_diameter_column"/> | ||
| <var name="max_hail_diameter_column_acc"/> | ||
| <var name="refl10cm_max"/> | ||
| <var name="refl10cm_compref_max"/> |
dustinswales
left a comment
dustinswales
left a comment
There was a problem hiding this comment.
@Jeff-Duda There are some issues that need to get resolved, but nothing major.
FWIW. For MPAS in the UWM, we already have a module where these (hourly) diagnostics are computed.
IF we need these (sub-hourly) diagnostics for MPAS in the UWM, we will need to compute them in this module.
| 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 |
There was a problem hiding this comment.
This l_diags block has two calculations. Reflectivity is accumulated, so you need to call it every time, but relhum is ephemeral and only needs to be called when l_diags is true.
Please split this block into two.
| 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') |
| character(len=16) :: desc | ||
| integer:: i,j,k,z | ||
| integer :: n_sub_windows, longest_window_n | ||
| integer, dimension(5) :: sub_limit, sub_window_len |
There was a problem hiding this comment.
Please parameterize magic number "5"
Also, I think it should be magic number "3"
| 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 |
There was a problem hiding this comment.
Just above, sub_limit, sub_window_len, n_dt_sub, and sub_window_time are defined with length of 5.
But here you hardcode n_sub_wondows = 3, and only the first three elements are used, for 1, 5 and 10min averaging.
Should n_sub_windows, and sub_window_len, be parameters?
| ! 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 |
| ! 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) |
| 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)) |
There was a problem hiding this comment.
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_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 |
| 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 |
There was a problem hiding this comment.
This doesn't appears to be used anywhere?
…dded for idealized testing.
|
I committed and pushed a set of changes that hopefully addresses this list of requested changes. Please check whether I addressed all that was requested and let me know what I forgot to address. Thanks. |
4 participants
Reference: #145
See the below examples for more information.
MPAS-Dev/MPAS#930
MPAS-Dev/MPAS#931
Test cases included the 1.25 km and 3 km realtime configurations being run by David Dowell on Jet. I copied the input files from the 1.25 km run initialized 18 June 2025 at 06Z and from the 3 km run initialized at 14 July 2025 at 00Z. The requested output files are present in the diagnc files for both runs, and futher MPASSIT runs show these arrays look good and appear to display the requested information.
Case comparison test results have run, but this webpage won't let me attach the file and it is very long to paste here. How shall I share the results?