From a827edaeb3a8924c570bc99a399f7867e8d8a867 Mon Sep 17 00:00:00 2001 From: Wei Wang Date: Thu, 13 Feb 2020 19:24:40 -0700 Subject: [PATCH 1/3] separating nwp diagnostics from the rest - in preparation for adding accumulated tendency changes --- phys/Makefile | 1 + phys/module_diag_misc.F | 720 +------------------------- phys/module_diag_nwp.F | 835 +++++++++++++++++++++++++++++++ phys/module_diagnostics_driver.F | 267 +++------- 4 files changed, 909 insertions(+), 914 deletions(-) create mode 100644 phys/module_diag_nwp.F diff --git a/phys/Makefile b/phys/Makefile index 4218640e8a..1591460797 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -207,6 +207,7 @@ DIAGNOSTIC_MODULES_EM = \ module_diag_functions.o \ module_diag_hailcast.o \ module_diag_misc.o \ + module_diag_nwp.o \ module_diag_rasm.o \ module_diag_pld.o \ module_diag_zld.o \ diff --git a/phys/module_diag_misc.F b/phys/module_diag_misc.F index bea48b87fc..24300d45ce 100644 --- a/phys/module_diag_misc.F +++ b/phys/module_diag_misc.F @@ -9,8 +9,6 @@ END MODULE module_diag_misc ! MODULE module_diag_misc - PRIVATE :: WGAMMA - PRIVATE :: GAMMLN CONTAINS SUBROUTINE diagnostic_output_calc( & ids,ide, jds,jde, kds,kde, & @@ -19,10 +17,9 @@ SUBROUTINE diagnostic_output_calc( & i_start,i_end,j_start,j_end,kts,kte,num_tiles & ,dpsdt,dmudt & ,p8w,pk1m,mu_2,mu_2m & - ,u,v, temp & ,raincv,rainncv,rainc,rainnc & ,i_rainc,i_rainnc & - ,hfx,sfcevp,lh & + ,hfx,sfcevp,lh,t2 & ,ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC & ! Optional ,ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC & ! Optional ,ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC & ! Optional @@ -31,52 +28,25 @@ SUBROUTINE diagnostic_output_calc( & ,I_ACSWUPB,I_ACSWUPBC,I_ACSWDNB,I_ACSWDNBC & ! Optional ,I_ACLWUPT,I_ACLWUPTC,I_ACLWDNT,I_ACLWDNTC & ! Optional ,I_ACLWUPB,I_ACLWUPBC,I_ACLWDNB,I_ACLWDNBC & ! Optional - ,dt,xtime,sbw,t2 & + ,dt,xtime & ,diag_print & ,bucket_mm, bucket_J & - ,mphysics_opt & - ,gsfcgce_hail, gsfcgce_2ice & - ,mpuse_hail & - ,nssl_cnoh, nssl_cnohl & - ,nssl_rho_qh, nssl_rho_qhl & - ,nssl_alphah, nssl_alphahl & ,prec_acc_c, prec_acc_nc, snow_acc_nc & ,snowncv, prec_acc_dt, curr_secs2 & - ,nwp_diagnostics, diagflag & ,history_interval & ,itimestep & - ,u10,v10,w & - ,wspd10max & - ,up_heli_max & - ,w_up_max,w_dn_max & - ,znw,w_colmean & - ,numcolpts,w_mean & - ,grpl_max,grpl_colint,refd_max,refl_10cm & - ,hail_maxk1,hail_max2d & - ,qg_curr & - ,ng_curr,qh_curr,nh_curr,qr_curr,nr_curr & ! Optional (gthompsn) - ,rho,ph,phb,g & - ,max_time_step,adaptive_ts & ) !---------------------------------------------------------------------- + USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval - USE module_state_description, ONLY : & - KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME, & - WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, & - MORR_TWO_MOMENT, GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, & - NSSL_2MOM, NSSL_2MOMG, NSSL_2MOMCCN, NSSL_1MOM, NSSL_1MOMLFO, & - MILBRANDT2MOM , CAMMGMPSCHEME, FAST_KHAIN_LYNN, FULL_KHAIN_LYNN, & - MORR_TM_AERO !TWG 2017 !,MILBRANDT3MOM, NSSL_3MOM - IMPLICIT NONE + IMPLICIT NONE !====================================================================== ! Definitions !----------- !-- DIAG_PRINT print control: 0 - no diagnostics; 1 - dmudt only; 2 - all !-- DT time step (second) !-- XTIME forecast time -!-- SBW specified boundary width - used later -! !-- P8W 3D pressure array at full eta levels !-- MU dry column hydrostatic pressure !-- RAINC cumulus scheme precipitation since hour 0 @@ -87,26 +57,11 @@ SUBROUTINE diagnostic_output_calc( & !-- HFX surface sensible heat flux !-- LH surface latent heat flux !-- SFCEVP total surface evaporation -!-- U u component of wind - to be used later to compute k.e. -!-- V v component of wind - to be used later to compute k.e. !-- PREC_ACC_C accumulated convective precip over accumulation time prec_acc_dt !-- PREC_ACC_NC accumulated explicit precip over accumulation time prec_acc_dt !-- SNOW_ACC_NC accumulated explicit snow precip over accumulation time prec_acc_dt !-- PREC_ACC_DT precip accumulation time, default is 60 min !-- CURR_SECS2 Time (s) since the beginning of the restart -!-- NWP_DIAGNOSTICS = 1, compute hourly maximum fields -!-- DIAGFLAG logical flag to indicate if this is a history output time -!-- U10, V10 10 m wind components -!-- WSPD10MAX 10 m max wind speed -!-- UP_HELI_MAX max updraft helicity -!-- W_UP_MAX max updraft vertical velocity -!-- W_DN_MAX max downdraft vertical velocity -!-- W_COLMEAN column mean vertical velocity -!-- NUMCOLPTS no of column points -!-- GRPL_MAX max column-integrated graupel -!-- GRPL_COLINT column-integrated graupel -!-- REF_MAX max derived radar reflectivity -!-- REFL_10CM model computed 3D reflectivity ! !-- ids start index for i in domain !-- ide end index for i in domain @@ -144,20 +99,11 @@ SUBROUTINE diagnostic_output_calc( & INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & & i_start,i_end,j_start,j_end - INTEGER, INTENT(IN ) :: diag_print + INTEGER, INTENT(IN ) :: diag_print REAL, INTENT(IN ) :: bucket_mm, bucket_J - INTEGER, INTENT(IN ) :: mphysics_opt - INTEGER, INTENT(IN) :: gsfcgce_hail, gsfcgce_2ice, mpuse_hail - REAL, INTENT(IN) :: nssl_cnoh, nssl_cnohl & - ,nssl_rho_qh, nssl_rho_qhl & - ,nssl_alphah, nssl_alphahl - INTEGER, INTENT(IN ) :: nwp_diagnostics - LOGICAL, INTENT(IN ) :: diagflag REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & - INTENT(IN ) :: u & - , v & - , p8w + INTENT(IN ) :: p8w REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: & MU_2 & @@ -178,7 +124,6 @@ SUBROUTINE diagnostic_output_calc( & , PK1M REAL, INTENT(IN ) :: DT, XTIME - INTEGER, INTENT(IN ) :: SBW INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: & I_RAINC, & I_RAINNC @@ -212,62 +157,8 @@ SUBROUTINE diagnostic_output_calc( & INTEGER, INTENT(IN) :: & history_interval,itimestep - REAL, DIMENSION( kms:kme ), INTENT(IN) :: & - znw - - REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: & - w & - ,temp & - ,qg_curr & - ,rho & - ,refl_10cm & - ,ph,phb - - REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: & - ng_curr, qh_curr, nh_curr & - ,qr_curr, nr_curr - - REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: & - u10 & - ,v10 - - REAL, INTENT(IN) :: g - - REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: & - wspd10max & - ,up_heli_max & - ,w_up_max,w_dn_max & - ,w_colmean,numcolpts,w_mean & - ,grpl_max,grpl_colint & - ,hail_maxk1,hail_max2d & - ,refd_max - - REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: temp_qg, temp_ng, temp_qr, temp_nr - INTEGER :: idump - REAL :: wind_vel - REAL :: depth - - DOUBLE PRECISION:: hail_max - REAL:: hail_max_sp - DOUBLE PRECISION, PARAMETER:: thresh_conc = 0.0005d0 ! number conc. of graupel/hail per cubic meter - LOGICAL:: scheme_has_graupel - INTEGER, PARAMETER:: ngbins=50 - DOUBLE PRECISION, DIMENSION(ngbins+1):: xxDx - DOUBLE PRECISION, DIMENSION(ngbins):: xxDg, xdtg - REAL:: xrho_g, xam_g, xbm_g, xmu_g - REAL, DIMENSION(3):: cge, cgg - DOUBLE PRECISION:: f_d, sum_ng, sum_t, lamg, ilamg, N0_g, lam_exp, N0exp - DOUBLE PRECISION:: lamr, N0min - REAL:: mvd_r, xslw1, ygra1, zans1 - INTEGER:: ng, n - - REAL :: time_from_output - INTEGER :: max_time_step - LOGICAL :: adaptive_ts - LOGICAL :: reset_arrays = .FALSE. - !----------------------------------------------------------------- ! Handle accumulations with buckets to prevent round-off truncation in long runs ! This is done every 360 minutes assuming time step fits exactly into 360 minutes @@ -440,567 +331,6 @@ SUBROUTINE diagnostic_output_calc( & ! !$OMP END PARALLEL DO ENDIF - - -! NSSL - - IF ( nwp_diagnostics .EQ. 1 ) THEN - - idump = (history_interval * 60.) / dt - -! IF ( MOD(itimestep, idump) .eq. 0 ) THEN -! WRITE(outstring,*) 'Computing PH0 for this domain with curr_secs2 = ', curr_secs2 -! CALL wrf_message ( TRIM(outstring) ) - - time_from_output = mod( xtime, REAL(history_interval) ) - -! print *,' history_interval = ', history_interval -! print *,' itimestep = ', itimestep -! print *,' idump = ', idump -! print *,' xtime = ', xtime -! print *,' time_from_output = ', time_from_output -! print *,' max_time_step = ', max_time_step - - IF ( adaptive_ts .EQV. .TRUE. ) THEN -! if timestep is adaptive, use new resetting method; -! also, we are multiplying max_time_step with 1.05; because of rounding error, -! the time_from_output can become slightly larger than max_time_step when -! adaptive time step is maximized, in which case the if condition below fails to detect -! that we just wrote an output - IF ( ( time_from_output .GT. 0 ) .AND. ( time_from_output .LE. ( ( max_time_step * 1.05 ) / 60. ) ) ) THEN - reset_arrays = .TRUE. - ENDIF - ELSE -! if timestep is fixed, use original resetting method - IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN - reset_arrays = .TRUE. - ENDIF - ENDIF - -! print *,' reset_arrays = ', reset_arrays - IF ( reset_arrays .EQV. .TRUE. ) THEN - -! print *,' Resetting NWP max arrays ' - - WRITE(outstring,*) 'NSSL Diagnostics: Resetting max arrays for domain with dt = ', dt - CALL wrf_debug ( 10,TRIM(outstring) ) - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO i=i_start(ij),i_end(ij) - wspd10max(i,j) = 0. - up_heli_max(i,j) = 0. - w_up_max(i,j) = 0. - w_dn_max(i,j) = 0. - w_mean(i,j) = 0. - grpl_max(i,j) = 0. - refd_max(i,j) = 0. - hail_maxk1(i,j) = 0. - hail_max2d(i,j) = 0. - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - reset_arrays = .FALSE. - ENDIF - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO i=i_start(ij),i_end(ij) - -! Zero some accounting arrays that will be used below - - w_colmean(i,j) = 0. - numcolpts(i,j) = 0. - grpl_colint(i,j) = 0. - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme - DO i=i_start(ij),i_end(ij) - -! Find vertical velocity max (up and down) below 400 mb - - IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .GT. w_up_max(i,j) ) THEN - w_up_max(i,j) = w(i,k,j) - ENDIF - - IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .LT. w_dn_max(i,j) ) THEN - w_dn_max(i,j) = w(i,k,j) - ENDIF - -! For the column mean vertical velocity calculation, first -! total the vertical velocity between sigma levels 0.5 and 0.8 - - IF ( znw(k) .GE. 0.5 .AND. znw(k) .LE. 0.8 ) THEN - w_colmean(i,j) = w_colmean(i,j) + w(i,k,j) - numcolpts(i,j) = numcolpts(i,j) + 1 - ENDIF - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - -! Calculate the column integrated graupel - - depth = ( ( ph(i,k+1,j) + phb(i,k+1,j) ) / g ) - & - ( ( ph(i,k ,j) + phb(i,k ,j) ) / g ) - grpl_colint(i,j) = grpl_colint(i,j) + qg_curr(i,k,j) * depth * rho(i,k,j) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO i=i_start(ij),i_end(ij) - -! Calculate the max 10 m wind speed between output times - - wind_vel = sqrt ( u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j) ) - IF ( wind_vel .GT. wspd10max(i,j) ) THEN - wspd10max(i,j) = wind_vel - ENDIF - -! Calculate the column mean vertical velocity between output times - - w_mean(i,j) = w_mean(i,j) + w_colmean(i,j) / numcolpts(i,j) - - IF ( MOD(itimestep, idump) .eq. 0 ) THEN - w_mean(i,j) = w_mean(i,j) / idump - ENDIF - -! Calculate the max column integrated graupel between output times - - IF ( grpl_colint(i,j) .gt. grpl_max(i,j) ) THEN - grpl_max(i,j) = grpl_colint(i,j) - ENDIF - - ! Calculate the max radar reflectivity between output times - - IF ( refl_10cm(i,kms,j) .GT. refd_max(i,j) ) THEN - refd_max(i,j) = refl_10cm(i,kms,j) - ENDIF - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - - - -!+---+-----------------------------------------------------------------+ -!+---+-----------------------------------------------------------------+ -!..Calculate a maximum hail diameter from the characteristics of the -!.. graupel category mixing ratio and number concentration (or hail, if -!.. available). This diagnostic uses the actual spectral distribution -!.. assumptions, calculated by breaking the distribution into 50 bins -!.. from 0.5mm to 7.5cm. Once a minimum number concentration of 0.01 -!.. particle per cubic meter of air is reached, from the upper size -!.. limit, then this bin is considered the max size. -!+---+-----------------------------------------------------------------+ - - WRITE(outstring,*) 'GT-Diagnostics, computing max-hail diameter' - CALL wrf_debug (100, TRIM(outstring)) - - - IF (PRESENT(qh_curr)) THEN - WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail mixing ratio' - CALL wrf_debug (150, TRIM(outstring)) -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - temp_qg(i,k,j) = MAX(1.E-12, qh_curr(i,k,j)*rho(i,k,j)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - ELSE -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - temp_qg(i,k,j) = MAX(1.E-12, qg_curr(i,k,j)*rho(i,k,j)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - ENDIF - - IF (PRESENT(nh_curr)) THEN - WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail number concentration' - CALL wrf_debug (150, TRIM(outstring)) -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - temp_ng(i,k,j) = MAX(1.E-8, nh_curr(i,k,j)*rho(i,k,j)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - ELSEIF (PRESENT(ng_curr)) THEN - WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel number concentration' -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - temp_ng(i,k,j) = MAX(1.E-8, ng_curr(i,k,j)*rho(i,k,j)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - ELSE -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - temp_ng(i,k,j) = 1.E-8 - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - ENDIF - - - ! Calculate the max hail size from graupel/hail parameters in microphysics scheme (gthompsn 12Mar2015) - ! First, we do know multiple schemes have assumed inverse-exponential distribution with constant - ! intercept parameter and particle density. Please leave next 4 settings alone for common - ! use and copy these and cge constants to re-assign per scheme if needed (e.g. NSSL). - - xrho_g = 500. - xam_g = 3.1415926536/6.0*xrho_g ! Assumed m(D) = a*D**b, where a=PI/6*rho_g and b=3 - xbm_g = 3. ! in other words, spherical graupel/hail - xmu_g = 0. - scheme_has_graupel = .false. - - !..Some constants. These *MUST* get changed below per scheme - !.. *IF* either xbm_g or xmu_g value is changed from 3 and zero, respectively. - - cge(1) = xbm_g + 1. - cge(2) = xmu_g + 1. - cge(3) = xbm_g + xmu_g + 1. - do n = 1, 3 - cgg(n) = WGAMMA(cge(n)) - enddo - - mp_select: SELECT CASE(mphysics_opt) - - CASE (KESSLERSCHEME) -! nothing to do here - - CASE (LINSCHEME) - scheme_has_graupel = .true. - xrho_g = 917. - xam_g = 3.1415926536/6.0*xrho_g - N0exp = 4.e4 -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kme-1, kms, -1 - DO i=i_start(ij),i_end(ij) - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) - temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - - CASE (WSM3SCHEME) -! nothing to do here - - CASE (WSM5SCHEME) -! nothing to do here - - CASE (WSM6SCHEME) - scheme_has_graupel = .true. - xrho_g = 500. - xam_g = 3.1415926536/6.0*xrho_g - N0exp = 4.e6 -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kme-1, kms, -1 - DO i=i_start(ij),i_end(ij) - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) - temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - - CASE (WDM5SCHEME) -! nothing to do here - - CASE (WDM6SCHEME) - scheme_has_graupel = .true. - xrho_g = 500. - N0exp = 4.e6 - if (mpuse_hail .eq. 1) then - xrho_g = 700. - N0exp = 4.e4 - endif - xam_g = 3.1415926536/6.0*xrho_g -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kme-1, kms, -1 - DO i=i_start(ij),i_end(ij) - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) - temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - - CASE (GSFCGCESCHEME) - if (gsfcgce_2ice.eq.0 .OR. gsfcgce_2ice.eq.2) then - scheme_has_graupel = .true. - if (gsfcgce_hail .eq. 1) then - xrho_g = 900. - else - xrho_g = 400. - endif - xam_g = 3.1415926536/6.0*xrho_g - N0exp = 4.e4 -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kme-1, kms, -1 - DO i=i_start(ij),i_end(ij) - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) - temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - endif - - CASE (SBU_YLINSCHEME) -! scheme_has_graupel = .true. ! Can be calculated from rime fraction variable. -! no time to implement - - CASE (ETAMPNEW) -! scheme_has_graupel = .true. ! Can be calculated from rime fraction variable. -! no time to implement - - CASE (THOMPSON, THOMPSONAERO) - - scheme_has_graupel = .true. - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO i=i_start(ij),i_end(ij) - DO k=kme-1, kms, -1 - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - zans1 = (2.5 + 2.5/7. * (ALOG10(temp_qg(i,k,j))+7.)) - zans1 = MAX(2., MIN(zans1, 7.)) - N0exp = 10.**zans1 - lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) - lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g) - N0_g = N0exp/(cgg(2)*lam_exp) * lamg**cge(2) - temp_ng(i,k,j) = N0_g*cgg(2)*lamg**(-cge(2)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - -! CASE (MORR_MILB_P3) -! scheme_has_graupel = .true. -! either Hugh or Jason need to implement. - - CASE (MORR_TWO_MOMENT, MORR_TM_AERO) - scheme_has_graupel = .true. - xrho_g = 400. - if (mpuse_hail .eq. 1) xrho_g = 900. - xam_g = 3.1415926536/6.0*xrho_g - - CASE (MILBRANDT2MOM) - WRITE(outstring,*) 'GT-Debug, using Milbrandt2mom, which has 2-moment hail' - CALL wrf_debug (150, TRIM(outstring)) - scheme_has_graupel = .true. - xrho_g = 900. - xam_g = 3.1415926536/6.0*xrho_g - -! CASE (MILBRANDT3MOM) -! coming in future? - - CASE (NSSL_1MOMLFO, NSSL_1MOM, NSSL_2MOM, NSSL_2MOMG, NSSL_2MOMCCN) - - scheme_has_graupel = .true. - xrho_g = nssl_rho_qh - N0exp = nssl_cnoh - if (PRESENT(qh_curr)) then - xrho_g = nssl_rho_qhl - N0exp = nssl_cnohl - endif - xam_g = 3.1415926536/6.0*xrho_g - - if (PRESENT(ng_curr)) xmu_g = nssl_alphah - if (PRESENT(nh_curr)) xmu_g = nssl_alphahl - - if (xmu_g .NE. 0.) then - cge(1) = xbm_g + 1. - cge(2) = xmu_g + 1. - cge(3) = xbm_g + xmu_g + 1. - do n = 1, 3 - cgg(n) = WGAMMA(cge(n)) - enddo - endif - - ! NSSL scheme has many options, but, if single-moment, just fill - ! in the number array for graupel from built-in assumptions. - - if (.NOT.(PRESENT(nh_curr).OR.PRESENT(ng_curr)) ) then -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kme-1, kms, -1 - DO i=i_start(ij),i_end(ij) - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) - temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - endif - -! CASE (NSSL_3MOM) -! coming in future? - - CASE (CAMMGMPSCHEME) -! nothing to do here - - CASE (FULL_KHAIN_LYNN) -! explicit bin scheme so code below not applicable -! scheme authors need to implement if desired. - - CASE (FAST_KHAIN_LYNN) -! explicit bin scheme so code below not applicable -! scheme authors need to implement if desired. - - CASE DEFAULT - - END SELECT mp_select - - - if (scheme_has_graupel) then - -!..Create bins of graupel/hail (from 500 microns up to 7.5 cm). - xxDx(1) = 500.D-6 - xxDx(ngbins+1) = 0.075d0 - do n = 2, ngbins - xxDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(ngbins) & - *DLOG(xxDx(ngbins+1)/xxDx(1)) +DLOG(xxDx(1))) - enddo - do n = 1, ngbins - xxDg(n) = DSQRT(xxDx(n)*xxDx(n+1)) - xdtg(n) = xxDx(n+1) - xxDx(n) - enddo - - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO k=kms,kme-1 - DO i=i_start(ij),i_end(ij) - if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - lamg = (xam_g*cgg(3)/cgg(2)*temp_ng(i,k,j)/temp_qg(i,k,j))**(1./xbm_g) - N0_g = temp_ng(i,k,j)/cgg(2)*lamg**cge(2) - sum_ng = 0.0d0 - sum_t = 0.0d0 - do ng = ngbins, 1, -1 - f_d = N0_g*xxDg(ng)**xmu_g * DEXP(-lamg*xxDg(ng))*xdtg(ng) - sum_ng = sum_ng + f_d - if (sum_ng .gt. thresh_conc) then - exit - endif - sum_t = sum_ng - enddo - if (ng .ge. ngbins) then - hail_max = xxDg(ngbins) - elseif (xxDg(ng+1) .gt. 1.E-3) then - hail_max = xxDg(ng) - (sum_ng-thresh_conc)/(sum_ng-sum_t) & - & * (xxDg(ng)-xxDg(ng+1)) - else - hail_max = 1.E-4 - endif - if (hail_max .gt. 1E-2) then - WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000. - CALL wrf_debug (350, TRIM(outstring)) - endif - hail_max_sp = hail_max - if (k.eq.kms) then - hail_maxk1(i,j) = MAX(hail_maxk1(i,j), hail_max_sp) - endif - hail_max2d(i,j) = MAX(hail_max2d(i,j), hail_max_sp) - ENDDO - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - - endif - - - ENDIF -! NSSL - if (diag_print .eq. 0 ) return IF ( xtime .ne. 0. ) THEN @@ -1180,43 +510,5 @@ SUBROUTINE diagnostic_output_calc( & END SUBROUTINE diagnostic_output_calc - -!+---+-----------------------------------------------------------------+ - REAL FUNCTION GAMMLN(XX) -! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0. - IMPLICIT NONE - REAL, INTENT(IN):: XX - DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0 - DOUBLE PRECISION, DIMENSION(6), PARAMETER:: & - COF = (/76.18009172947146D0, -86.50532032941677D0, & - 24.01409824083091D0, -1.231739572450155D0, & - .1208650973866179D-2, -.5395239384953D-5/) - DOUBLE PRECISION:: SER,TMP,X,Y - INTEGER:: J - - X=XX - Y=X - TMP=X+5.5D0 - TMP=(X+0.5D0)*LOG(TMP)-TMP - SER=1.000000000190015D0 - DO 11 J=1,6 - Y=Y+1.D0 - SER=SER+COF(J)/Y -11 CONTINUE - GAMMLN=TMP+LOG(STP*SER/X) - END FUNCTION GAMMLN -! (C) Copr. 1986-92 Numerical Recipes Software 2.02 -!+---+-----------------------------------------------------------------+ - REAL FUNCTION WGAMMA(y) - - IMPLICIT NONE - REAL, INTENT(IN):: y - - WGAMMA = EXP(GAMMLN(y)) - - END FUNCTION WGAMMA -!+---+-----------------------------------------------------------------+ - - END MODULE module_diag_misc #endif diff --git a/phys/module_diag_nwp.F b/phys/module_diag_nwp.F new file mode 100644 index 0000000000..3d3c96fd0f --- /dev/null +++ b/phys/module_diag_nwp.F @@ -0,0 +1,835 @@ +#if ( NMM_CORE == 1) +MODULE module_diag_nwp +CONTAINS + SUBROUTINE diag_nwp_stub + END SUBROUTINE diag_nwp_stub +END MODULE module_diag_nwp +#else +!WRF:MEDIATION_LAYER:PHYSICS +! + +MODULE module_diag_nwp + PRIVATE :: WGAMMA + PRIVATE :: GAMMLN +CONTAINS + SUBROUTINE diagnostic_output_nwp( & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + ips,ipe, jps,jpe, kps,kpe, & ! patch dims + i_start,i_end,j_start,j_end,kts,kte,num_tiles & + ,u,v,temp,p8w & + ,dt,xtime,sbw & + ,mphysics_opt & + ,gsfcgce_hail, gsfcgce_2ice & + ,mpuse_hail & + ,nssl_cnoh, nssl_cnohl & + ,nssl_rho_qh, nssl_rho_qhl & + ,nssl_alphah, nssl_alphahl & + ,curr_secs2 & + ,nwp_diagnostics, diagflag & + ,history_interval & + ,itimestep & + ,u10,v10,w & + ,wspd10max & + ,up_heli_max & + ,w_up_max,w_dn_max & + ,znw,w_colmean & + ,numcolpts,w_mean & + ,grpl_max,grpl_colint,refd_max,refl_10cm & + ,hail_maxk1,hail_max2d & + ,qg_curr & + ,ng_curr,qh_curr,nh_curr,qr_curr,nr_curr & ! Optional (gthompsn) + ,rho,ph,phb,g & + ,max_time_step,adaptive_ts & + ) +!---------------------------------------------------------------------- + + USE module_state_description, ONLY : & + KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME, & + WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, & + MORR_TWO_MOMENT, GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, & + NSSL_2MOM, NSSL_2MOMG, NSSL_2MOMCCN, NSSL_1MOM, NSSL_1MOMLFO, & + MILBRANDT2MOM , CAMMGMPSCHEME, FAST_KHAIN_LYNN, FULL_KHAIN_LYNN, & + MORR_TM_AERO !TWG 2017 !,MILBRANDT3MOM, NSSL_3MOM + + IMPLICIT NONE +!====================================================================== +! Definitions +!----------- +!-- DIAG_PRINT print control: 0 - no diagnostics; 1 - dmudt only; 2 - all +!-- DT time step (second) +!-- XTIME forecast time +!-- SBW specified boundary width - used later +! +!-- P8W 3D pressure array at full eta levels +!-- U u component of wind - to be used later to compute k.e. +!-- V v component of wind - to be used later to compute k.e. +!-- CURR_SECS2 Time (s) since the beginning of the restart +!-- NWP_DIAGNOSTICS = 1, compute hourly maximum fields +!-- DIAGFLAG logical flag to indicate if this is a history output time +!-- U10, V10 10 m wind components +!-- WSPD10MAX 10 m max wind speed +!-- UP_HELI_MAX max updraft helicity +!-- W_UP_MAX max updraft vertical velocity +!-- W_DN_MAX max downdraft vertical velocity +!-- W_COLMEAN column mean vertical velocity +!-- NUMCOLPTS no of column points +!-- GRPL_MAX max column-integrated graupel +!-- GRPL_COLINT column-integrated graupel +!-- REF_MAX max derived radar reflectivity +!-- REFL_10CM model computed 3D reflectivity +! +!-- ids start index for i in domain +!-- ide end index for i in domain +!-- jds start index for j in domain +!-- jde end index for j in domain +!-- kds start index for k in domain +!-- kde end index for k in domain +!-- ims start index for i in memory +!-- ime end index for i in memory +!-- jms start index for j in memory +!-- jme end index for j in memory +!-- ips start index for i in patch +!-- ipe end index for i in patch +!-- jps start index for j in patch +!-- jpe end index for j in patch +!-- kms start index for k in memory +!-- kme end index for k in memory +!-- i_start start indices for i in tile +!-- i_end end indices for i in tile +!-- j_start start indices for j in tile +!-- j_end end indices for j in tile +!-- kts start index for k in tile +!-- kte end index for k in tile +!-- num_tiles number of tiles +! +!====================================================================== + + INTEGER, INTENT(IN ) :: & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + ips,ipe, jps,jpe, kps,kpe, & + kts,kte, & + num_tiles + + INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & + & i_start,i_end,j_start,j_end + + INTEGER, INTENT(IN ) :: mphysics_opt + INTEGER, INTENT(IN) :: gsfcgce_hail, gsfcgce_2ice, mpuse_hail + REAL, INTENT(IN) :: nssl_cnoh, nssl_cnohl & + ,nssl_rho_qh, nssl_rho_qhl & + ,nssl_alphah, nssl_alphahl + INTEGER, INTENT(IN ) :: nwp_diagnostics + LOGICAL, INTENT(IN ) :: diagflag + + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & + INTENT(IN ) :: u & + , v & + , p8w + + REAL, INTENT(IN ) :: DT, XTIME + INTEGER, INTENT(IN ) :: SBW + + REAL, OPTIONAL, INTENT(IN):: CURR_SECS2 + + INTEGER :: i,j,k,its,ite,jts,jte,ij + INTEGER :: idp,jdp,irc,jrc,irnc,jrnc,isnh,jsnh + + REAL :: no_points + REAL :: dpsdt_sum, dmudt_sum, dardt_sum, drcdt_sum, drndt_sum + REAL :: hfx_sum, lh_sum, sfcevp_sum, rainc_sum, rainnc_sum, raint_sum + REAL :: dmumax, raincmax, rainncmax, snowhmax + LOGICAL, EXTERNAL :: wrf_dm_on_monitor + CHARACTER*256 :: outstring + CHARACTER*6 :: grid_str + + INTEGER, INTENT(IN) :: & + history_interval,itimestep + + REAL, DIMENSION( kms:kme ), INTENT(IN) :: & + znw + + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: & + w & + ,temp & + ,qg_curr & + ,rho & + ,refl_10cm & + ,ph,phb + + REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: & + ng_curr, qh_curr, nh_curr & + ,qr_curr, nr_curr + + REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: & + u10 & + ,v10 + + REAL, INTENT(IN) :: g + + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: & + wspd10max & + ,up_heli_max & + ,w_up_max,w_dn_max & + ,w_colmean,numcolpts,w_mean & + ,grpl_max,grpl_colint & + ,hail_maxk1,hail_max2d & + ,refd_max + + REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: temp_qg, temp_ng, temp_qr, temp_nr + + INTEGER :: idump + + REAL :: wind_vel + REAL :: depth + + DOUBLE PRECISION:: hail_max + REAL:: hail_max_sp + DOUBLE PRECISION, PARAMETER:: thresh_conc = 0.0005d0 ! number conc. of graupel/hail per cubic meter + LOGICAL:: scheme_has_graupel + INTEGER, PARAMETER:: ngbins=50 + DOUBLE PRECISION, DIMENSION(ngbins+1):: xxDx + DOUBLE PRECISION, DIMENSION(ngbins):: xxDg, xdtg + REAL:: xrho_g, xam_g, xbm_g, xmu_g + REAL, DIMENSION(3):: cge, cgg + DOUBLE PRECISION:: f_d, sum_ng, sum_t, lamg, ilamg, N0_g, lam_exp, N0exp + DOUBLE PRECISION:: lamr, N0min + REAL:: mvd_r, xslw1, ygra1, zans1 + INTEGER:: ng, n + + REAL :: time_from_output + INTEGER :: max_time_step + LOGICAL :: adaptive_ts + LOGICAL :: reset_arrays = .FALSE. + +!----------------------------------------------------------------- + +! idump = (history_interval * 60.) / dt +! IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN +! IF (mod(curr_secs2, 60.*60.) == 0.) THEN +! WRITE(outstring,*) 'Resetting max refl array for domain with dt = ', dt +! CALL wrf_debug ( 0,TRIM(outstring) ) + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + refd_max(i,j) = 0. + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO +! ENDIF + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) +! Calculate the max radar reflectivity at this output times + IF ( refl_10cm(i,kms,j) .GT. refd_max(i,j) ) THEN + refd_max(i,j) = refl_10cm(i,kms,j) + ENDIF + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + + idump = (history_interval * 60.) / dt + +! IF ( MOD(itimestep, idump) .eq. 0 ) THEN +! WRITE(outstring,*) 'Computing PH0 for this domain with curr_secs2 = ', curr_secs2 +! CALL wrf_message ( TRIM(outstring) ) + + time_from_output = mod( xtime, REAL(history_interval) ) + +! print *,' history_interval = ', history_interval +! print *,' itimestep = ', itimestep +! print *,' idump = ', idump +! print *,' xtime = ', xtime +! print *,' time_from_output = ', time_from_output +! print *,' max_time_step = ', max_time_step + + IF ( adaptive_ts .EQV. .TRUE. ) THEN +! if timestep is adaptive, use new resetting method; +! also, we are multiplying max_time_step with 1.05; because of rounding error, +! the time_from_output can become slightly larger than max_time_step when +! adaptive time step is maximized, in which case the if condition below fails to detect +! that we just wrote an output + IF ( ( time_from_output .GT. 0 ) .AND. ( time_from_output .LE. ( ( max_time_step * 1.05 ) / 60. ) ) ) THEN + reset_arrays = .TRUE. + ENDIF + ELSE +! if timestep is fixed, use original resetting method + IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN + reset_arrays = .TRUE. + ENDIF + ENDIF + +! print *,' reset_arrays = ', reset_arrays + IF ( reset_arrays .EQV. .TRUE. ) THEN + +! print *,' Resetting NWP max arrays ' + +! IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN + WRITE(outstring,*) 'NSSL Diagnostics: Resetting max arrays for domain with dt = ', dt + CALL wrf_debug ( 10,TRIM(outstring) ) + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + wspd10max(i,j) = 0. + up_heli_max(i,j) = 0. + w_up_max(i,j) = 0. + w_dn_max(i,j) = 0. + w_mean(i,j) = 0. + grpl_max(i,j) = 0. + refd_max(i,j) = 0. + hail_maxk1(i,j) = 0. + hail_max2d(i,j) = 0. + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + reset_arrays = .FALSE. + ENDIF + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + +! Zero some accounting arrays that will be used below + + w_colmean(i,j) = 0. + numcolpts(i,j) = 0. + grpl_colint(i,j) = 0. + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme + DO i=i_start(ij),i_end(ij) + +! Find vertical velocity max (up and down) below 400 mb + + IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .GT. w_up_max(i,j) ) THEN + w_up_max(i,j) = w(i,k,j) + ENDIF + + IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .LT. w_dn_max(i,j) ) THEN + w_dn_max(i,j) = w(i,k,j) + ENDIF + +! For the column mean vertical velocity calculation, first +! total the vertical velocity between sigma levels 0.5 and 0.8 + + IF ( znw(k) .GE. 0.5 .AND. znw(k) .LE. 0.8 ) THEN + w_colmean(i,j) = w_colmean(i,j) + w(i,k,j) + numcolpts(i,j) = numcolpts(i,j) + 1 + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + +! Calculate the column integrated graupel + + depth = ( ( ph(i,k+1,j) + phb(i,k+1,j) ) / g ) - & + ( ( ph(i,k ,j) + phb(i,k ,j) ) / g ) + grpl_colint(i,j) = grpl_colint(i,j) + qg_curr(i,k,j) * depth * rho(i,k,j) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + +! Calculate the max 10 m wind speed between output times + + wind_vel = sqrt ( u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j) ) + IF ( wind_vel .GT. wspd10max(i,j) ) THEN + wspd10max(i,j) = wind_vel + ENDIF + +! Calculate the column mean vertical velocity between output times + + w_mean(i,j) = w_mean(i,j) + w_colmean(i,j) / numcolpts(i,j) + + IF ( MOD(itimestep, idump) .eq. 0 ) THEN + w_mean(i,j) = w_mean(i,j) / idump + ENDIF + +! Calculate the max column integrated graupel between output times + + IF ( grpl_colint(i,j) .gt. grpl_max(i,j) ) THEN + grpl_max(i,j) = grpl_colint(i,j) + ENDIF + + ! Calculate the max radar reflectivity between output times + + IF ( refl_10cm(i,kms,j) .GT. refd_max(i,j) ) THEN + refd_max(i,j) = refl_10cm(i,kms,j) + ENDIF + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ +!..Calculate a maximum hail diameter from the characteristics of the +!.. graupel category mixing ratio and number concentration (or hail, if +!.. available). This diagnostic uses the actual spectral distribution +!.. assumptions, calculated by breaking the distribution into 50 bins +!.. from 0.5mm to 7.5cm. Once a minimum number concentration of 0.01 +!.. particle per cubic meter of air is reached, from the upper size +!.. limit, then this bin is considered the max size. +!+---+-----------------------------------------------------------------+ + + WRITE(outstring,*) 'GT-Diagnostics, computing max-hail diameter' + CALL wrf_debug (100, TRIM(outstring)) + + IF (PRESENT(qh_curr)) THEN + WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail mixing ratio' + CALL wrf_debug (150, TRIM(outstring)) +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + temp_qg(i,k,j) = MAX(1.E-12, qh_curr(i,k,j)*rho(i,k,j)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + ELSE +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + temp_qg(i,k,j) = MAX(1.E-12, qg_curr(i,k,j)*rho(i,k,j)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + ENDIF + + IF (PRESENT(nh_curr)) THEN + WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail number concentration' + CALL wrf_debug (150, TRIM(outstring)) +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + temp_ng(i,k,j) = MAX(1.E-8, nh_curr(i,k,j)*rho(i,k,j)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + ELSEIF (PRESENT(ng_curr)) THEN + WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel number concentration' +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + temp_ng(i,k,j) = MAX(1.E-8, ng_curr(i,k,j)*rho(i,k,j)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + ELSE +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + temp_ng(i,k,j) = 1.E-8 + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + ENDIF + + ! Calculate the max hail size from graupel/hail parameters in microphysics scheme (gthompsn 12Mar2015) + ! First, we do know multiple schemes have assumed inverse-exponential distribution with constant + ! intercept parameter and particle density. Please leave next 4 settings alone for common + ! use and copy these and cge constants to re-assign per scheme if needed (e.g. NSSL). + + xrho_g = 500. + xam_g = 3.1415926536/6.0*xrho_g ! Assumed m(D) = a*D**b, where a=PI/6*rho_g and b=3 + xbm_g = 3. ! in other words, spherical graupel/hail + xmu_g = 0. + scheme_has_graupel = .false. + + !..Some constants. These *MUST* get changed below per scheme + !.. *IF* either xbm_g or xmu_g value is changed from 3 and zero, respectively. + + cge(1) = xbm_g + 1. + cge(2) = xmu_g + 1. + cge(3) = xbm_g + xmu_g + 1. + do n = 1, 3 + cgg(n) = WGAMMA(cge(n)) + enddo + + mp_select: SELECT CASE(mphysics_opt) + + CASE (KESSLERSCHEME) +! nothing to do here + + CASE (LINSCHEME) + scheme_has_graupel = .true. + xrho_g = 917. + xam_g = 3.1415926536/6.0*xrho_g + N0exp = 4.e4 +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kme-1, kms, -1 + DO i=i_start(ij),i_end(ij) + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) + temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + + CASE (WSM3SCHEME) +! nothing to do here + + CASE (WSM5SCHEME) +! nothing to do here + + CASE (WSM6SCHEME) + scheme_has_graupel = .true. + xrho_g = 500. + xam_g = 3.1415926536/6.0*xrho_g + N0exp = 4.e6 +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kme-1, kms, -1 + DO i=i_start(ij),i_end(ij) + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) + temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + + CASE (WDM5SCHEME) +! nothing to do here + + CASE (WDM6SCHEME) + scheme_has_graupel = .true. + xrho_g = 500. + N0exp = 4.e6 + if (mpuse_hail .eq. 1) then + xrho_g = 700. + N0exp = 4.e4 + endif + xam_g = 3.1415926536/6.0*xrho_g +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kme-1, kms, -1 + DO i=i_start(ij),i_end(ij) + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) + temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + + CASE (GSFCGCESCHEME) + if (gsfcgce_2ice.eq.0 .OR. gsfcgce_2ice.eq.2) then + scheme_has_graupel = .true. + if (gsfcgce_hail .eq. 1) then + xrho_g = 900. + else + xrho_g = 400. + endif + xam_g = 3.1415926536/6.0*xrho_g + N0exp = 4.e4 +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kme-1, kms, -1 + DO i=i_start(ij),i_end(ij) + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) + temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + endif + + CASE (SBU_YLINSCHEME) +! scheme_has_graupel = .true. ! Can be calculated from rime fraction variable. +! no time to implement + + CASE (ETAMPNEW) +! scheme_has_graupel = .true. ! Can be calculated from rime fraction variable. +! no time to implement + + CASE (THOMPSON, THOMPSONAERO) + + scheme_has_graupel = .true. + xmu_g = 1. + cge(1) = xbm_g + 1. + cge(2) = xmu_g + 1. + cge(3) = xbm_g + xmu_g + 1. + do n = 1, 3 + cgg(n) = WGAMMA(cge(n)) + enddo + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + DO k=kme-1, kms, -1 + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + zans1 = (2.5 + 2./7. * (ALOG10(temp_qg(i,k,j))+7.)) + zans1 = MAX(2., MIN(zans1, 7.)) + N0exp = 10.**zans1 + lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) + lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g) + N0_g = N0exp/(cgg(2)*lam_exp) * lamg**cge(2) + temp_ng(i,k,j) = N0_g*cgg(2)*lamg**(-cge(2)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + + +! CASE (MORR_MILB_P3) +! scheme_has_graupel = .true. +! either Hugh or Jason need to implement. + + CASE (MORR_TWO_MOMENT, MORR_TM_AERO) + scheme_has_graupel = .true. + xrho_g = 400. + if (mpuse_hail .eq. 1) xrho_g = 900. + xam_g = 3.1415926536/6.0*xrho_g + + CASE (MILBRANDT2MOM) + WRITE(outstring,*) 'GT-Debug, using Milbrandt2mom, which has 2-moment hail' + CALL wrf_debug (150, TRIM(outstring)) + scheme_has_graupel = .true. + xrho_g = 900. + xam_g = 3.1415926536/6.0*xrho_g + +! CASE (MILBRANDT3MOM) +! coming in future? + + CASE (NSSL_1MOMLFO, NSSL_1MOM, NSSL_2MOM, NSSL_2MOMG, NSSL_2MOMCCN) + + scheme_has_graupel = .true. + xrho_g = nssl_rho_qh + N0exp = nssl_cnoh + if (PRESENT(qh_curr)) then + xrho_g = nssl_rho_qhl + N0exp = nssl_cnohl + endif + xam_g = 3.1415926536/6.0*xrho_g + + if (PRESENT(ng_curr)) xmu_g = nssl_alphah + if (PRESENT(nh_curr)) xmu_g = nssl_alphahl + + if (xmu_g .NE. 0.) then + cge(1) = xbm_g + 1. + cge(2) = xmu_g + 1. + cge(3) = xbm_g + xmu_g + 1. + do n = 1, 3 + cgg(n) = WGAMMA(cge(n)) + enddo + endif + + ! NSSL scheme has many options, but, if single-moment, just fill + ! in the number array for graupel from built-in assumptions. + + if (.NOT.(PRESENT(nh_curr).OR.PRESENT(ng_curr)) ) then +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kme-1, kms, -1 + DO i=i_start(ij),i_end(ij) + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) + temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + endif + +! CASE (NSSL_3MOM) +! coming in future? + + CASE (CAMMGMPSCHEME) +! nothing to do here + + CASE (FULL_KHAIN_LYNN) +! explicit bin scheme so code below not applicable +! scheme authors need to implement if desired. + + CASE (FAST_KHAIN_LYNN) +! explicit bin scheme so code below not applicable +! scheme authors need to implement if desired. + + CASE DEFAULT + + END SELECT mp_select + + + if (scheme_has_graupel) then + +!..Create bins of graupel/hail (from 500 microns up to 7.5 cm). + xxDx(1) = 500.D-6 + xxDx(ngbins+1) = 0.075d0 + do n = 2, ngbins + xxDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(ngbins) & + *DLOG(xxDx(ngbins+1)/xxDx(1)) +DLOG(xxDx(1))) + enddo + do n = 1, ngbins + xxDg(n) = DSQRT(xxDx(n)*xxDx(n+1)) + xdtg(n) = xxDx(n+1) - xxDx(n) + enddo + + +! !$OMP PARALLEL DO & +! !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + DO j=j_start(ij),j_end(ij) + DO k=kms,kme-1 + DO i=i_start(ij),i_end(ij) + if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE + lamg = (xam_g*cgg(3)/cgg(2)*temp_ng(i,k,j)/temp_qg(i,k,j))**(1./xbm_g) + N0_g = temp_ng(i,k,j)/cgg(2)*lamg**cge(2) + sum_ng = 0.0d0 + sum_t = 0.0d0 + do ng = ngbins, 1, -1 + f_d = N0_g*xxDg(ng)**xmu_g * DEXP(-lamg*xxDg(ng))*xdtg(ng) + sum_ng = sum_ng + f_d + if (sum_ng .gt. thresh_conc) then + exit + endif + sum_t = sum_ng + enddo + if (ng .ge. ngbins) then + hail_max = xxDg(ngbins) + elseif (xxDg(ng+1) .gt. 1.E-3) then + hail_max = xxDg(ng) - (sum_ng-thresh_conc)/(sum_ng-sum_t) & + & * (xxDg(ng)-xxDg(ng+1)) + else + hail_max = 1.E-4 + endif + if (hail_max .gt. 1E-2) then + WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000. + CALL wrf_debug (350, TRIM(outstring)) + endif + hail_max_sp = hail_max + if (k.eq.kms) then + hail_maxk1(i,j) = MAX(hail_maxk1(i,j), hail_max_sp) + endif + hail_max2d(i,j) = MAX(hail_max2d(i,j), hail_max_sp) + ENDDO + ENDDO + ENDDO + ENDDO +! !$OMP END PARALLEL DO + + endif + + END SUBROUTINE diagnostic_output_nwp + +!+---+-----------------------------------------------------------------+ + REAL FUNCTION GAMMLN(XX) +! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0. + IMPLICIT NONE + REAL, INTENT(IN):: XX + DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0 + DOUBLE PRECISION, DIMENSION(6), PARAMETER:: & + COF = (/76.18009172947146D0, -86.50532032941677D0, & + 24.01409824083091D0, -1.231739572450155D0, & + .1208650973866179D-2, -.5395239384953D-5/) + DOUBLE PRECISION:: SER,TMP,X,Y + INTEGER:: J + + X=XX + Y=X + TMP=X+5.5D0 + TMP=(X+0.5D0)*LOG(TMP)-TMP + SER=1.000000000190015D0 + DO 11 J=1,6 + Y=Y+1.D0 + SER=SER+COF(J)/Y +11 CONTINUE + GAMMLN=TMP+LOG(STP*SER/X) + END FUNCTION GAMMLN +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION WGAMMA(y) + + IMPLICIT NONE + REAL, INTENT(IN):: y + + WGAMMA = EXP(GAMMLN(y)) + + END FUNCTION WGAMMA +!+---+-----------------------------------------------------------------+ + +END MODULE module_diag_nwp +#endif diff --git a/phys/module_diagnostics_driver.F b/phys/module_diagnostics_driver.F index 867c5dae71..bee460fead 100644 --- a/phys/module_diagnostics_driver.F +++ b/phys/module_diagnostics_driver.F @@ -74,6 +74,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & USE module_lightning_driver, ONLY : lightning_driver USE module_diag_misc, ONLY : diagnostic_output_calc + USE module_diag_nwp, ONLY : diagnostic_output_nwp USE module_diag_cl, ONLY : clwrf_output_calc USE module_diag_pld, ONLY : pld USE module_diag_zld, ONLY : zld @@ -337,29 +338,21 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & !$OMP END PARALLEL DO END IF TRADITIONAL_FIELDS - ! Mostly surface values, precip, column integrated quantities. - CALL wrf_debug ( 100 , '--> CALL DIAGNOSTICS PACKAGE: NWP DIAGNOSTICS' ) - - - - mp_select: SELECT CASE(config_flags%mp_physics) - - CASE (LINSCHEME, WSM6SCHEME, WDM6SCHEME, GSFCGCESCHEME, NSSL_1MOMLFO) + CALL wrf_debug ( 100 , '--> CALL DIAGNOSTICS PACKAGE: ACCUMULATED AND BUCKET DIAGNOSTICS' ) CALL diagnostic_output_calc( & DPSDT=grid%dpsdt ,DMUDT=grid%dmudt & ,P8W=p8w ,PK1M=grid%pk1m & ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m & - ,U=grid%u_2 ,V=grid%v_2 & - ,TEMP=grid%t_phy & ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv & ,RAINC=grid%rainc ,RAINNC=grid%rainnc & ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc & ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh & - ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & - ,XTIME=grid%xtime ,T2=grid%t2 & + ,T2=grid%t2 & + ,DT=grid%dt & + ,XTIME=grid%xtime & ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc & @@ -380,6 +373,36 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,DIAG_PRINT=config_flags%diag_print & ,BUCKET_MM=config_flags%bucket_mm & ,BUCKET_J =config_flags%bucket_J & + ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & + ,PREC_ACC_C=grid%prec_acc_c & + ,PREC_ACC_NC=grid%prec_acc_nc & + ,PREC_ACC_DT=config_flags%prec_acc_dt & + ,CURR_SECS2=curr_secs2 & + ,HISTORY_INTERVAL=grid%history_interval & + ,ITIMESTEP=grid%itimestep & + ! Dimension arguments + ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & + ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & + ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe & + ,I_START=grid%i_start,I_END=min(grid%i_end, ide-1) & + ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) & + ,KTS=k_start, KTE=min(k_end,kde-1) & + ,NUM_TILES=grid%num_tiles & + ) + + NWPDIAGS: IF ( config_flags%nwp_diagnostics == 1 ) THEN + CALL wrf_debug ( 100 , '--> CALL DIAGNOSTICS PACKAGE: NWP DIAGNOSTICS' ) + + mp_select: SELECT CASE(config_flags%mp_physics) + + CASE (LINSCHEME, WSM6SCHEME, WDM6SCHEME, GSFCGCESCHEME, NSSL_1MOMLFO) + + CALL diagnostic_output_nwp( & + U=grid%u_2 ,V=grid%v_2 & + ,TEMP=grid%t_phy ,P8W=p8w & + ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & + ,XTIME=grid%xtime & + ! Selection flag ,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn ,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn ,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn @@ -390,10 +413,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn ,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn ,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn - ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & - ,PREC_ACC_C=grid%prec_acc_c & - ,PREC_ACC_NC=grid%prec_acc_nc & - ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & @@ -420,43 +439,17 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,KTS=k_start, KTE=min(k_end,kde-1) & ,NUM_TILES=grid%num_tiles & ,MAX_TIME_STEP=grid%max_time_step & - ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & + ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & ) CASE (THOMPSON, THOMPSONAERO) - CALL diagnostic_output_calc( & - DPSDT=grid%dpsdt ,DMUDT=grid%dmudt & - ,P8W=p8w ,PK1M=grid%pk1m & - ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m & - ,U=grid%u_2 ,V=grid%v_2 & - ,TEMP=grid%t_phy & - ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv & - ,RAINC=grid%rainc ,RAINNC=grid%rainnc & - ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc & - ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh & + CALL diagnostic_output_nwp( & + U=grid%u_2 ,V=grid%v_2 & + ,TEMP=grid%t_phy ,P8W=p8w & ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & - ,XTIME=grid%xtime ,T2=grid%t2 & - ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & - ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & - ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc & - ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc & - ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc & - ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc & - ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc & - ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc & - ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc & - ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc & - ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc & - ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc & - ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc & - ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc & - ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc & - ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc & + ,XTIME=grid%xtime & ! Selection flag - ,DIAG_PRINT=config_flags%diag_print & - ,BUCKET_MM=config_flags%bucket_mm & - ,BUCKET_J =config_flags%bucket_J & ,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn ,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn ,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn @@ -467,10 +460,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn ,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn ,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn - ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & - ,PREC_ACC_C=grid%prec_acc_c & - ,PREC_ACC_NC=grid%prec_acc_nc & - ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & @@ -498,44 +487,18 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) & ,KTS=k_start, KTE=min(k_end,kde-1) & ,NUM_TILES=grid%num_tiles & - ,MAX_TIME_STEP=grid%max_time_step & - ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & + ,MAX_TIME_STEP=grid%max_time_step & + ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & ) CASE (MORR_TWO_MOMENT, MORR_TM_AERO) ! TWG add - CALL diagnostic_output_calc( & - DPSDT=grid%dpsdt ,DMUDT=grid%dmudt & - ,P8W=p8w ,PK1M=grid%pk1m & - ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m & - ,U=grid%u_2 ,V=grid%v_2 & - ,TEMP=grid%t_phy & - ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv & - ,RAINC=grid%rainc ,RAINNC=grid%rainnc & - ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc & - ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh & + CALL diagnostic_output_nwp( & + U=grid%u_2 ,V=grid%v_2 & + ,TEMP=grid%t_phy ,P8W=p8w & ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & - ,XTIME=grid%xtime ,T2=grid%t2 & - ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & - ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & - ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc & - ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc & - ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc & - ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc & - ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc & - ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc & - ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc & - ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc & - ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc & - ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc & - ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc & - ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc & - ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc & - ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc & + ,XTIME=grid%xtime & ! Selection flag - ,DIAG_PRINT=config_flags%diag_print & - ,BUCKET_MM=config_flags%bucket_mm & - ,BUCKET_J =config_flags%bucket_J & ,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn ,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn ,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn @@ -546,10 +509,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn ,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn ,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn - ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & - ,PREC_ACC_C=grid%prec_acc_c & - ,PREC_ACC_NC=grid%prec_acc_nc & - ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & @@ -576,44 +535,18 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) & ,KTS=k_start, KTE=min(k_end,kde-1) & ,NUM_TILES=grid%num_tiles & - ,MAX_TIME_STEP=grid%max_time_step & - ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & + ,MAX_TIME_STEP=grid%max_time_step & + ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & ) CASE (NSSL_1MOM) - CALL diagnostic_output_calc( & - DPSDT=grid%dpsdt ,DMUDT=grid%dmudt & - ,P8W=p8w ,PK1M=grid%pk1m & - ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m & - ,U=grid%u_2 ,V=grid%v_2 & - ,TEMP=grid%t_phy & - ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv & - ,RAINC=grid%rainc ,RAINNC=grid%rainnc & - ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc & - ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh & + CALL diagnostic_output_nwp( & + U=grid%u_2 ,V=grid%v_2 & + ,TEMP=grid%t_phy ,P8W=p8w & ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & - ,XTIME=grid%xtime ,T2=grid%t2 & - ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & - ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & - ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc & - ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc & - ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc & - ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc & - ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc & - ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc & - ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc & - ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc & - ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc & - ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc & - ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc & - ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc & - ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc & - ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc & + ,XTIME=grid%xtime & ! Selection flag - ,DIAG_PRINT=config_flags%diag_print & - ,BUCKET_MM=config_flags%bucket_mm & - ,BUCKET_J =config_flags%bucket_J & ,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn ,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn ,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn @@ -624,10 +557,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn ,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn ,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn - ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & - ,PREC_ACC_C=grid%prec_acc_c & - ,PREC_ACC_NC=grid%prec_acc_nc & - ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & @@ -654,44 +583,18 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) & ,KTS=k_start, KTE=min(k_end,kde-1) & ,NUM_TILES=grid%num_tiles & - ,MAX_TIME_STEP=grid%max_time_step & - ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & + ,MAX_TIME_STEP=grid%max_time_step & + ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & ) CASE (MILBRANDT2MOM, NSSL_2MOM, NSSL_2MOMCCN) - CALL diagnostic_output_calc( & - DPSDT=grid%dpsdt ,DMUDT=grid%dmudt & - ,P8W=p8w ,PK1M=grid%pk1m & - ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m & - ,U=grid%u_2 ,V=grid%v_2 & - ,TEMP=grid%t_phy & - ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv & - ,RAINC=grid%rainc ,RAINNC=grid%rainnc & - ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc & - ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh & + CALL diagnostic_output_nwp( & + U=grid%u_2 ,V=grid%v_2 & + ,TEMP=grid%t_phy ,P8W=p8w & ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & - ,XTIME=grid%xtime ,T2=grid%t2 & - ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & - ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & - ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc & - ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc & - ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc & - ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc & - ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc & - ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc & - ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc & - ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc & - ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc & - ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc & - ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc & - ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc & - ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc & - ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc & + ,XTIME=grid%xtime & ! Selection flag - ,DIAG_PRINT=config_flags%diag_print & - ,BUCKET_MM=config_flags%bucket_mm & - ,BUCKET_J =config_flags%bucket_J & ,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn ,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn ,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn @@ -702,10 +605,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn ,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn ,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn - ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & - ,PREC_ACC_C=grid%prec_acc_c & - ,PREC_ACC_NC=grid%prec_acc_nc & - ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & @@ -734,11 +633,10 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) & ,KTS=k_start, KTE=min(k_end,kde-1) & ,NUM_TILES=grid%num_tiles & - ,MAX_TIME_STEP=grid%max_time_step & - ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & + ,MAX_TIME_STEP=grid%max_time_step & + ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & ) - !..The remaining microphysics schemes do not have graupel, but !..P_QG will just be empty and the remaining NWP-diagnostics can !..still be computed, so go ahead, under DEFAULT, not their own. @@ -769,38 +667,12 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & CASE DEFAULT - CALL diagnostic_output_calc( & - DPSDT=grid%dpsdt ,DMUDT=grid%dmudt & - ,P8W=p8w ,PK1M=grid%pk1m & - ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m & - ,U=grid%u_2 ,V=grid%v_2 & - ,TEMP=grid%t_phy & - ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv & - ,RAINC=grid%rainc ,RAINNC=grid%rainnc & - ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc & - ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh & + CALL diagnostic_output_nwp( & + U=grid%u_2 ,V=grid%v_2 & + ,TEMP=grid%t_phy ,P8W=p8w & ,DT=grid%dt ,SBW=config_flags%spec_bdy_width & - ,XTIME=grid%xtime ,T2=grid%t2 & - ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & - ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & - ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc & - ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc & - ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc & - ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc & - ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc & - ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc & - ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc & - ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc & - ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc & - ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc & - ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc & - ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc & - ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc & - ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc & + ,XTIME=grid%xtime & ! Selection flag - ,DIAG_PRINT=config_flags%diag_print & - ,BUCKET_MM=config_flags%bucket_mm & - ,BUCKET_J =config_flags%bucket_J & ,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn ,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn ,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn @@ -811,10 +683,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn ,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn ,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn - ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc & - ,PREC_ACC_C=grid%prec_acc_c & - ,PREC_ACC_NC=grid%prec_acc_nc & - ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & @@ -840,13 +708,12 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) & ,KTS=k_start, KTE=min(k_end,kde-1) & ,NUM_TILES=grid%num_tiles & - ,MAX_TIME_STEP=grid%max_time_step & - ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & + ,MAX_TIME_STEP=grid%max_time_step & + ,ADAPTIVE_TS=config_flags%use_adaptive_time_step & ) - - END SELECT mp_select + END IF NWPDIAGS ! Climate-oriented diagnostic quantities. From 701f82647ce5cf98bb7073fdc98ba1a757245256 Mon Sep 17 00:00:00 2001 From: Wei Wang Date: Thu, 13 Feb 2020 19:25:11 -0700 Subject: [PATCH 2/3] add dependency because of new file --- main/depend.common | 1 + 1 file changed, 1 insertion(+) diff --git a/main/depend.common b/main/depend.common index 1dd9e79b2b..143d66de00 100644 --- a/main/depend.common +++ b/main/depend.common @@ -760,6 +760,7 @@ module_sf_ocean_driver.o : \ module_diagnostics_driver.o: \ module_lightning_driver.o \ module_diag_misc.o \ + module_diag_nwp.o \ module_diag_cl.o \ module_diag_pld.o \ module_diag_zld.o \ From 67b856e00a05e7c04a10bde615ad907a781a3fc3 Mon Sep 17 00:00:00 2001 From: Wei Wang Date: Fri, 14 Feb 2020 10:47:33 -0700 Subject: [PATCH 3/3] remove refd_max calculation near the beginning, which is not part of the original code --- phys/module_diag_nwp.F | 32 -------------------------------- 1 file changed, 32 deletions(-) diff --git a/phys/module_diag_nwp.F b/phys/module_diag_nwp.F index 3d3c96fd0f..d15d8698ff 100644 --- a/phys/module_diag_nwp.F +++ b/phys/module_diag_nwp.F @@ -205,38 +205,6 @@ SUBROUTINE diagnostic_output_nwp( & !----------------------------------------------------------------- -! idump = (history_interval * 60.) / dt -! IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN -! IF (mod(curr_secs2, 60.*60.) == 0.) THEN -! WRITE(outstring,*) 'Resetting max refl array for domain with dt = ', dt -! CALL wrf_debug ( 0,TRIM(outstring) ) - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO i=i_start(ij),i_end(ij) - refd_max(i,j) = 0. - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO -! ENDIF - -! !$OMP PARALLEL DO & -! !$OMP PRIVATE ( ij ) - DO ij = 1 , num_tiles - DO j=j_start(ij),j_end(ij) - DO i=i_start(ij),i_end(ij) -! Calculate the max radar reflectivity at this output times - IF ( refl_10cm(i,kms,j) .GT. refd_max(i,j) ) THEN - refd_max(i,j) = refl_10cm(i,kms,j) - ENDIF - ENDDO - ENDDO - ENDDO -! !$OMP END PARALLEL DO - idump = (history_interval * 60.) / dt ! IF ( MOD(itimestep, idump) .eq. 0 ) THEN