#if ( NMM_CORE == 1)
MODULE module_diag_misc 1
CONTAINS
SUBROUTINE diag_misc_stub
END SUBROUTINE diag_misc_stub
END MODULE module_diag_misc
#else
!WRF:MEDIATION_LAYER:PHYSICS
!
MODULE module_diag_misc 1
PRIVATE :: WGAMMA
PRIVATE :: GAMMLN
CONTAINS
SUBROUTINE diagnostic_output_calc( & 6,22
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 &
,dpsdt,dmudt &
,p8w,pk1m,mu_2,mu_2m &
,u,v, temp &
,raincv,rainncv,rainc,rainnc &
,i_rainc,i_rainnc &
,hfx,sfcevp,lh &
,ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC & ! Optional
,ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC & ! Optional
,ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC & ! Optional
,ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC & ! Optional
,I_ACSWUPT,I_ACSWUPTC,I_ACSWDNT,I_ACSWDNTC & ! Optional
,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 &
,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 &
)
!----------------------------------------------------------------------
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 !,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
!-- MU dry column hydrostatic pressure
!-- RAINC cumulus scheme precipitation since hour 0
!-- RAINCV cumulus scheme precipitation in one time step (mm)
!-- RAINNC explicit scheme precipitation since hour 0
!-- RAINNCV explicit scheme precipitation in one time step (mm)
!-- SNOWNCV explicit scheme snow in one time step (mm)
!-- 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
!-- 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 ) :: 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
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: &
MU_2 &
, RAINNCV &
, RAINCV &
, SNOWNCV &
, HFX &
, LH &
, SFCEVP &
, T2
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: DPSDT &
, DMUDT &
, RAINNC &
, RAINC &
, MU_2M &
, PK1M
REAL, INTENT(IN ) :: DT, XTIME
INTEGER, INTENT(IN ) :: SBW
INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: &
I_RAINC, &
I_RAINNC
REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, &
ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, &
ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, &
ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC
INTEGER, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
I_ACSWUPT,I_ACSWUPTC,I_ACSWDNT,I_ACSWDNTC, &
I_ACSWUPB,I_ACSWUPBC,I_ACSWDNB,I_ACSWDNBC, &
I_ACLWUPT,I_ACLWUPTC,I_ACLWDNT,I_ACLWDNTC, &
I_ACLWUPB,I_ACLWUPBC,I_ACLWDNB,I_ACLWDNBC
REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
PREC_ACC_C, PREC_ACC_NC, SNOW_ACC_NC
REAL, OPTIONAL, INTENT(IN):: PREC_ACC_DT, 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.001d0 ! 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, lamg, ilamg, N0_g, lam_exp, N0exp
DOUBLE PRECISION:: lamr, N0min
REAL:: mvd_r, xslw1, ygra1, zans1
INTEGER:: ng, n
!-----------------------------------------------------------------
! 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
IF(bucket_mm .gt. 0. .AND. MOD(NINT(XTIME),360) .EQ. 0)THEN
! SET START AND END POINTS FOR TILES
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
DO ij = 1 , num_tiles
IF (xtime .eq. 0.0)THEN
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
i_rainnc(i,j) = 0
i_rainc(i,j) = 0
ENDDO
ENDDO
ENDIF
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
IF(rainnc(i,j) .gt. bucket_mm)THEN
rainnc(i,j) = rainnc(i,j) - bucket_mm
i_rainnc(i,j) = i_rainnc(i,j) + 1
ENDIF
IF(rainc(i,j) .gt. bucket_mm)THEN
rainc(i,j) = rainc(i,j) - bucket_mm
i_rainc(i,j) = i_rainc(i,j) + 1
ENDIF
ENDDO
ENDDO
IF (xtime .eq. 0.0 .and. bucket_J .gt. 0.0 .and. PRESENT(ACSWUPT))THEN
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
i_acswupt(i,j) = 0
i_acswuptc(i,j) = 0
i_acswdnt(i,j) = 0
i_acswdntc(i,j) = 0
i_acswupb(i,j) = 0
i_acswupbc(i,j) = 0
i_acswdnb(i,j) = 0
i_acswdnbc(i,j) = 0
ENDDO
ENDDO
ENDIF
IF (xtime .eq. 0.0 .and. bucket_J .gt. 0.0 .and. PRESENT(ACLWUPT))THEN
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
i_aclwupt(i,j) = 0
i_aclwuptc(i,j) = 0
i_aclwdnt(i,j) = 0
i_aclwdntc(i,j) = 0
i_aclwupb(i,j) = 0
i_aclwupbc(i,j) = 0
i_aclwdnb(i,j) = 0
i_aclwdnbc(i,j) = 0
ENDDO
ENDDO
ENDIF
IF (PRESENT(ACSWUPT) .and. bucket_J .gt. 0.0)THEN
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
IF(acswupt(i,j) .gt. bucket_J)THEN
acswupt(i,j) = acswupt(i,j) - bucket_J
i_acswupt(i,j) = i_acswupt(i,j) + 1
ENDIF
IF(acswuptc(i,j) .gt. bucket_J)THEN
acswuptc(i,j) = acswuptc(i,j) - bucket_J
i_acswuptc(i,j) = i_acswuptc(i,j) + 1
ENDIF
IF(acswdnt(i,j) .gt. bucket_J)THEN
acswdnt(i,j) = acswdnt(i,j) - bucket_J
i_acswdnt(i,j) = i_acswdnt(i,j) + 1
ENDIF
IF(acswdntc(i,j) .gt. bucket_J)THEN
acswdntc(i,j) = acswdntc(i,j) - bucket_J
i_acswdntc(i,j) = i_acswdntc(i,j) + 1
ENDIF
IF(acswupb(i,j) .gt. bucket_J)THEN
acswupb(i,j) = acswupb(i,j) - bucket_J
i_acswupb(i,j) = i_acswupb(i,j) + 1
ENDIF
IF(acswupbc(i,j) .gt. bucket_J)THEN
acswupbc(i,j) = acswupbc(i,j) - bucket_J
i_acswupbc(i,j) = i_acswupbc(i,j) + 1
ENDIF
IF(acswdnb(i,j) .gt. bucket_J)THEN
acswdnb(i,j) = acswdnb(i,j) - bucket_J
i_acswdnb(i,j) = i_acswdnb(i,j) + 1
ENDIF
IF(acswdnbc(i,j) .gt. bucket_J)THEN
acswdnbc(i,j) = acswdnbc(i,j) - bucket_J
i_acswdnbc(i,j) = i_acswdnbc(i,j) + 1
ENDIF
ENDDO
ENDDO
ENDIF
IF (PRESENT(ACLWUPT) .and. bucket_J .gt. 0.0)THEN
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
IF(aclwupt(i,j) .gt. bucket_J)THEN
aclwupt(i,j) = aclwupt(i,j) - bucket_J
i_aclwupt(i,j) = i_aclwupt(i,j) + 1
ENDIF
IF(aclwuptc(i,j) .gt. bucket_J)THEN
aclwuptc(i,j) = aclwuptc(i,j) - bucket_J
i_aclwuptc(i,j) = i_aclwuptc(i,j) + 1
ENDIF
IF(aclwdnt(i,j) .gt. bucket_J)THEN
aclwdnt(i,j) = aclwdnt(i,j) - bucket_J
i_aclwdnt(i,j) = i_aclwdnt(i,j) + 1
ENDIF
IF(aclwdntc(i,j) .gt. bucket_J)THEN
aclwdntc(i,j) = aclwdntc(i,j) - bucket_J
i_aclwdntc(i,j) = i_aclwdntc(i,j) + 1
ENDIF
IF(aclwupb(i,j) .gt. bucket_J)THEN
aclwupb(i,j) = aclwupb(i,j) - bucket_J
i_aclwupb(i,j) = i_aclwupb(i,j) + 1
ENDIF
IF(aclwupbc(i,j) .gt. bucket_J)THEN
aclwupbc(i,j) = aclwupbc(i,j) - bucket_J
i_aclwupbc(i,j) = i_aclwupbc(i,j) + 1
ENDIF
IF(aclwdnb(i,j) .gt. bucket_J)THEN
aclwdnb(i,j) = aclwdnb(i,j) - bucket_J
i_aclwdnb(i,j) = i_aclwdnb(i,j) + 1
ENDIF
IF(aclwdnbc(i,j) .gt. bucket_J)THEN
aclwdnbc(i,j) = aclwdnbc(i,j) - bucket_J
i_aclwdnbc(i,j) = i_aclwdnbc(i,j) + 1
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
! !$OMP END PARALLEL DO
ENDIF
! Compute precipitation accumulation in a given time window: prec_acc_dt
IF (prec_acc_dt .gt. 0.) THEN
! !$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)
IF (mod(curr_secs2, 60.* prec_acc_dt) == 0.) THEN
prec_acc_c(i,j) = 0.
prec_acc_nc(i,j) = 0.
snow_acc_nc(i,j) = 0.
ENDIF
prec_acc_c(i,j) = prec_acc_c(i,j) + RAINCV(i,j)
prec_acc_nc(i,j) = prec_acc_nc(i,j) + RAINNCV(i,j)
prec_acc_c(i,j) = MAX (prec_acc_c(i,j), 0.0)
prec_acc_nc(i,j) = MAX (prec_acc_nc(i,j), 0.0)
snow_acc_nc(i,j) = snow_acc_nc(i,j) + SNOWNCV(I,J)
! add convective precip to snow bucket if t2 < 273.15
IF ( t2(i,j) .lt. 273.15 ) THEN
snow_acc_nc(i,j) = snow_acc_nc(i,j) + RAINCV(i,j)
snow_acc_nc(i,j) = MAX (snow_acc_nc(i,j), 0.0)
ENDIF
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
ENDIF
! NSSL
IF ( nwp_diagnostics .EQ. 1 ) THEN
idump = (history_interval * 60.) / dt
! print *,' history_interval = ', history_interval
! print *,' itimestep = ', itimestep
! print *,' idump = ', idump
! print *,' xtime = ', xtime
! IF ( MOD(itimestep, idump) .eq. 0 ) THEN
! WRITE(outstring,*) 'Computing PH0 for this domain with curr_secs2 = ', curr_secs2
! CALL wrf_message ( TRIM(outstring) )
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
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 k=kms,kme-1
DO i=i_start(ij),i_end(ij)
temp_qr(i,k,j) = MAX(1.E-10, qr_curr(i,k,j)*rho(i,k,j))
temp_nr(i,k,j) = MAX(1.E-8, nr_curr(i,k,j)*rho(i,k,j))
ENDDO
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
! Technically, the equation below for lambda-r should have constants
! for the rain mass-diameter relation and mu, but these are currently
! the same as graupel, so we are cheating to avoid passing more
! constants from mp_thompson.
! !$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)
N0min = 1.E6
DO k=kme-1, kms, -1
if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
lamr = (1000.0*3.1415926536/6.0*cgg(3)/cgg(2)*temp_nr(i,k,j)/temp_qr(i,k,j))**(1./3.)
mvd_r = 3.672 / lamr ! Technically this should have (+mu_r)
mvd_r = MAX(100.E-6, MIN(mvd_r, 2.5E-3))
if (temp(i,k,j).lt.270.65 .and. mvd_r.gt.100.E-6) then
xslw1 = 4.01 + alog10(mvd_r)
else
xslw1 = 0.01
endif
ygra1 = 4.31 + alog10(max(5.E-5, temp_qg(i,k,j)))
zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))
N0exp = 10.**(zans1)
N0exp = MAX(DBLE(1.E4), MIN(N0exp, DBLE(1.E6)))
N0min = MIN(N0exp, N0min)
N0exp = N0min
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))
if (N0exp .ge. 1.E4 .AND. N0exp.le.1.E5) then
WRITE(outstring,*) 'GT-Debug, ', N0exp, temp_qr(i,k,j)*1000., temp_ng(i,k,j)
CALL wrf_debug (850, TRIM(outstring))
endif
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)
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
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
enddo
if (ng .ge. ngbins) then
hail_max = xxDg(ngbins)
elseif (xxDg(ng+1) .gt. 1.E-3) then
hail_max = 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
! COMPUTE THE NUMBER OF MASS GRID POINTS
no_points = float((ide-ids)*(jde-jds))
! SET START AND END POINTS FOR TILES
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
dmumax = 0.
DO ij = 1 , num_tiles
! print *, i_start(ij),i_end(ij),j_start(ij),j_end(ij)
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
dpsdt(i,j)=(p8w(i,kms,j)-pk1m(i,j))/dt
dmudt(i,j)=(mu_2(i,j)-mu_2m(i,j))/dt
if(abs(dmudt(i,j)*dt).gt.dmumax)then
dmumax=abs(dmudt(i,j)*dt)
idp=i
jdp=j
endif
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
! convert DMUMAX from (PA) to (bars) per time step
dmumax = dmumax*1.e-5
! compute global MAX
CALL wrf_dm_maxval ( dmumax, idp, jdp )
! print *, 'p8w(30,1,30),pk1m(30,30) : ', p8w(30,1,30),pk1m(30,30)
! print *, 'mu_2(30,30),mu_2m(30,30) : ', mu_2(30,30),mu_2m(30,30)
dpsdt_sum = 0.
dmudt_sum = 0.
DO j = jps, min(jpe,jde-1)
DO i = ips, min(ipe,ide-1)
dpsdt_sum = dpsdt_sum + abs(dpsdt(i,j))
dmudt_sum = dmudt_sum + abs(dmudt(i,j))
ENDDO
ENDDO
! compute global sum
dpsdt_sum = wrf_dm_sum_real ( dpsdt_sum )
dmudt_sum = wrf_dm_sum_real ( dmudt_sum )
! print *, 'dpsdt, dmudt : ', dpsdt_sum, dmudt_sum
IF ( diag_print .eq. 2 ) THEN
dardt_sum = 0.
drcdt_sum = 0.
drndt_sum = 0.
rainc_sum = 0.
raint_sum = 0.
rainnc_sum = 0.
sfcevp_sum = 0.
hfx_sum = 0.
lh_sum = 0.
raincmax = 0.
rainncmax = 0.
DO j = jps, min(jpe,jde-1)
DO i = ips, min(ipe,ide-1)
drcdt_sum = drcdt_sum + abs(raincv(i,j))
drndt_sum = drndt_sum + abs(rainncv(i,j))
dardt_sum = dardt_sum + abs(raincv(i,j)) + abs(rainncv(i,j))
rainc_sum = rainc_sum + abs(rainc(i,j))
! MAX for accumulated conv precip
IF(rainc(i,j).gt.raincmax)then
raincmax=rainc(i,j)
irc=i
jrc=j
ENDIF
rainnc_sum = rainnc_sum + abs(rainnc(i,j))
! MAX for accumulated resolved precip
IF(rainnc(i,j).gt.rainncmax)then
rainncmax=rainnc(i,j)
irnc=i
jrnc=j
ENDIF
raint_sum = raint_sum + abs(rainc(i,j)) + abs(rainnc(i,j))
sfcevp_sum = sfcevp_sum + abs(sfcevp(i,j))
hfx_sum = hfx_sum + abs(hfx(i,j))
lh_sum = lh_sum + abs(lh(i,j))
ENDDO
ENDDO
! compute global MAX
CALL wrf_dm_maxval ( raincmax, irc, jrc )
CALL wrf_dm_maxval ( rainncmax, irnc, jrnc )
! compute global sum
drcdt_sum = wrf_dm_sum_real ( drcdt_sum )
drndt_sum = wrf_dm_sum_real ( drndt_sum )
dardt_sum = wrf_dm_sum_real ( dardt_sum )
rainc_sum = wrf_dm_sum_real ( rainc_sum )
rainnc_sum = wrf_dm_sum_real ( rainnc_sum )
raint_sum = wrf_dm_sum_real ( raint_sum )
sfcevp_sum = wrf_dm_sum_real ( sfcevp_sum )
hfx_sum = wrf_dm_sum_real ( hfx_sum )
lh_sum = wrf_dm_sum_real ( lh_sum )
ENDIF
! print out the average values
CALL get_current_grid_name( grid_str )
#ifdef DM_PARALLEL
IF ( wrf_dm_on_monitor() ) THEN
#endif
WRITE(outstring,*) grid_str,'Domain average of dpsdt, dmudt (mb/3h): ', xtime, &
dpsdt_sum/no_points*108., &
dmudt_sum/no_points*108.
CALL wrf_message ( TRIM(outstring) )
WRITE(outstring,*) grid_str,'Max mu change time step: ', idp,jdp,dmumax
CALL wrf_message ( TRIM(outstring) )
IF ( diag_print .eq. 2) THEN
WRITE(outstring,*) grid_str,'Domain average of dardt, drcdt, drndt (mm/sec): ', xtime, &
dardt_sum/dt/no_points, &
drcdt_sum/dt/no_points, &
drndt_sum/dt/no_points
CALL wrf_message ( TRIM(outstring) )
WRITE(outstring,*) grid_str,'Domain average of rt_sum, rc_sum, rnc_sum (mm): ', xtime, &
raint_sum/no_points, &
rainc_sum/no_points, &
rainnc_sum/no_points
CALL wrf_message ( TRIM(outstring) )
WRITE(outstring,*) grid_str,'Max Accum Resolved Precip, I,J (mm): ' ,&
rainncmax,irnc,jrnc
CALL wrf_message ( TRIM(outstring) )
WRITE(outstring,*) grid_str,'Max Accum Convective Precip, I,J (mm): ' ,&
raincmax,irc,jrc
CALL wrf_message ( TRIM(outstring) )
WRITE(outstring,*) grid_str,'Domain average of sfcevp, hfx, lh: ', xtime, &
sfcevp_sum/no_points, &
hfx_sum/no_points, &
lh_sum/no_points
CALL wrf_message ( TRIM(outstring) )
ENDIF
#ifdef DM_PARALLEL
ENDIF
#endif
ENDIF
! save values at this time step
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij,i,j )
DO ij = 1 , num_tiles
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
pk1m(i,j)=p8w(i,kms,j)
mu_2m(i,j)=mu_2(i,j)
ENDDO
ENDDO
IF ( xtime .lt. 0.0001 ) THEN
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
dpsdt(i,j)=0.
dmudt(i,j)=0.
ENDDO
ENDDO
ENDIF
ENDDO
!$OMP END PARALLEL DO
END SUBROUTINE diagnostic_output_calc
!+---+-----------------------------------------------------------------+
REAL FUNCTION GAMMLN(XX) 4
! --- 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) 20
IMPLICIT NONE
REAL, INTENT(IN):: y
WGAMMA = EXP(GAMMLN(y))
END FUNCTION WGAMMA
!+---+-----------------------------------------------------------------+
END MODULE module_diag_misc
#endif