diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index f469dd438a..731699a5ef 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -682,6 +682,7 @@ + @@ -1142,6 +1143,7 @@ + @@ -1149,6 +1151,7 @@ + @@ -3082,6 +3085,10 @@ description="10 cm maximum radar reflectivity" packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in;mp_nssl2m_in;mp_tempo_in"/> + + @@ -3114,6 +3121,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"/> + + + + + + @@ -3126,6 +3145,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 3f628a0cce..523877ff85 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -982,13 +982,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 diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index d78d56bec3..a32c6f6271 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) ) @@ -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) 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: - 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. & @@ -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 +! 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') @@ -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 + 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 :: 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 !----------------------------------------------------------------------------------------------------------------- @@ -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. + 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 @@ -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(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 +! IF (allocated(rolling_precip_1min)) THEN + 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) + 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 + 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) 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) 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 + real(kind=RKIND),dimension(:),pointer:: refl10cm_1km_max, refl10cm_1km, refl10cm_max !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_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 - 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 35d577efab..8303e1469c 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 @@ -1171,7 +1173,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 @@ -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 b379086243..5fe236d5e6 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F @@ -389,7 +389,8 @@ module mpas_atmphys_vars graupelnc_p, &! graupelncv_p, &! sr_p, &! - nwfa2d_p + nwfa2d_p, & + prate_p real(kind=RKIND),dimension(:,:),allocatable:: & frainnc_p @@ -410,7 +411,8 @@ module mpas_atmphys_vars !... for Thompson cloud microphysics parameterization, including aerosol-aware option: 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 !=================================================================================================================