#if ( NMM_CORE == 1) MODULE module_diag_gsd CONTAINS SUBROUTINE diag_gsd_stub END SUBROUTINE diag_gsd_stub END MODULE module_diag_gsd #else !WRF:MEDIATION_LAYER:PHYSICS ! MODULE module_diag_gsd USE module_model_constants PRIVATE :: WGAMMA PRIVATE :: GAMMLN CONTAINS SUBROUTINE gsd_refl_first_time_step( & ims,ime, jms,jme, kms,kme, kts,kte & ,i_start,i_end,j_start,j_end,num_tiles & ,mp_physics, th, pii, p, rho_curr, qv_curr & ,qc_curr, qnr_curr, qr_curr, qs_curr, qg_curr & ! ,refl_10cm, vt_dbz_wt & ,refl_10cm & ) !---------------------------------------------------------------------- USE module_state_description, ONLY : THOMPSON,THOMPSONAERO USE module_mp_thompson, ONLY : calc_refl10cm IMPLICIT NONE !====================================================================== ! dcd subroutine to compute reflectivity and fall speed for wrfout at t=0 ! currently works for Thompson microphysics scheme only ! ! Definitions !----------- ! Input: !-- 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 !-- kms start index for k in memory !-- kme end index for k in memory !-- kts !-- kte !-- 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 !-- num_tiles number of tiles !-- mp_physics index of precipitation microphysics scheme !-- th !-- pii !-- p !-- rho_curr air density (kg m-3) !-- qv_curr water vapor mixing ratio (kg kg-1) !-- qc_curr cloud water mixing ratio (kg kg-1) !-- qnr_curr rain number concentration (kg-1) !-- qr_curr rain water mixing ratio (kg kg-1) !-- qs_curr snow mixing ratio (kg kg-1) !-- qg_curr graupel mixing ratio (kg kg-1) ! ! Output: !-- refl_10cm radar reflectivity factor for 10-cm radar (dBZ) !-- vt_dbz_wt reflectivity-weighted precipitation fall speed (m s-1) !====================================================================== INTEGER, INTENT(IN ) :: & ims,ime, jms,jme, kms,kme, kts,kte, & num_tiles INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & & i_start,i_end,j_start,j_end INTEGER, INTENT(IN) :: mp_physics REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: & qv_curr,qr_curr,qs_curr,qg_curr & ,qc_curr,qnr_curr & ,rho_curr,th,pii,p REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT) :: refl_10cm ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT) :: vt_dbz_wt REAL :: rand1 INTEGER :: i,j,k,ij REAL, DIMENSION(kts:kte) :: qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d, dbz, fallspd LOGICAL :: qnr_initialized IF ( mp_physics .EQ. THOMPSON .or. mp_physics .EQ. THOMPSONAERO ) THEN qnr_initialized = .true. !$OMP PARALLEL DO & !$OMP PRIVATE ( ij,qnr_initialized,i,j,k,rand1,qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d, dbz, fallspd ) refl_loop: DO ij = 1 , num_tiles rand1 = 0.0 DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) DO k = kts, kte IF ( (qnr_curr(i,k,j) .le. 0.0) .and. (qr_curr(i,k,j) .ge. 1.e-12) ) THEN ! Rain number concentration was not initialized properly, but we would still ! like a reflectivity estimate. Let's arbitrarily assume the drops have a ! Marshall-Palmer distribution with intercept parameter 8000000 m**-4. ! The desired units of qnr are kg**-1. nr1d(k) = & (8.e6 / rho_curr(i,k,j)) * & (rho_curr(i,k,j) * qr_curr(i,k,j) / (3.14159 * 1000.0 * 8.e6) )**0.25 qnr_initialized = .false. ELSE nr1d(k) = qnr_curr(i,k,j) ENDIF t1d(k) = th(i,k,j)*pii(i,k,j) p1d(k) = p(i,k,j) qv1d(k) = qv_curr(i,k,j) qc1d(k) = qc_curr(i,k,j) qr1d(k) = qr_curr(i,k,j) qs1d(k) = qs_curr(i,k,j) qg1d(k) = qg_curr(i,k,j) ENDDO call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & t1d, p1d, dbz, fallspd, rand1, kts, kte, i, j, .true.) DO k = kts, kte refl_10cm(i,k,j) = MAX(-35., dbz(k)) !vt_dbz_wt(i,k,j) = fallspd(k) ENDDO ENDDO ENDDO ENDDO refl_loop IF (.not. qnr_initialized) THEN print *, 'QNR not initialized. Marshall-Palmer distribution will be assumed for reflectivity calculation.' ENDIF ELSE print*, 'mp_physics = ', mp_physics CALL wrf_error_fatal('This microphysics scheme is not supported in refl_first_time_step.') ENDIF ! IF ( mp_physics .EQ. THOMPSON ) END SUBROUTINE gsd_refl_first_time_step SUBROUTINE gsd_diagnostic_driver( & 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 & ,gsd_diagnostics & ,dpsdt,dmudt & ,p8w,pk1m,mu_2,mu_2m & ,u,v, temp & ,raincv,rainncv,rainc,rainnc,frain & ,i_rainc,i_rainnc & ! Optional - bucketr_opt =1 from bucket_mm ,hfx,sfcevp,lh & ,ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC & ! Optional - SW radiation package ,ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC & ! Optional - SW radiation package ,ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC & ! Optional - LW radiation package ,ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC & ! Optional - LW radiation package ,I_ACSWUPT,I_ACSWUPTC,I_ACSWDNT,I_ACSWDNTC & ! Optional - bucketf_opt =1 from bucket_J ,I_ACSWUPB,I_ACSWUPBC,I_ACSWDNB,I_ACSWDNBC & ! Optional - bucketf_opt =1 from bucket_J ,I_ACLWUPT,I_ACLWUPTC,I_ACLWDNT,I_ACLWDNTC & ! Optional - bucketf_opt =1 from bucket_J ,I_ACLWUPB,I_ACLWUPBC,I_ACLWDNB,I_ACLWDNBC & ! Optional - bucketf_opt =1 from bucket_J ,dt,xtime,sbw,t2 & ,diag_print & ! Option to print out time-series ,dfi_stage & ,bucket_mm, bucket_J & ,mphysics_opt & ,gsfcgce_hail, gsfcgce_2ice & ,mpuse_hail & ,nssl_cnoh, nssl_cnohl & !namelist parameters ,nssl_rho_qh, nssl_rho_qhl & !namelist parameters ,nssl_alphah, nssl_alphahl & !namelist parameters ,prec_acc_c, prec_acc_nc, snow_acc_nc & ! Optional - prec_acc_opt =1 ,prec_acc_c1, prec_acc_nc1, snow_acc_nc1 & ! Optional - prec_acc_opt =1 ,graup_acc_nc, graup_acc_nc1, graupelncv & ,acsnow, acgraup, acrunoff,acfrain & ,sfcrunoff, udrunoff & ,snowncv, prec_acc_dt, prec_acc_dt1, curr_secs2 & ,diagflag & ,history_interval & ,reset_interval1 & ,itimestep & ,ntimesteps & ,t,u10,v10,w,u_phy,v_phy,p,pb & !JOE ,wspd10max,wspd10umax,wspd10vmax & ! Optional - gsd_diagnostics =1 ,wspd80max,wspd80umax,wspd80vmax & ,wspd10,wspd80,LWP,IWP,maxcldfra & !JOE: Optional - gsd_diagnostics =1 ,up_heli_max,up_heli_min & ! Optional - gsd_diagnostics =1 ,up_heli_max16,up_heli_min16 & ,up_heli_max02,up_heli_min02 & ,up_heli_max03,up_heli_min03 & ,rel_vort_max,rel_vort_max01 & ,w_up_max,w_dn_max & ! Optional - gsd_diagnostics =1 ,znw,w_colmean & ,numcolpts,w_mean_sum,w_mean & ,grpl_max,grpl_colint,refd_max,refl_10cm & ,refdm10c_max,refdm10c_calc & ,ltg1_max,ltg2_max,ltg3_max & ,ltg1,ltg2,ltg3 & ,ltg1_calc & ,totice_colint,msft & ,refl_10cm_1km,refl_10cm_4km,composite_refl_10cm & ! dcd ,ht & ! dcd ,hail_maxk1,hail_max2d & ! Optional - gsd_diagnostics =1 ,qg_curr,qs_curr & ,ng_curr,qh_curr,nh_curr,qr_curr,nr_curr & ! Optional (gthompsn) ,qc_curr,qi_curr & !JOE ,rho,ph,phb,g & ,cldfra,icloud_bl,cldfra_bl,qc_bl & !JOE ) !---------------------------------------------------------------------- USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval USE module_wrf_error 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 !-- GRAUP_ACC_NC accumulated explicit graupel precip over accumulation time prec_acc_dt !-- PREC_ACC_C1 accumulated convective precip over accumulation time prec_acc_dt1 !-- PREC_ACC_NC1 accumulated explicit precip over accumulation time prec_acc_dt1 !-- SNOW_ACC_NC1 accumulated explicit snow precip over accumulation time prec_acc_dt1 !-- GRAUP_ACC_NC1 accumulated explicit graupel precip over accumulation time prec_acc_dt1 !-- PREC_ACC_DT precip accumulation time, default is 60 min !-- PREC_ACC_DT1 precip accumulation time, default is 15 min !-- ACSNOW run-total accumulated snow !-- ACGRAUP run-total accumulated graupel !-- ACFRAIN run-total accumulated freezing rain (T<=273) !-- CURR_SECS2 Time (s) since the beginning of the restart !-- NWP_DIAGNOSTICS = 1, compute hourly maximum fields *** JOE: removed - compute always !-- 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 !-- REFL_10CM_1KM reflectivity at 1 km AGL (dBZ) !-- REFL_10CM_4KM reflectivity at 4 km AGL (dBZ) !-- COMPOSITE_REFL_10CM maximum reflectivity in column (dBZ) !-- HT ht terrain height (m) ! !-- 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 INTEGER, INTENT(IN ) :: dfi_stage 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 ) :: icloud_bl INTEGER, INTENT(IN ) :: gsd_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 & , FRAIN & , SNOWNCV & ,GRAUPELNCV & , HFX & , LH & , SFCEVP & , T2 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, & INTENT(IN ) :: CLDFRA, CLDFRA_BL, QC_BL 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 ! these optional radiation variables should only be used when ! gsd_diagnostics == 2 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 & ,PREC_ACC_C1, PREC_ACC_NC1, SNOW_ACC_NC1 & , GRAUP_ACC_NC & , GRAUP_ACC_NC1& , ACSNOW & , ACFRAIN & , ACGRAUP & , ACRUNOFF & , SFCRUNOFF & , UDRUNOFF REAL, OPTIONAL, INTENT(IN):: PREC_ACC_DT, PREC_ACC_DT1, CURR_SECS2 INTEGER :: i,j,k,its,ite,jts,jte,ij INTEGER :: idp,jdp,irc,jrc,irnc,jrnc,isnh,jsnh INTEGER :: prfreq REAL, PARAMETER :: dbzmin=-10.0 ! dcd INTEGER :: k1, k2 ! dcd REAL :: h_agl ! dcd REAL :: cuprate ! dcd REAL :: ze, ze_conv, dbz_sum ! dcd LOGICAL :: found_1km, found_4km ! dcd 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,reset_interval1,itimestep REAL, DIMENSION( kms:kme ), INTENT(IN) :: & znw REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: & w & ,temp & ,qs_curr,qg_curr & ,qc_curr,qi_curr & ,rho & ,ph,phb & ,p,pb & ,t,u_phy,v_phy REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: & !dcd refl_10cm 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 & ,msft REAL, INTENT(IN) :: g REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: & wspd10max,wspd10umax,wspd10vmax & ,wspd80max,wspd80umax,wspd80vmax & ,up_heli_max,up_heli_min & ,up_heli_max16,up_heli_min16 & ,up_heli_max02,up_heli_min02 & ,up_heli_max03,up_heli_min03 & ,rel_vort_max,rel_vort_max01 & ,w_up_max,w_dn_max & ,w_colmean,numcolpts,w_mean_sum,w_mean & ,grpl_max,grpl_colint & ,hail_maxk1,hail_max2d & ,refd_max & ,refdm10c_max & ,ltg1_max,ltg2_max,ltg3_max & ,ltg1,ltg2,ltg3 & ,totice_colint & ,refl_10cm_1km,refl_10cm_4km,composite_refl_10cm & ! dcd ,ht & ! dcd ,wspd10,wspd80,LWP,IWP & ,maxcldfra INTEGER, INTENT(INOUT) :: ntimesteps REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: temp_qg, temp_ng, temp_qr, temp_nr INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: & ltg1_calc,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 INTEGER :: ithr1mx,jthr1mx,ithr2mx,jthr2mx INTEGER :: ithr1,jthr1,ithr2,jthr2 ! Per-thread extrema REAL :: t_thr1mx(num_tiles), t_thr2mx(num_tiles), t_thr1(num_tiles), t_thr2(num_tiles) INTEGER :: t_ithr1mx(num_tiles),t_jthr1mx(num_tiles),t_ithr2mx(num_tiles),t_jthr2mx(num_tiles) INTEGER :: t_ithr1(num_tiles),t_jthr1(num_tiles),t_ithr2(num_tiles),t_jthr2(num_tiles) ! Derived composite radar reflectivity every 5 minutes REAL, PARAMETER :: refddt=5. ! Constants used in ltg threat calculations REAL, PARAMETER :: coef1=0.042*1000.*1.22 REAL, PARAMETER :: coef2=0.20*1.22 REAL, PARAMETER :: clim1=1.50 REAL, PARAMETER :: clim2=0.40*1.22 REAL, PARAMETER :: clim3=0.02*1.22 REAL, PARAMETER :: capa=0.28589641e0 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 :: dp,sum1,sum2,sum3,zm,zm0,wt1,ter,u100,v100,t100, & qctotal,qitotal,dumy1,dumy2 !JOE #ifdef MINIMAL_DEBUGGING LOGICAL, PARAMETER :: print10flag=.false. LOGICAL, PARAMETER :: print100flag=.false. LOGICAL, PARAMETER :: print150flag=.false. LOGICAL, PARAMETER :: print350flag=.false. #else LOGICAL :: print10flag,print100flag,print150flag,print350flag print10flag=wrf_at_debug_level(10) print100flag=wrf_at_debug_level(100) print150flag=wrf_at_debug_level(150) print350flag=wrf_at_debug_level(350) #endif !----------------------------------------------------------------- ! 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 PRIVATE ( ij, i, j ) 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 !LIMIT THESE RADIATION DIAGS TO NON-REALTIME CONFIGURATION IF (gsd_diagnostics .eq. 1) THEN !only set them == 0 for realtime application !and avoid using them again 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) = 1 i_acswuptc(i,j) = 1 i_acswdnt(i,j) = 1 i_acswdntc(i,j) = 1 i_acswupb(i,j) = 1 i_acswupbc(i,j) = 1 i_acswdnb(i,j) = 1 i_acswdnbc(i,j) = 1 acswupt(i,j) = 0.0 acswuptc(i,j) = 0.0 acswdnt(i,j) = 0.0 acswdntc(i,j) = 0.0 acswupb(i,j) = 0.0 acswupbc(i,j) = 0.0 acswdnb(i,j) = 0.0 acswdnbc(i,j) = 0.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) = 1 i_aclwuptc(i,j) = 1 i_aclwdnt(i,j) = 1 i_aclwdntc(i,j) = 1 i_aclwupb(i,j) = 1 i_aclwupbc(i,j) = 1 i_aclwdnb(i,j) = 1 i_aclwdnbc(i,j) = 1 aclwupt(i,j) = 0.0 aclwuptc(i,j) = 0.0 aclwdnt(i,j) = 0.0 aclwdntc(i,j) = 0.0 aclwupb(i,j) = 0.0 aclwupbc(i,j) = 0.0 aclwdnb(i,j) = 0.0 aclwdnbc(i,j) = 0.0 ENDDO ENDDO ENDIF ELSEIF (gsd_diagnostics .eq. 2) THEN !For research (gsd_diagnostics=2), calculate these fields 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 ENDIF !gsd_diagnostic check ENDDO ENDIF ! Compute precipitation accumulation in a given time window: prec_acc_dt1 IF (prec_acc_dt1 .gt. 0.) THEN !$OMP PARALLEL DO 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) IF (curr_secs2 == 0 .or. mod(curr_secs2, 60.* prec_acc_dt1) == 0.) THEN !print *,'zero out buckets for dt=15min', curr_secs2 prec_acc_c1(i,j) = 0. prec_acc_nc1(i,j) = 0. snow_acc_nc1(i,j) = 0. graup_acc_nc1(i,j) = 0. ENDIF prec_acc_c1(i,j) = prec_acc_c1(i,j) + RAINCV(i,j) prec_acc_nc1(i,j) = prec_acc_nc1(i,j) + RAINNCV(i,j) prec_acc_c1(i,j) = MAX (prec_acc_c1(i,j), 0.0) prec_acc_nc1(i,j) = MAX (prec_acc_nc1(i,j), 0.0) snow_acc_nc1(i,j) = snow_acc_nc1(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_nc1(i,j) = snow_acc_nc1(i,j) + RAINCV(i,j) ENDIF graup_acc_nc1(i,j) = graup_acc_nc1(i,j) + graupelncv(i,j) snow_acc_nc1(i,j) = MAX (snow_acc_nc1(i,j), 0.0) graup_acc_nc1(i,j) = MAX (graup_acc_nc1(i,j), 0.0) ENDDO ENDDO ENDDO ENDIF ! IF (prec_acc_dt .gt. 0.) THEN !$OMP PARALLEL DO 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) ! IF (mod(curr_secs2, 60.* prec_acc_dt) == 0.) THEN IF (curr_secs2 == 0 .or. mod(curr_secs2, 60.* prec_acc_dt) == 0.) THEN ! print *,'zero out buckets', curr_secs2 prec_acc_c(i,j) = 0. prec_acc_nc(i,j) = 0. snow_acc_nc(i,j) = 0. graup_acc_nc(i,j) = 0. sfcrunoff (i,j) = 0. udrunoff (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) acsnow(i,j) = acsnow(i,j) + SNOWNCV(I,J) graup_acc_nc(i,j) = graup_acc_nc(i,j) + graupelncv(i,j) acgraup(i,j) = acgraup(i,j) + graupelncv(i,j) acfrain(i,j) = acfrain(i,j) + frain(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) acsnow(i,j) = acsnow(i,j) + RAINCV(i,j) ENDIF snow_acc_nc(i,j) = MAX (snow_acc_nc(i,j), 0.0) graup_acc_nc(i,j) = MAX (graup_acc_nc(i,j), 0.0) acsnow(i,j) = MAX (acsnow(i,j), 0.0) acgraup(i,j)= MAX (acgraup(i,j), 0.0) ENDDO ENDDO ENDDO ENDIF ! NSSL !JOE: removed IF ( nwp_diagnostics .EQ. 1 ) THEN idump = (history_interval * 60.) / dt time_minutes = MOD(xtime,60.) ! 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 IF ( MOD(curr_secs2,60.*reset_interval1) .eq. 0. ) THEN if(print10flag) then WRITE(outstring,*) 'Resetting max arrays for domain with dt = ', dt WRITE(outstring,*) 'Diagnostics: curr_sec2 = ', curr_secs2 CALL wrf_debug ( 10,TRIM(outstring) ) endif ntimesteps = 0 !$OMP PARALLEL DO 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) wspd10max(i,j) = 0. wspd10umax(i,j) = 0. wspd10vmax(i,j) = 0. wspd80max(i,j) = 0. wspd80umax(i,j) = 0. wspd80vmax(i,j) = 0. up_heli_max(i,j) = 0. up_heli_min(i,j) = 0. up_heli_max16(i,j) = 0. up_heli_min16(i,j) = 0. up_heli_max02(i,j) = 0. up_heli_min02(i,j) = 0. up_heli_max03(i,j) = 0. up_heli_min03(i,j) = 0. rel_vort_max(i,j) = 0. rel_vort_max01(i,j) = 0. w_up_max(i,j) = 0. w_dn_max(i,j) = 0. w_mean_sum(i,j) = 0. w_mean(i,j) = 0. grpl_max(i,j) = 0. refd_max(i,j) = 0. refdm10c_max(i,j) = 0. ltg1_max(i,j) = 0. ltg2_max(i,j) = 0. ltg3_max(i,j) = 0. hail_maxk1(i,j) = 0. hail_max2d(i,j) = 0. WSPD10(i,j)=-99. WSPD80(i,j)=-99. LWP(i,j)=0. IWP(i,j)=0. MAXCLDFRA(i,j)=0. ENDDO ENDDO ENDDO ENDIF ntimesteps = ntimesteps + 1 !$OMP PARALLEL DO 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) ! Zero some accounting arrays that will be used below w_colmean(i,j) = 0. numcolpts(i,j) = 0. grpl_colint(i,j) = 0. totice_colint(i,j) = 0. refdm10c_calc(i,j) = 0 ltg1_calc(i,j) = 0 ENDDO ENDDO ENDDO !$OMP PARALLEL DO PRIVATE ( ij, i, j, k ) 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 PARALLEL DO PRIVATE ( ij, i, j, k, depth, tval1, tval2 ) 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) totice_colint(i,j) = totice_colint(i,j) + ( qi_curr(i,k,j) + \ qs_curr(i,k,j) + qg_curr(i,k,j) ) * \ depth * rho(i,k,j) IF ( ltg1_calc(i,j) .EQ. 0 ) THEN 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+tval2)*0.5) .LT. 258.15 ) THEN ltg1_calc(i,j) = 1 ltg1(i,j) = coef1*w(i,k,j)* \ ((qg_curr(i,k,j)+qg_curr(i,k+1,j))*0.5)/msft(i,j) IF ( ltg1(i,j) .LT. clim1 ) ltg1(i,j) = 0. IF ( ltg1(i,j) .GT. ltg1_max(i,j) ) THEN ltg1_max(i,j) = ltg1(i,j) ENDIF ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO !$OMP PARALLEL DO PRIVATE ( ij, i, j, k, wind_vel, dbzr, dbzs, dbzg, dbz ) 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 wspd10umax(i,j) = u10(i,j) wspd10vmax(i,j) = v10(i,j) ENDIF ! Calculate the column mean vertical velocity between output times w_mean_sum(i,j) = w_mean_sum(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 IF ( ntimesteps .gt. 0 ) THEN w_mean(i,j) = w_mean_sum(i,j) / ntimesteps 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 lightning threat 2 ltg2(i,j) = coef2 * totice_colint(i,j) / msft(i,j) IF ( ltg2(i,j) .LT. clim2 ) ltg2(i,j) = 0. IF ( ltg2(i,j) .GT. ltg2_max(i,j) ) THEN ltg2_max(i,j) = ltg2(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 ! Use k=kms for now, but can be reset to whatever level seems best dbzr = 0. dbzs = 0. dbzg = 0. IF ( qr_curr(i,kms,j) .GT. 0. ) then dbzr = (( qr_curr(i,kms,j) * rho(i,kms,j) ) ** 1.75 )* & 3.630803e-9 * 1.e18 ENDIF IF ( qs_curr(i,kms,j) .GT. 0. ) then dbzs = (( qs_curr(i,kms,j) * rho(i,kms,j) ) ** 1.75 )* & 2.18500e-10 * 1.e18 ENDIF IF ( qg_curr(i,kms,j) .GT. 0. ) then dbzg = (( qg_curr(i,kms,j) * rho(i,kms,j) ) ** 1.75 )* & 1.033267e-9 * 1.e18 ENDIF dbz = dbzr + dbzs + dbzg IF ( dbz .GT. 0. ) dbz = 10.0 * LOG10 ( dbz ) dbz = MAX ( dbzmin, dbz ) IF ( dbz .GT. refd_max(i,j) ) THEN refd_max(i,j) = dbz ENDIF ENDDO ENDDO ENDDO ! Calculate the derived radar reflectivity at the -10 deg C ! level between output times -- used in CI algorithm !$OMP PARALLEL DO PRIVATE ( ij, i, j, k, tval1, tval2, afrac, qr_val, qs_val, qg_val, rho_val, dbzr, dbzs, dbzg, dbz ) 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 ( j.eq.j_start(ij) .and. i.eq.i_start(ij) )then ! print *,' k, tval1, tval2 = ', k,tval1,tval2,refdm10c_calc(i,j) ! endif 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 ! Use k=kms for now, but can be reset to whatever level seems best dbzr = 0. dbzs = 0. dbzg = 0. 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 ) dBz = MAX ( dbzmin, dbz ) IF ( dbz .GT. refdm10c_max(i,j) ) THEN refdm10c_max(i,j) = dbz ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO ! Calculate the max lightning threat 3 -- max values for ! lightning threats 1 and 2 must be determined first ! Make the calculation for the instantaneous lightning ! threat fields as well call tiles_max_loc_ij(ltg1_max, thr1mx, ithr1mx, jthr1mx ) CALL tiles_max_loc_ij(ltg2_max, thr2mx, ithr2mx, jthr2mx ) CALL tiles_max_loc_ij(ltg1, thr1, ithr1, jthr1 ) CALL tiles_max_loc_ij(ltg2, thr2, ithr2, jthr2 ) ! Now calculate the max and instantaneous lightning threat 3 !$OMP PARALLEL DO 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) IF ( thr2mx .GT. 0. ) THEN ltg3_max(i,j) = 0.95 * ltg1_max(i,j) + & 0.05 * ( thr1mx / thr2mx ) * ltg2_max(i,j) ELSE ltg3_max(i,j) = 0.95 * ltg1_max(i,j) ENDIF IF ( ltg3_max(i,j) .LT. clim3 ) ltg3_max(i,j) = 0. IF ( thr2 .GT. 0. ) THEN ltg3(i,j) = 0.95 * ltg1(i,j) + & 0.05 * ( thr1 / thr2 ) * ltg2(i,j) ELSE ltg3(i,j) = 0.95 * ltg1(i,j) ENDIF IF ( ltg3(i,j) .LT. clim3 ) ltg3(i,j) = 0. ENDDO ENDDO ENDDO ! dcd reflectivity diagnostics for RAP and HRRR IF (diagflag) THEN k1 = kms k2 = kme-1 !$OMP PARALLEL DO PRIVATE ( ij, i, j, ze_conv, cuprate, k, h_agl, ze, dbz_sum, found_1km, found_4km ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) composite_refl_10cm(i,j) = dbzmin refl_10cm_1km(i,j) = dbzmin refl_10cm_4km(i,j) = dbzmin found_1km = .false. found_4km = .false. ! reflectivity parameterization for parameterized convection (reference: Unipost MDLFLD.f) ze_conv = 0.0 cuprate = raincv(i,j) * 3600.0 / dt ! cu precip rate (mm/h) ze_conv = 300.0 * cuprate**1.4 DO k=k1,k2 ! approximation: use base-state geopotential to determine height h_agl = phb(i,k,j) / g - ht(i,j) ze = 10.0 ** (0.1 * refl_10cm(i,k,j)) dbz_sum = max(dbzmin, 10.0 * log10(ze + ze_conv)) refl_10cm(i,k,j) = dbz_sum IF ( (.not. found_1km) .and. (h_agl.gt.1000.0) ) THEN refl_10cm_1km(i,j) = dbz_sum found_1km = .true. ENDIF IF ( (.not. found_4km) .and. (h_agl.gt.4000.0) ) THEN refl_10cm_4km(i,j) = dbz_sum found_4km = .true. ENDIF IF ( dbz_sum .GT. composite_refl_10cm(i,j) ) THEN composite_refl_10cm(i,j) = dbz_sum ENDIF ENDDO ENDDO ENDDO ENDDO ENDIF !Diagnostics for PBL/Cloud Development IF (diagflag) THEN !$OMP PARALLEL DO PRIVATE ( ij, i, j, ter, zm, sum1, sum2, sum3, k, zm0, depth, wt1, dp, dumy1, dumy2, qctotal, qitotal ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) ter = ( ph(i,1,j) + phb(i,1,j) ) / g zm=0. sum1=0. sum2=0. sum3=0. DO k=kms,kme-1 !Note: ph & phb are defined at the layer interface, so we average to ! get a height (AMSL) at the middle of the model layers zm0=zm zm = 0.5*((ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j)) / g ) depth = zm-zm0 !10 m WIND SPEED IF (zm-ter > 10. .AND. WSPD10(i,j)==-99.) THEN IF (k==kms) THEN WSPD10(i,j)=SQRT(u_phy(i,k,j)**2 + v_phy(i,k,j)**2) ELSE wt1=(zm-ter - 10.)/depth ! portion of layer above 10 m wt1=MIN(MAX(0.,wt1),1.0) WSPD10(i,j)=(1.-wt1)*SQRT(u_phy(i,k,j)**2 + v_phy(i,k,j)**2) + & wt1*SQRT(u_phy(i,k-1,j)**2 + v_phy(i,k-1,j)**2) !print*,"z-ter=",zm-ter," wt1=",wt1 !print*,i,j," WSPD10=",WSPD10(i,j) ENDIF ENDIF !80 m WIND SPEED IF (zm-ter > 80. .AND. WSPD80(i,j)==-99.) THEN IF (k==kms) THEN WSPD80(i,j)=SQRT(u_phy(i,k,j)**2 + v_phy(i,k,j)**2) ELSE wt1=(zm-ter - 80.)/depth ! portion of layer above 80 m wt1=MIN(MAX(0.,wt1),1.0) WSPD80(i,j)=(1.-wt1)*SQRT(u_phy(i,k,j)**2 + v_phy(i,k,j)**2) + & wt1*SQRT(u_phy(i,k-1,j)**2 + v_phy(i,k-1,j)**2) !print*,"z-ter=",zm-ter," wt1=",wt1 !print*,i,j," WSPD80=",WSPD80(i,j) ENDIF ENDIF !JOE -add LWP & IWP & MAXCLDFRA; use temp (abs temperature in kelvin) dp=p8w(i,k,j)-p8w(i,k+1,j) IF (icloud_bl>0 .AND. PRESENT(QC_BL) .AND. & PRESENT(CLDFRA_BL) .AND. qc_curr(i,k,j)<1E-6 .AND. & qi_curr(i,k,j)<1E-8 .AND. CLDFRA_BL(i,k,j)>0.01) THEN !guard against fpes by truncating small values dumy1=QC_BL(i,k,j) IF (dumy1 < 1E-9)dumy1=0. dumy2=CLDFRA_BL(i,k,j) IF (dumy2 < 1E-2)dumy2=0. !use temperature-dependent sorting of BL clouds qctotal = dumy1*(MIN(1., MAX(0., (temp(i,k,j)-254.)/15.)))*dumy2 qitotal = dumy1*(1. - MIN(1., MAX(0., (temp(i,k,j)-254.)/15.)))*dumy2 ELSE qctotal = QC_CURR(i,k,j) qitotal = QI_CURR(i,k,j) ENDIF sum1 = sum1 + MAX((dp/g)*qctotal,0.0) sum2 = sum2 + MAX((dp/g)*qitotal,0.0) sum3 = MAX(sum3,CLDFRA_BL(i,k,j)) !isolate BL clouds only ENDDO LWP(i,j)=sum1*1000. ! kg m-2 --> g m-2 IWP(i,j)=sum2*1000. MAXCLDFRA(i,j)=sum3 ENDDO ENDDO ENDDO ENDIF !check for diag flag !End extra diags for cloud/PBL development !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ !..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. !+---+-----------------------------------------------------------------+ if(print100flag) then WRITE(outstring,*) 'GT-Diagnostics, computing max-hail diameter' CALL wrf_debug (100, TRIM(outstring)) endif IF (PRESENT(qh_curr)) THEN if(print150flag) then WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail mixing ratio' CALL wrf_debug (150, TRIM(outstring)) endif !$OMP PARALLEL DO PRIVATE ( ij, i, j, k ) 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 ELSE !$OMP PARALLEL DO PRIVATE ( ij, i, j, k ) 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 ENDIF IF (PRESENT(nh_curr)) THEN if(print150flag) then WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail number concentration' CALL wrf_debug (150, TRIM(outstring)) endif !$OMP PARALLEL DO PRIVATE ( ij, i, j, k ) 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 ELSEIF (PRESENT(ng_curr)) THEN if(print150flag) then WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel number concentration' endif !$OMP PARALLEL DO PRIVATE ( ij, i, j, k ) 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 ELSE !$OMP PARALLEL DO PRIVATE ( ij, i, j, k ) 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 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 PRIVATE ( ij, i, j, k, lam_exp ) 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 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 PRIVATE ( ij, i, j, k, lam_exp ) 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 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 PRIVATE ( ij, i, j, k, lam_exp ) 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 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 PRIVATE ( ij, i, j, k, lam_exp ) 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 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 PRIVATE ( ij, i, j, k, zans1, N0exp, lam_exp, lamg, N0_g ) 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 ! 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) if(print150flag) then WRITE(outstring,*) 'GT-Debug, using Milbrandt2mom, which has 2-moment hail' CALL wrf_debug (150, TRIM(outstring)) endif 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 PRIVATE ( ij, i, j, k, lam_exp ) 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 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 PRIVATE ( ij, i, j, k, lamg, N0_g, sum_ng, sum_t, f_d, hail_max, hail_max_sp, ng ) 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 #ifdef MODULE_DIAG_GSD_EXTREME_DEBUGGING if (hail_max .gt. 1E-2 .and. print350flag ) then !$OMP CRITICAL WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000. CALL wrf_debug (350, TRIM(outstring)) !$OMP END CRITICAL endif #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 endif !JOE:removed nwp_diagnostics ENDIF ! NSSL 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 dmumax = 0. !$OMP PARALLEL DO PRIVATE ( ij ) 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 !$OMP CRITICAL if(abs(dmudt(i,j)*dt).gt.dmumax)then dmumax=abs(dmudt(i,j)*dt) idp=i jdp=j endif !$OMP END CRITICAL endif ENDDO ENDDO ENDDO ! 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. !$OMP PARALLEL DO PRIVATE(ij, i, j) REDUCTION(+:drcdt_sum,drndt_sum,dardt_sum,rainc_sum,raint_sum,sfcevp_sum,hfx_sum,lh_sum) do ij=1,num_tiles DO j = j_start(ij), min(j_end(ij),jde-1) DO i = i_start(ij), min(i_end(ij),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 !$OMP CRITICAL IF(rainc(i,j).gt.raincmax)then raincmax=rainc(i,j) irc=i jrc=j ENDIF !$OMP END CRITICAL ENDIF rainnc_sum = rainnc_sum + abs(rainnc(i,j)) ! MAX for accumulated resolved precip IF(rainnc(i,j).gt.rainncmax)then !$OMP CRITICAL IF(rainnc(i,j).gt.rainncmax)then rainncmax=rainnc(i,j) irnc=i jrnc=j ENDIF !$OMP END CRITICAL 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 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 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 CONTAINS SUBROUTINE tiles_max_loc_ij(var,vmax,imax,jmax) real, intent(in) :: var(ims:ime,jms:jme) integer, intent(inout) :: imax,jmax real, intent(inout) :: vmax integer :: i, j, ij real :: t_vmax(num_tiles) integer :: t_imax(num_tiles), t_jmax(num_tiles) !$OMP PARALLEL DO PRIVATE(ij) do ij=1,num_tiles t_vmax(ij)=var(i_start(ij),j_start(ij)) t_imax(ij)=i_start(ij) t_jmax(ij)=j_start(ij) do j=j_start(ij),j_end(ij) do i=i_start(ij),i_end(ij) if(var(i,j)>t_vmax(ij)) then t_vmax(ij)=var(i,j) t_imax=i t_jmax=j endif enddo enddo enddo vmax=t_vmax(1) imax=t_imax(1) jmax=t_jmax(1) do ij=2,num_tiles if(t_vmax(ij)>vmax) then vmax=t_vmax(ij) imax=t_imax(ij) jmax=t_jmax(ij) endif enddo CALL wrf_dm_maxval(vmax, imax, jmax) END SUBROUTINE tiles_max_loc_ij END SUBROUTINE gsd_diagnostic_driver !+---+-----------------------------------------------------------------+ 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_gsd #endif