#if ( NMM_CORE == 1)
MODULE module_diag_misc
CONTAINS
SUBROUTINE diag_misc_stub
END SUBROUTINE diag_misc_stub
END MODULE module_diag_misc
#else
!WRF:MEDIATION_LAYER:PHYSICS
!
MODULE module_diag_misc
CONTAINS
SUBROUTINE diagnostic_output_calc( &
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 &
,prec_acc_c, prec_acc_nc, snow_acc_nc &
,snowncv, prec_acc_dt, curr_secs2 &
,nwp_diagnostics, diagflag &
,history_interval &
,itimestep &
,q2,psfc &
,t2max,t2min,rh2max,rh2min &
,t,u10,v10,w,p,pb &
,wspd10max &
,up_heli_max25 &
,up_heli_min25 &
,up_heli_max03 &
,up_heli_min03 &
,w_up_max,w_dn_max &
,znw,w_colmean &
,numcolpts,w_mean &
,grpl_max,grpl_colint,refd_max,refl_10cm &
,refdm10c_max,refdm10c_calc &
,qg_curr,qs_curr,qr_curr &
,rho,ph,phb,g &
)
!----------------------------------------------------------------------
USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval
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 ) :: 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 &
, Q2 &
, PSFC
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: DPSDT &
, DMUDT &
, RAINNC &
, RAINC &
, MU_2M &
, PK1M &
, T2MAX &
, T2MIN &
, RH2MAX &
, RH2MIN
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
INTEGER :: prfreq
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
REAL :: rh2,brack1,brack2,brack3,brack,vpres,satvpres
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 &
,qr_curr &
,qs_curr &
,rho &
,refl_10cm &
,t,ph,phb &
,p,pb
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_max25 &
,up_heli_min25 &
,up_heli_max03 &
,up_heli_min03 &
,w_up_max,w_dn_max &
,w_colmean,numcolpts,w_mean &
,grpl_max,grpl_colint &
,refd_max,refdm10c_max
INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
refdm10c_calc
INTEGER :: idump
REAL :: time_minutes
REAL :: afrac
REAL :: uh, vh, wind_vel
REAL :: depth
REAL :: tval2,tval1
REAL :: qr_val, qs_val, qg_val, rho_val
REAL :: thr1mx, thr2mx, thr1, thr2
REAL :: dbz,dbzr,dbzs,dbzg
REAL, PARAMETER :: tm10c=263.15
REAL, PARAMETER :: capa=0.28589641e0
! REAL, PARAMETER :: t0 = 273.155
REAL, PARAMETER :: t0 = 300.00 ! should be t0base value
REAL, PARAMETER :: dbzmin=-10.0 ! dcd
!-----------------------------------------------------------------
! 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_max25(i,j) = 0.
up_heli_min25(i,j) = 0.
up_heli_max03(i,j) = 0.
up_heli_min03(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.
refdm10c_max(i,j) = 0.
t2max(i,j) = -9999.
t2min(i,j) = 9999.
rh2max(i,j) = -9999.
rh2min(i,j) = 9999.
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 100 mb
IF ( p8w(i,k,j) .GT. 10000. .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. 10000. .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 max/min 2 m temperatures
IF ( t2(i,j) .GT. t2max(i,j) ) THEN
t2max(i,j) = t2(i,j)
ENDIF
IF ( t2(i,j) .LT. t2min(i,j) ) THEN
t2min(i,j) = t2(i,j)
ENDIF
! Calculate the max/min 2 m RH
vpres=(q2(I,J)/(q2(I,J)+0.622))*psfc(I,J)
brack1=2.5e6/461
brack2=1/273.15
brack3=1/t2(i,j)
brack=brack1*(brack2-brack3)
satvpres=611.*exp(brack)
rh2=vpres/satvpres
rh2=min(rh2,0.99)
rh2=max(rh2,0.01)
if (rh2 .lt. 0 .or. rh2 .gt. 1.05) then
write(0,*) 'i,j,rh2: ',i,j, rh2
call wrf_error_fatal("bad rh2 value")
endif
IF ( rh2 .GT. rh2max(i,j) ) THEN
rh2max(i,j) = rh2
ENDIF
IF ( rh2 .LT. rh2min(i,j) ) THEN
rh2min(i,j) = rh2
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,10,j) .GT. refd_max(i,j) ) THEN
refd_max(i,j) = refl_10cm(i,10,j)
ENDIF
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
! Calculate the derived radar reflectivity at the -10 deg C
! level between output times -- used in CI algorithm
! !$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,-1
DO i=i_start(ij),i_end(ij)
tval1 = ( t(i,k,j) + t0 ) * &
(( p(i,k,j) + pb(i,k,j) ) * 1.e-5 ) ** capa
tval2 = ( t(i,k-1,j) + t0 ) * &
(( p(i,k-1,j) + pb(i,k-1,j) ) * 1.e-5 ) ** capa
IF ( tval1 .LE. tm10c .AND. tval2 .GE. tm10c .AND. &
refdm10c_calc(i,j) .EQ. 0 ) THEN
! refdm10c_calc(i,j) = 1
! afrac = ( tm10c - tval2 ) / ( tval1 - tval2 )
! qr_val = qr_curr(i,k-1,j) + afrac * ( qr_curr(i,k,j) - qr_curr(i,k-1,j))
! qs_val = qs_curr(i,k-1,j) + afrac * ( qs_curr(i,k,j) - qs_curr(i,k-1,j))
! qg_val = qg_curr(i,k-1,j) + afrac * ( qg_curr(i,k,j) - qg_curr(i,k-1,j))
! rho_val = rho(i,k-1,j) + afrac * ( rho(i,k,j) - rho(i,k-1,j) )
! Derived radar reflectivity between output times
! dbzr = 0.
! dbzs = 0.
! dbzg = 0.
! guessing these values are Thompson specific - modify for HRW WSM6??
! IF ( qr_val .GT. 0. ) then
! dbzr = (( qr_val * rho_val ) ** 1.75 )* &
! 3.630803e-9 * 1.e18
! ENDIF
! IF ( qs_val .GT. 0. ) then
! dbzs = (( qs_val * rho_val ) ** 1.75 )* &
! 2.18500e-10 * 1.e18
! ENDIF
! IF ( qg_val .GT. 0. ) then
! dbzg = (( qg_val * rho_val ) ** 1.75 )* &
! 1.033267e-9 * 1.e18
! ENDIF
! dbz = dbzr + dbzs + dbzg
! IF ( dbz .GT. 0. ) dbz = 10.0 * LOG10 ( dbz )
! just use nearest level of reflectivity??
! taking max of two nearest levels to mimic what is done with NMMB, and to help
! make the max > instantaneous.
dbz = max(refl_10cm(i,k,j),refl_10cm(i,k-1,j))
dBz = MAX ( dbzmin, dbz )
IF ( dbz .GT. refdm10c_max(i,j) ) THEN
refdm10c_max(i,j) = dbz
ENDIF
ENDIF
! NSSL
ENDDO
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
ENDIF
if (diag_print .eq. 0 ) return
IF ( xtime .ne. 0. ) THEN
if(diag_print.eq.1) then
prfreq = dt
! prfreq = max(2,int(dt/60.)) ! in min
else
prfreq=10 ! in min
endif
IF (MOD(nint(dt),prfreq) == 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 ! print frequency
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
END MODULE module_diag_misc
#endif