Add prate_max and other related precip variables to MPAS output #151
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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) ) | ||
|
Collaborator
There was a problem hiding this comment. This isn't used -- delete |
||
|
|
||
| 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) ) | ||
|
|
@@ -302,6 +304,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) | ||
|
Collaborator
There was a problem hiding this comment. This isn't used -- delete |
||
|
|
||
| microp_select: select case(trim(microp_scheme)) | ||
| case ("mp_thompson","mp_thompson_aerosols","mp_wsm6") | ||
|
|
@@ -378,6 +381,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) | ||
|
|
@@ -838,7 +842,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: | ||
| !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 | ||
|
Collaborator
There was a problem hiding this comment. undo this change. l_diags is likely used by others.
Author
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.
clark-evans marked this conversation as resolved.
Outdated
|
||
| if(trim(microp_scheme) == "mp_wsm6" .or. & | ||
| trim(microp_scheme) == "mp_thompson" .or. & | ||
| trim(microp_scheme) == "mp_thompson_aerosols" .or. & | ||
|
|
@@ -855,12 +861,12 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten | |
| !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') | ||
|
|
@@ -926,6 +932,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 | ||
|
|
||
|
|
@@ -1010,12 +1017,23 @@ 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 | ||
|
Collaborator
There was a problem hiding this comment. remove commented line |
||
| 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, longest_window_n | ||
| integer, dimension(5) :: sub_limit, sub_window_len | ||
|
clark-evans marked this conversation as resolved.
Outdated
|
||
| ! logical, dimension(5), save :: init_rolling_precip_array = .false. | ||
|
Collaborator
There was a problem hiding this comment. remove commented line |
||
| 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 | ||
|
|
||
| !----------------------------------------------------------------------------------------------------------------- | ||
|
|
@@ -1029,10 +1047,46 @@ 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) | ||
|
|
||
| calc_prate_max = .false. | ||
|
Collaborator
There was a problem hiding this comment. move calc_prate_max to the registry so that users can set the value to true or 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. | ||
|
Collaborator
There was a problem hiding this comment. delete line 1060 |
||
|
|
||
| if (calc_prate_max) THEN | ||
| n_sub_windows = 3 | ||
|
clark-evans marked this conversation as resolved.
Outdated
|
||
| 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 | ||
|
|
@@ -1050,15 +1104,50 @@ 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 (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 | ||
|
clark-evans marked this conversation as resolved.
Outdated
|
||
| DO k = longest_window_n-1,1,-1 | ||
| ! 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) | ||
|
clark-evans marked this conversation as resolved.
Outdated
|
||
| 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)) | ||
|
clark-evans marked this conversation as resolved.
|
||
| 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 | ||
|
clark-evans marked this conversation as resolved.
Outdated
|
||
| ! 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 | ||
|
Collaborator
There was a problem hiding this comment. delete commented lines |
||
| ! ! 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. & | ||
|
Collaborator
There was a problem hiding this comment. move 1146-1150 back to above the calc_prate_max code |
||
| rainnc(i).gt.config_bucket_rainnc) then | ||
| i_rainnc(i) = i_rainnc(i) + 1 | ||
| rainnc(i) = rainnc(i) - config_bucket_rainnc | ||
| endif | ||
|
|
||
| enddo | ||
| enddo | ||
|
|
@@ -1174,6 +1263,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite) | |
| enddo | ||
|
|
||
| refl10cm_max(i) = maxval(dBZ1d(:)) | ||
| refl10cm_max_p(i,j) = refl10cm_max(i) | ||
|
Collaborator
There was a problem hiding this comment. Is this used? you can delete this |
||
| w1 = (zp(kp+1)-1000.)/(zp(kp+1)-zp(kp)) | ||
| w2 = 1.0 - w1 | ||
| refl10cm_1km(i) = w1*dBZ1d(kp) + w2*dBZ1d(kp+1) | ||
|
|
@@ -1207,6 +1297,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite) | |
| enddo | ||
|
|
||
| refl10cm_max(i) = maxval(dBZ1d(:)) | ||
| refl10cm_max_p(i,j) = refl10cm_max(i) | ||
|
Collaborator
There was a problem hiding this comment. delete this -- I don't think it is used |
||
| w1 = (zp(kp+1)-1000.)/(zp(kp+1)-zp(kp)) | ||
| w2 = 1.0 - w1 | ||
| refl10cm_1km(i) = w1*dBZ1d(kp) + w2*dBZ1d(kp+1) | ||
|
|
@@ -1340,7 +1431,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, refl10cm_max | ||
|
clark-evans marked this conversation as resolved.
Outdated
|
||
|
|
||
| !local variables and arrays: | ||
| integer:: i,j,k,kp | ||
|
|
@@ -1352,7 +1443,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_1km',refl10cm_1km) | ||
| call mpas_pool_get_array(diag_physics,'refl10cm_1km_max',refl10cm_1km_max) | ||
|
|
||
|
|
@@ -1398,7 +1489,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 | ||
|
Collaborator
There was a problem hiding this comment. add this missing line back |
||
| case default | ||
|
|
||
| end select microp_select | ||
|
|
||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fix spacing
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@Jeff-Duda could you address this comment from Anders? It won't show up in the terminal, but it does here because of how Github handles tabs vs. spaces. Could you please match whichever mechanism of creating the spacing is used in the lines above?