diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 4d133c9bfa..b8bb880ef8 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -1594,6 +1594,13 @@ state real GRAUPELNCV ij misc 1 - r "G state real HAILNCV ij misc 1 - r "HAILNCV" "TIME-STEP NONCONVECTIVE HAIL" "mm" state real refl_10cm ikj dyn_em 1 - hdu "refl_10cm" "Radar reflectivity (lamda = 10 cm)" "dBZ" state real mskf_refl_10cm ikj dyn_em 1 - hdu "mskf_refl_10cm" "Full Radar reflectivity (lamda = 10 cm)" "dBZ" +! jdduda +state real MAX_PRATE ij misc 1 - rh "MAX_PRATE" "Time-maximum instantaneous precipitation rate" "mm s-1" +state real MAX_PRATE_1MIN ij misc 1 - rh "MAX_PRATE_1MIN" "Time-maximum 1-minute averaged precipitation rate" "mm s-1" +state real MAX_PRATE_5MIN ij misc 1 - rh "MAX_PRATE_5MIN" "Time-maximum 5-minute averaged precipitation rate" "mm s-1" +state real MAX_PRATE_10MIN ij misc 1 - rh "MAX_PRATE_10MIN" "Time-maximum 10-minute averaged precipitation rate" "mm s-1" +state real MAX_COMPREF ij misc 1 - rh "MAX_COMPREF" "Time-maximum composite reflectivity" "dBZ" +! jdduda state real th_old ikj misc 1 - rusd "TH_OLD" "Old Value of Th" "K" state real qv_old ikj misc 1 - rusd "QV_OLD" "Old Value of qv" "kg kg-1" state real vmi3d ikj misc 1 - hdu "v_ice" "Mass-weighted ice fallspeed cat 1" "m s-1" @@ -3006,6 +3013,71 @@ rconfig real windfarm_tke_factor namelist,physics 1 0 # # Ideal case selection rconfig integer ideal_case namelist,ideal 1 0 rh "ideal_case" "" "" +# options that impact quarter_ss only +rconfig logical linear_bubble namelist,ideal max_domains .false. - "linear_bubble" "" "" +rconfig real line_fract namelist,ideal max_domains 1.0 - "line_fract" "" "" +rconfig integer line_seed namelist,ideal max_domains -1 - "line_seed" "" "" +rconfig logical init_vv namelist,ideal max_domains .false. - "init_vv" "superimpose an updraft core at the perturbation location" "" +rconfig real init_w_core namelist,ideal max_domains 0 - "init_w_core" "maximum vertical velocity at perturbation center" "" +rconfig real vv_lx namelist,ideal max_domains 10000. - "vv_lx" "horizontal length scale of vertical velocity perturbation" "" +rconfig real vv_ly namelist,ideal max_domains 10000. - "vv_ly" "horizontal length scale of vertical velocity perturbation" "" +rconfig real vv_lz namelist,ideal max_domains 2000. - "vv_lz" "vertical length scale of vertical velocity perturbation" "" +rconfig logical init_div namelist,ideal max_domains .false. - "init_max_div" "superimpose a divergent region over the perturbation" "" +rconfig real init_max_div namelist,ideal max_domains 0. - "init_max_div" "maximum divergence at perturbation center" "" +rconfig real div_lx namelist,ideal max_domains 10000. - "div_lx" "horizontal scale of divergence perturbation" "" +rconfig real div_ly namelist,ideal max_domains 10000. - "div_ly" "horizontal length scale of divergence perturbation" "" +rconfig real div_depth namelist,ideal max_domains 2000. - "div_depth" "depth (from surface) of forced divergence region" "" +rconfig integer num_perts namelist,ideal max_domains 0 - "num_perts" "" "" +rconfig integer pert_1_x namelist,ideal max_domains 0 - "pert_1_x" "" "" +rconfig integer pert_1_y namelist,ideal max_domains 0 - "pert_1_y" "" "" +rconfig real pert_1_z namelist,ideal max_domains 1500. - "pert_1_z" "" "" +rconfig real pert_1_delt namelist,ideal max_domains 3. - "pert_1_delt" "" "" +rconfig real pert_1_width namelist,ideal max_domains 10000. - "pert_1_width" "" "" +rconfig integer pert_2_x namelist,ideal max_domains 0 - "pert_2_x" "" "" +rconfig integer pert_2_y namelist,ideal max_domains 0 - "pert_2_y" "" "" +rconfig real pert_2_z namelist,ideal max_domains 1500. - "pert_2_z" "" "" +rconfig real pert_2_width namelist,ideal max_domains 10000. - "pert_2_width" "" "" +rconfig real pert_2_delt namelist,ideal max_domains 3. - "pert_2_delt" "" "" +rconfig integer pert_3_x namelist,ideal max_domains 0 - "pert_3_x" "" "" +rconfig integer pert_3_y namelist,ideal max_domains 0 - "pert_3_y" "" "" +rconfig real pert_3_z namelist,ideal max_domains 1500. - "pert_3_z" "" "" +rconfig real pert_3_delt namelist,ideal max_domains 3. - "pert_3_delt" "" "" +rconfig real pert_3_width namelist,ideal max_domains 10000. - "pert_3_width" "" "" +rconfig integer pert_4_x namelist,ideal max_domains 0 - "pert_4_x" "" "" +rconfig integer pert_4_y namelist,ideal max_domains 0 - "pert_4_y" "" "" +rconfig real pert_4_z namelist,ideal max_domains 1500. - "pert_4_z" "" "" +rconfig real pert_4_delt namelist,ideal max_domains 3. - "pert_4_delt" "" "" +rconfig real pert_4_width namelist,ideal max_domains 10000. - "pert_4_width" "" "" +rconfig integer pert_5_x namelist,ideal max_domains 0 - "pert_5_x" "" "" +rconfig integer pert_5_y namelist,ideal max_domains 0 - "pert_5_y" "" "" +rconfig real pert_5_z namelist,ideal max_domains 1500. - "pert_5_z" "" "" +rconfig real pert_5_delt namelist,ideal max_domains 3. - "pert_5_delt" "" "" +rconfig real pert_5_width namelist,ideal max_domains 10000. - "pert_5_width" "" "" +rconfig integer pert_6_x namelist,ideal max_domains 0 - "pert_6_x" "" "" +rconfig integer pert_6_y namelist,ideal max_domains 0 - "pert_6_y" "" "" +rconfig real pert_6_z namelist,ideal max_domains 1500. - "pert_6_z" "" "" +rconfig real pert_6_delt namelist,ideal max_domains 3. - "pert_6_delt" "" "" +rconfig real pert_6_width namelist,ideal max_domains 10000. - "pert_6_width" "" "" +rconfig integer pert_7_x namelist,ideal max_domains 0 - "pert_7_x" "" "" +rconfig integer pert_7_y namelist,ideal max_domains 0 - "pert_7_y" "" "" +rconfig real pert_7_z namelist,ideal max_domains 1500. - "pert_7_z" "" "" +rconfig real pert_7_delt namelist,ideal max_domains 3. - "pert_7_delt" "" "" +rconfig real pert_7_width namelist,ideal max_domains 10000. - "pert_7_width" "" "" +rconfig integer pert_8_x namelist,ideal max_domains 0 - "pert_8_x" "" "" +rconfig integer pert_8_y namelist,ideal max_domains 0 - "pert_8_y" "" "" +rconfig real pert_8_z namelist,ideal max_domains 1500. - "pert_8_z" "" "" +rconfig real pert_8_delt namelist,ideal max_domains 3. - "pert_8_delt" "" "" +rconfig real pert_8_width namelist,ideal max_domains 10000. - "pert_8_width" "" "" +rconfig integer pert_9_x namelist,ideal max_domains 0 - "pert_9_x" "" "" +rconfig integer pert_9_y namelist,ideal max_domains 0 - "pert_9_y" "" "" +rconfig real pert_9_z namelist,ideal max_domains 1500. - "pert_9_z" "" "" +rconfig real pert_9_delt namelist,ideal max_domains 3. - "pert_9_delt" "" "" +rconfig real pert_9_width namelist,ideal max_domains 10000. - "pert_9_width" "" "" +rconfig integer pert_10_x namelist,ideal max_domains 0 - "pert_10_x" "" "" +rconfig integer pert_10_y namelist,ideal max_domains 0 - "pert_10_y" "" "" +rconfig real pert_10_z namelist,ideal max_domains 1500. - "pert_10_z" "" "" +rconfig real pert_10_delt namelist,ideal max_domains 3. - "pert_10_delt" "" "" +rconfig real pert_10_width namelist,ideal max_domains 10000. - "pert_10_width" "" "" #--------------------------------------------------------------------------------------------------------------------------------------- # Package Declarations # diff --git a/dyn_em/module_initialize_ideal.F b/dyn_em/module_initialize_ideal.F index 058b299233..b607fe2e7e 100644 --- a/dyn_em/module_initialize_ideal.F +++ b/dyn_em/module_initialize_ideal.F @@ -86,7 +86,7 @@ SUBROUTINE init_domain_rk ( grid & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & - i, j, k, kk + a, i, j, k, kk ! Local data @@ -96,11 +96,13 @@ SUBROUTINE init_domain_rk ( grid & INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc, lm + REAL :: div_mag, u_pert, v_pert REAL :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u - REAL :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2 + REAL :: z_scale, xrad, yrad, zrad, hgt, rad, delt, cof1, cof2 ! REAL, EXTERNAL :: interp_0 REAL :: hm, xa - REAL :: pi, rnd + REAL :: pi, rnd, scaled_rnd + INTEGER :: scaled_rand_int ! stuff from original initialization that has been dropped from the Registry REAL :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt @@ -113,6 +115,17 @@ SUBROUTINE init_domain_rk ( grid & LOGICAL :: stretch_grid, dry_sounding character (len=256) :: mminlu2 +! Multiple perturbation bubbles in em_quarter_ss + REAL, DIMENSION(10) :: pert_x, pert_y, pert_z, pert_dt, pert_width + INTEGER :: end_loop_a, internal_i, internal_j, isize,jsize + REAL, DIMENSION(2) :: fct1,fct2,fract_to_bdy + REAL :: qv_internal + + INTEGER :: seed_size + INTEGER, DIMENSION(:), ALLOCATABLE :: rng_seed + INTEGER :: k_cent + REAL :: min_zrad + REAL :: xa1, xal1,pii,hm1 ! data for intercomparison setup from dale ! space for initial jet in b_wave @@ -178,6 +191,9 @@ SUBROUTINE init_domain_rk ( grid & CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) + isize = ide-ids ! unstaggered size + jsize = jde-jds ! unstaggered size + ideal_constants: SELECT CASE ( model_config_rec%ideal_case ) CASE ( hill2d_x ) hm = 100. @@ -205,13 +221,89 @@ SUBROUTINE init_domain_rk ( grid & CASE ( quarter_ss, squall2d_x, squall2d_y ) stretch_grid = .true. - delt = 3. +! delt = 3. ! z_scale = .40 z_scale = 8000./config_flags%ztop pi = 2.*asin(1.0) write(6,*) ' pi is ',pi - nxc = (ide-ids)/2 - nyc = (jde-jds)/2 + IF (config_flags%num_perts > 0) THEN + ! The point here is to convert the namelist variables into arrays + pert_x(1) = config_flags%pert_1_x + pert_x(2) = config_flags%pert_2_x + pert_x(3) = config_flags%pert_3_x + pert_x(4) = config_flags%pert_4_x + pert_x(5) = config_flags%pert_5_x + pert_x(6) = config_flags%pert_6_x + pert_x(7) = config_flags%pert_7_x + pert_x(8) = config_flags%pert_8_x + pert_x(9) = config_flags%pert_9_x + pert_x(10) = config_flags%pert_10_x + pert_y(1) = config_flags%pert_1_y + pert_y(2) = config_flags%pert_2_y + pert_y(3) = config_flags%pert_3_y + pert_y(4) = config_flags%pert_4_y + pert_y(5) = config_flags%pert_5_y + pert_y(6) = config_flags%pert_6_y + pert_y(7) = config_flags%pert_7_y + pert_y(8) = config_flags%pert_8_y + pert_y(9) = config_flags%pert_9_y + pert_y(10) = config_flags%pert_10_y + pert_z(1) = config_flags%pert_1_z + pert_z(2) = config_flags%pert_2_z + pert_z(3) = config_flags%pert_3_z + pert_z(4) = config_flags%pert_4_z + pert_z(5) = config_flags%pert_5_z + pert_z(6) = config_flags%pert_6_z + pert_z(7) = config_flags%pert_7_z + pert_z(8) = config_flags%pert_8_z + pert_z(9) = config_flags%pert_9_z + pert_z(10) = config_flags%pert_10_z + pert_dt(1) = config_flags%pert_1_delT + pert_dt(2) = config_flags%pert_2_delT + pert_dt(3) = config_flags%pert_3_delT + pert_dt(4) = config_flags%pert_4_delT + pert_dt(5) = config_flags%pert_5_delT + pert_dt(6) = config_flags%pert_6_delT + pert_dt(7) = config_flags%pert_7_delT + pert_dt(8) = config_flags%pert_8_delT + pert_dt(9) = config_flags%pert_9_delT + pert_dt(10) = config_flags%pert_10_delT + pert_width(1) = config_flags%pert_1_width + pert_width(2) = config_flags%pert_2_width + pert_width(3) = config_flags%pert_3_width + pert_width(4) = config_flags%pert_4_width + pert_width(5) = config_flags%pert_5_width + pert_width(6) = config_flags%pert_6_width + pert_width(7) = config_flags%pert_7_width + pert_width(8) = config_flags%pert_8_width + pert_width(9) = config_flags%pert_9_width + pert_width(10) = config_flags%pert_10_width + end_loop_a = config_flags%num_perts + ELSE + end_loop_a = 1 + IF (config_flags%pert_1_x .gt. 0) THEN + pert_x(1) = INT(config_flags%pert_1_x) + ELSE + pert_x(1) = (ide-ids)/2 + ! nxc = (ide-ids)/2 + END IF + IF (config_flags%pert_1_y .gt. 0) THEN + pert_y(1) = INT(config_flags%pert_1_y) + ELSE + pert_y(1) = (jde-jds)/2 + ! nyc = (jde-jds)/2 + END IF + END IF + + ! Initialize random number generator using a fixed sequence that will be repeated every time the program runs + CALL random_seed(size=seed_size) + ALLOCATE(rng_seed(seed_size)) + IF (config_flags%line_seed < 0 ) THEN + rng_seed = 15571 + ELSE + rng_seed = config_flags%line_seed + END IF + CALL random_seed(put=rng_seed) CASE (grav2d_x) @@ -1028,7 +1120,7 @@ SUBROUTINE init_domain_rk ( grid & ! first from the top of the model to the top pressure kk = kte-1 ! top level k=kk+1 - qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk,j,P_QV)) + qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,kk,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%p(i,kk,j) = - 0.5*((grid%c1f(k)*grid%Mu_1(i,j))+qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/grid%rdnw(kk)/qvf2 @@ -1089,37 +1181,93 @@ SUBROUTINE init_domain_rk ( grid & CASE (quarter_ss) ! QSS thermal perturbation to kick off convection - write(6,*) ' nxc, nyc for perturbation ',nxc,nyc - write(6,*) ' delt for perturbation ',delt - - DO J = jts, min(jde-1,jte) - yrad = config_flags%dy*float(j-nyc)/10000. -! yrad = 0. - DO I = its, min(ide-1,ite) - xrad = config_flags%dx*float(i-nxc)/10000. -! xrad = 0. - DO K = 1, kte-1 - -! put in preturbation theta (bubble) and recalc density. note, + DO A = 1,end_loop_a + nyc = pert_y(a) + nxc = pert_x(a) + delt = pert_dt(a) + write(6,*) ' nxc, nyc for perturbation #',a,nxc,nyc + write(6,*) ' delt for perturbation #',a,delt + DO J = jts, min(jde-1,jte) + IF (config_flags%linear_bubble) THEN + yrad = 1.0 + IF (ABS(config_flags%dy*float(j-nyc)) .le. (0.5*config_flags%line_fract*((jde-jds)+1)*config_flags%dy)) yrad = 0.0 + ELSE + yrad = (config_flags%dy*float(j-nyc)) / pert_width(a) ! bubble width is pert_width(a) [m] + END IF + DO I = its, min(ide-1,ite) + IF (config_flags%linear_bubble) THEN + IF (config_flags%line_seed == 0) THEN + scaled_rnd = 0.0 + ELSE + CALL RANDOM_NUMBER(rnd) + scaled_rnd = 0.25*pert_width(a)*(rnd-0.5) + END IF + xrad = config_flags%dx*float(i-nxc)/(pert_width(a) + scaled_rnd) + ELSE + xrad = config_flags%dx*float(i-nxc)/pert_width(a) + END IF + k_cent = 0 + min_zrad = 9.E+10 + DO K = 1, kte-1 + +! put in perturbation theta (bubble) and recalc density. note, ! the mass in the column is not changing, so when theta changes, ! we recompute density and geopotential - zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j) & + hgt = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j) & +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g - zrad = (zrad-1500.)/1500. + IF (pert_z(a) >= 0.) THEN + zrad = (hgt-pert_z(a))/1500. ! Bubble has 1.5 km vertical radius, centered pert_z above ground + ELSE + zrad = (hgt-1500)/1500. ! Bubble has 1.5 km vertical radius, centered 1.5 km above ground + END IF + IF (ABS(zrad) .lt. min_zrad) THEN + min_zrad = ABS(zrad) + k_cent = k + END IF RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad) IF(RAD <= 1.) THEN - grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2 + IF (config_flags%linear_bubble) THEN + IF (config_flags%line_seed == 0) THEN + scaled_rnd = 0.0 + ELSE + CALL RANDOM_NUMBER(rnd) + scaled_rnd = 0.5*rnd - 0.5 + ENDIF + grid%t_1(i,k,j)=grid%t_1(i,k,j)+(delt+scaled_rnd)*COS(.5*PI*RAD)**2 + ELSE + grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2 + END IF grid%t_2(i,k,j)=grid%t_1(i,k,j) qvf = 1. + rvovrd*moist(i,k,j,P_QV) grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) - ENDIF - ENDDO + END IF + END DO ! k + + DO k = 1,kte-1 + ! Input an updraft core extending from the height of the bubble to the bottom of the model + ! Input a Gaussian shaped updraft embryo in the lower troposphere + IF (config_flags%init_vv) THEN + hgt = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j) + grid%phb(i,k,j)+grid%phb(i,k+1,j))/g + zrad = 0 + ! Set a maximum altitude of the vertical velocity adjustment to 10 km + IF (hgt <= 10000.) zrad = exp(-((hgt-pert_z(1))/config_flags%vv_lz)**2) + rad = exp(-(config_flags%dx*(i-pert_x(1))/config_flags%vv_ly)**2) * exp(-(config_flags%dy*(j-pert_y(1))/config_flags%vv_ly)**2) * zrad + grid%w_1(i,k,j) = config_flags%init_w_core*rad + grid%w_2(i,k,j) = config_flags%init_w_core*rad + grid%ww(i,k,j) = config_flags%init_w_core*rad + END IF + END DO ! k -! rebalance hydrostatically + END DO ! i + END DO ! j + END DO ! A - number of perturbation bubbles to create +! rebalance hydrostatically + DO J = jts, min(jde-1,jte) + DO I = its, min(ide-1,ite) DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & @@ -1131,6 +1279,7 @@ SUBROUTINE init_domain_rk ( grid & ENDDO ENDDO + CASE (squall2d_x) ! QSS thermal perturbation to kick off convection @@ -1554,13 +1703,37 @@ SUBROUTINE init_domain_rk ( grid & ENDDO ENDIF + ! Modify initial wind field to include a convergence component + IF (model_config_rec%ideal_case == quarter_ss .and. config_flags%init_div) THEN + div_mag = -0.5*config_flags%init_max_div * 1.0/(1./config_flags%div_lx**2 + 1./config_flags%div_ly**2) + DO K = kts,kte + DO J = jts,min(jde-1,jte) + DO I = its,min(ide-1,ite) + rad = & + exp(-(config_flags%dx*(i-pert_x(1))/config_flags%div_lx)**2) * & + exp(-(config_flags%dy*(j-pert_y(1))/config_flags%div_ly)**2) + hgt = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j) + grid%phb(i,k,j)+grid%phb(i,k+1,j))/g + zrad = MAX((config_flags%div_depth - hgt)/config_flags%div_depth,0.) + u_pert = -2*zrad*div_mag*(config_flags%dx*(i-pert_x(1))/config_flags%div_lx**2) * rad + v_pert = -2*zrad*div_mag*(config_flags%dy*(i-pert_y(1))/config_flags%div_ly**2) * rad + grid%u_1(i,k,j) = grid%u_1(i,k,j) + u_pert + grid%v_1(i,k,j) = grid%v_1(i,k,j) + v_pert + grid%u_2(i,k,j) = grid%u_1(i,k,j) + grid%v_2(i,k,j) = grid%v_1(i,k,j) + END DO + END DO + END DO + END IF + ! set w DO J = jts, min(jde-1,jte) DO K = kts, kte DO I = its, min(ide-1,ite) - grid%w_1(i,k,j) = 0. - grid%w_2(i,k,j) = 0. + IF (.not. config_flags%init_vv) THEN + grid%w_2(i,k,j) = 0. + grid%w_1(i,k,j) = 0. + END IF ENDDO ENDDO ENDDO diff --git a/dyn_em/module_stoch.F b/dyn_em/module_stoch.F index 8f1a6f13c2..d08a59bb24 100644 --- a/dyn_em/module_stoch.F +++ b/dyn_em/module_stoch.F @@ -168,7 +168,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_sppt, & grid%lengthscale_sppt, & grid%timescale_sppt, & @@ -186,7 +186,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_sppt, & grid%lengthscale_sppt, & grid%timescale_sppt, & @@ -210,7 +210,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_sppt, & grid%lengthscale_sppt, & grid%timescale_sppt, & @@ -234,7 +234,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_rand_pert, & grid%lengthscale_rand_pert, & grid%timescale_rand_pert, & @@ -396,7 +396,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%stepsp * grid%time_step,grid%DX,grid%DY, & + grid%stepsp * grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_mult3d(n), & grid%lengthscale_mult3d(n), & grid%timescale_mult3d(n), & @@ -460,7 +460,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_spp_conv, & grid%lengthscale_spp_conv, & grid%timescale_spp_conv, & @@ -483,7 +483,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_spp_pbl, & grid%lengthscale_spp_pbl, & grid%timescale_spp_pbl, & @@ -506,7 +506,7 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & - grid%time_step,grid%DX,grid%DY, & + grid%dt,grid%DX,grid%DY, & grid%gridpt_stddev_spp_lsm, & grid%lengthscale_spp_lsm, & grid%timescale_spp_lsm, & @@ -552,7 +552,9 @@ subroutine SETUP_RAND_PERTURB( variable_in,& INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte - INTEGER :: IER,IK,IL,I,J,itime_step,skebs_vertstruc, & + REAL, INTENT(IN) :: itime_step ! jdduda + ! jdduda ALL USES OF ITIME_STEP HAD BEEN ENCLOSED IN FLOAT(). THESE WERE REMOVED. + INTEGER :: IER,IK,IL,I,J,skebs_vertstruc, & ! jdduda removed itime_step from this line KMINFORCT,LMINFORCT,KMAXFORCT,LMAXFORCT,KMAXFORCTH,LMAXFORCTH, & KMAX,LMAX,LENSAV,ILEV REAL :: DX,DY,RY,RX,ALPH,RHOKLMAX,ZREF,RHOKL,EPS @@ -617,7 +619,8 @@ subroutine SETUP_RAND_PERTURB( variable_in,& ! -------------------------------------------------------------------------------------- ! ---------- INITIALIZE STOCHASTIC KINETIC-ENERGY BACKSCATTER SCHEME (SKEBS) ---------- ! -------------------------------------------------------------------------------------- - ALPH = float(itime_step)/ZTAU ! approximation of 1.-exp(-itime_step/ZTAU_PSI) +! ALPH = float(itime_step)/ZTAU ! approximation of 1.-exp(-itime_step/ZTAU_PSI) + ALPH = itime_step/ZTAU ! approximation of 1.-exp(-itime_step/ZTAU_PSI) ZSIGMA2=1./(12.0*ALPH) if (is_print) then @@ -633,7 +636,7 @@ subroutine SETUP_RAND_PERTURB( variable_in,& WRITE(*,'('' Minimal wavenumber of streamfunction forcing, KMINFORC ='',I10)') KMINFORCT WRITE(*,'('' Maximal wavenumber of streamfunction forcing, KMAXFORC ='',I10)') KMAXFORCT WRITE(*,'('' skebs_vertstruc '',I10)') skebs_vertstruc - WRITE(*,'('' Time step: itime_step='',I10)') itime_step + WRITE(*,'('' Time step: itime_step='',F10.3)') itime_step ! jdduda WRITE(*,'('' Decorrelation time of noise, ZTAU_PSI ='',E12.5)') ZTAU WRITE(*,'('' Variance of noise, ZSIGMA2_EPS ='',E12.5)') ZSIGMA2 WRITE(*,'('' Autoregressive parameter 1-ALPH_PSI ='',E12.5)') 1.-ALPH @@ -660,7 +663,7 @@ subroutine SETUP_RAND_PERTURB( variable_in,& IF ((variable == 'P') .or. (variable == 'R') .or. (variable == 'S') .or. (variable == 'Q') .or. (variable == 'O')) then kappat= L_rand_perturb**2 ! L^2= kappa*T, where L is a length scale in m; set to for L=100km - phi = exp (-float(itime_step)/tau_rand_perturb) + phi = exp (-itime_step/tau_rand_perturb) ! jdduda alph = 1.-phi endif @@ -761,9 +764,9 @@ subroutine SETUP_RAND_PERTURB( variable_in,& enddo ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled if (variable == 'W') then - ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT/(float(itime_step)*ZSIGMA2*ZGAMMAN))/(2*RPI) + ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT/(itime_step*ZSIGMA2*ZGAMMAN))/(2*RPI) elseif (variable == 'T') then - ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT/(float(itime_step)*cp*ZSIGMA2*ZGAMMAN)) + ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT/(itime_step*cp*ZSIGMA2*ZGAMMAN)) elseif ((variable == 'P') .or. (variable == 'R') .or. (variable == 'S') .or. (variable == 'O') .or. (variable == 'Q ')) then ZCONSTF0= gridpt_stddev_rand_perturb*sqrt((1.-phi**2)/(2.*ZGAMMAN)) endif diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 0761df036f..b8f4ab889f 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -3814,6 +3814,18 @@ END SUBROUTINE CMAQ_DRIVER ! Optional & , RAINNC=grid%rainnc, RAINNCV=grid%rainncv & & , SNOWNC=grid%snownc, SNOWNCV=grid%snowncv & + ! jdduda + , MAX_PRATE=grid%max_prate & + , MAX_PRATE_1MIN=grid%max_prate_1min & + , MAX_PRATE_5MIN=grid%max_prate_5min & + , MAX_PRATE_10MIN=grid%max_prate_10min & + , MAX_COMPREF=grid%max_compref & + , HISTORY_INTERVAL=config_flags%history_interval & + , HISTORY_INTERVAL_M=config_flags%history_interval_m & + , HISTORY_INTERVAL_S=config_flags%history_interval_s & + , ADAPTIVE_TS=config_flags%use_adaptive_time_step & + , MAX_TIME_STEP=grid%max_time_step & + ! jdduda & , GRAUPELNC=grid%graupelnc, GRAUPELNCV=grid%graupelncv & ! for milbrandt2mom & , HAILNC=grid%hailnc, HAILNCV=grid%hailncv & & , HAIL_MAXK1=grid%hail_maxk1,HAIL_MAX2D=grid%hail_max2d & diff --git a/phys/module_diag_misc.F b/phys/module_diag_misc.F index 54ba69e42d..d5817bc7e8 100644 --- a/phys/module_diag_misc.F +++ b/phys/module_diag_misc.F @@ -57,7 +57,7 @@ SUBROUTINE diagnostic_output_calc( & !----------- !-- DIAG_PRINT print control: 0 - no diagnostics; 1 - dmudt only; 2 - all !-- DT time step (second) -!-- XTIME forecast time +!-- XTIME forecast time (in minutes) !-- ACSWUPT !-- ACSWUPTC !-- ACSWDNT @@ -213,7 +213,8 @@ SUBROUTINE diagnostic_output_calc( & CHARACTER*6 :: grid_str INTEGER, INTENT(IN) :: & - history_interval,itimestep + history_interval, & ! units of sec + itimestep INTEGER :: idump @@ -226,7 +227,7 @@ SUBROUTINE diagnostic_output_calc( & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles - IF (mod(curr_secs2, history_interval*60.) == 0.) THEN + IF (mod(curr_secs2, FLOAT(history_interval)) == 0.) THEN WRITE(outstring,*) 'Reseting accumulation to 0' CALL wrf_debug ( 10, TRIM(outstring) ) DO j=j_start(ij),j_end(ij) @@ -299,7 +300,7 @@ SUBROUTINE diagnostic_output_calc( & DO ij = 1 , num_tiles - IF (xtime .eq. 0.0)THEN + IF (curr_secs2 .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 @@ -320,7 +321,7 @@ SUBROUTINE diagnostic_output_calc( & ENDDO ENDDO - IF (xtime .eq. 0.0 .and. bucket_J .gt. 0.0 .and. PRESENT(ACSWUPT))THEN + IF (curr_secs2 .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 @@ -334,7 +335,7 @@ SUBROUTINE diagnostic_output_calc( & ENDDO ENDDO ENDIF - IF (xtime .eq. 0.0 .and. bucket_J .gt. 0.0 .and. PRESENT(ACLWUPT))THEN + IF (curr_secs2 .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 @@ -463,7 +464,7 @@ SUBROUTINE diagnostic_output_calc( & if (diag_print .eq. 0 ) return - IF ( xtime .ne. 0. ) THEN + IF ( curr_secs2 .ne. 0. ) THEN ! COMPUTE THE NUMBER OF MASS GRID POINTS no_points = float((ide-ids)*(jde-jds)) @@ -626,7 +627,7 @@ SUBROUTINE diagnostic_output_calc( & ENDDO ENDDO - IF ( xtime .lt. 0.0001 ) THEN + IF ( curr_secs2 .lt. 0.0001 ) THEN DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) dpsdt(i,j)=0. diff --git a/phys/module_diag_nwp.F b/phys/module_diag_nwp.F index 336b0cd372..e67d9eb12e 100644 --- a/phys/module_diag_nwp.F +++ b/phys/module_diag_nwp.F @@ -62,7 +62,7 @@ SUBROUTINE diagnostic_output_nwp( & !----------- !-- DIAG_PRINT print control: 0 - no diagnostics; 1 - dmudt only; 2 - all !-- DT time step (second) -!-- XTIME forecast time +!-- XTIME forecast time (in minutes) !-- SBW specified boundary width - used later ! !-- P8W 3D pressure array at full eta levels @@ -153,7 +153,8 @@ SUBROUTINE diagnostic_output_nwp( & CHARACTER*6 :: grid_str INTEGER, INTENT(IN) :: & - history_interval,itimestep + history_interval, & ! units of sec + itimestep REAL, DIMENSION( kms:kme ), INTENT(IN) :: & znw @@ -217,13 +218,13 @@ SUBROUTINE diagnostic_output_nwp( & !----------------------------------------------------------------- - idump = (history_interval * 60.) / dt + idump = history_interval / dt ! IF ( MOD(itimestep, idump) .eq. 0 ) THEN ! WRITE(outstring,*) 'Computing PH0 for this domain with curr_secs2 = ', curr_secs2 ! CALL wrf_message ( TRIM(outstring) ) - time_from_output = mod( xtime, REAL(history_interval) ) + time_from_output = mod( xtime*60, REAL(history_interval) ) ! convert xtime from min to sec ! print *,' history_interval = ', history_interval ! print *,' itimestep = ', itimestep diff --git a/phys/module_diagnostics_driver.F b/phys/module_diagnostics_driver.F index aa583b505f..0357e9d098 100644 --- a/phys/module_diagnostics_driver.F +++ b/phys/module_diagnostics_driver.F @@ -171,6 +171,10 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & TYPE(WRFU_Time) :: currentTime + ! history interval converted to seconds + + INTEGER :: total_history_interval_sec + !============================================================= ! Start of executable code !============================================================= @@ -182,6 +186,16 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & k_start = kps k_end = kpe-1 + IF (config_flags%history_interval .NE. 0 .and. & + (config_flags%history_interval_d .ne. 0 .or. & + config_flags%history_interval_h .ne. 0 .or. & + config_flags%history_interval_m .ne. 0 .or. & + config_flags%history_interval_s .ne. 0 ) ) THEN + CALL wrf_message("NOTICE: *both* 'history_interval' and a specific history_interval type (e.g.,d,h,m,s) have been set to non-zero values.") + CALL wrf_message("NOTICE (cont): The physics diagnostics considers the specific types, such as 'history_interval_m' and 'history_interval_s', to take precedence over 'history_interval'") + END IF + total_history_interval_sec = 86400*config_flags%history_interval_d + 3600*config_flags%history_interval_h + 60*config_flags%history_interval_m + config_flags%history_interval_s + ! There are some fields that were defined in the first RK loop for ! physics and are now a time-step old. @@ -235,9 +249,11 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & !WRF-HAILCAST diagnostic - hail size prediction HAILCAST: IF ( config_flags%hailcast_opt /= 0 ) THEN - IF ( ( config_flags%history_interval == 0 ) ) THEN +! IF ( ( config_flags%history_interval == 0 ) ) THEN + IF ( total_history_interval_sec == 0) THEN WRITE (diag_message , * ) & - "HAILCAST Error : No 'history_interval' defined in namelist" +! "HAILCAST Error : No 'history_interval' defined in namelist" + "HAILCAST Error : total history_interval is 0" CALL wrf_error_fatal ( diag_message ) END IF @@ -389,7 +405,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,PREC_ACC_NC=grid%prec_acc_nc & ,PREC_ACC_DT=config_flags%prec_acc_dt & ,CURR_SECS2=curr_secs2 & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,cu_used=grid%cu_used & ,shcu_used=grid%shcu_used & @@ -431,7 +447,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,U10=grid%u10,V10=grid%v10,W=grid%w_2 & ,WSPD10MAX=grid%wspd10max & @@ -479,7 +495,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,U10=grid%u10,V10=grid%v10,W=grid%w_2 & ,WSPD10MAX=grid%wspd10max & @@ -529,7 +545,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,U10=grid%u10,V10=grid%v10,W=grid%w_2 & ,WSPD10MAX=grid%wspd10max & @@ -581,7 +597,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,U10=grid%u10,V10=grid%v10,W=grid%w_2 & ,WSPD10MAX=grid%wspd10max & @@ -631,7 +647,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,U10=grid%u10,V10=grid%v10,W=grid%w_2 & ,WSPD10MAX=grid%wspd10max & @@ -708,7 +724,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,CURR_SECS2=curr_secs2 & ,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics & ,DIAGFLAG=diag_flag & - ,HISTORY_INTERVAL=grid%history_interval & + ,HISTORY_INTERVAL=total_history_interval_sec & ,ITIMESTEP=grid%itimestep & ,U10=grid%u10,V10=grid%v10,W=grid%w_2 & ,WSPD10MAX=grid%wspd10max & @@ -948,9 +964,11 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & AFWA_DIAGS : IF ( config_flags%afwa_diag_opt == 1 ) THEN - IF ( ( config_flags%history_interval == 0 ) ) THEN +! IF ( ( config_flags%history_interval == 0 ) ) THEN + IF ( total_history_interval_sec == 0 ) THEN WRITE (diag_message , * ) & - "AFWA Diagnostics Error : No 'history_interval' defined in namelist" +! "AFWA Diagnostics Error : No 'history_interval' defined in namelist" + "AFWA Diagnostics Error : total history_interval is 0" CALL wrf_error_fatal ( diag_message ) END IF diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index 94004191cb..d3245d810c 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -100,6 +100,7 @@ SUBROUTINE microphysics_driver( & !NUWRF JJS 20110525 ^^^^^ ! ,ccntype & ! for mp_milbrandt2mom ,u,v,w,z & + ,max_prate,max_prate_1min,max_prate_5min,max_prate_10min,max_compref & ,rainnc, rainncv & ,snownc, snowncv & ,hailnc, hailncv & @@ -162,8 +163,10 @@ SUBROUTINE microphysics_driver( & ,perts_qsnow, perts_ni & ,pert_thom_qv,pert_thom_qc,pert_thom_qi & ,pert_thom_qs,pert_thom_ni & - ,cloudnc & - ) + ,cloudnc & + ,history_interval, history_interval_m, history_interval_s & + ,adaptive_ts, max_time_step & + ) ! Framework USE module_state_description, ONLY : & @@ -689,8 +692,13 @@ SUBROUTINE microphysics_driver( & ,GRAUPELNCV & ,HAILNC & ,HAILNCV & - ,CLOUDNC & - ,hail_maxk1, hail_max2d + ,CLOUDNC & + ,hail_maxk1, hail_max2d & + ,max_prate & + ,max_prate_1min & + ,max_prate_5min & + ,max_prate_10min & + ,max_compref #if ( WRF_CHEM == 1) ! NUWRF JJS 20110525 vvvvv @@ -773,7 +781,7 @@ SUBROUTINE microphysics_driver( & ! LOCAL VAR - INTEGER :: i,j,k,its,ite,jts,jte,ij,sz,n + INTEGER :: i,j,k,its,ite,jts,jte,ij,sz,n,i_mid,j_mid LOGICAL :: channel LOGICAL :: nssl_progn = .false. REAL :: z0, z1, z2, w1, w2 @@ -785,17 +793,26 @@ SUBROUTINE microphysics_driver( & REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, INTENT(IN):: IRRIGATION !ARI REAL, OPTIONAL, INTENT(IN):: irr_daily_amount, julian_in, xtime, gmt + INTEGER, OPTIONAL, INTENT(IN) :: history_interval, history_interval_m, history_interval_s INTEGER, OPTIONAL, INTENT(IN ):: sf_surf_irr_scheme, irr_start_hour, irr_num_hours, & irr_start_julianday,irr_end_julianday,irr_freq,irr_ph REAL, PARAMETER :: PI_GRECO=3.14159 - INTEGER :: end_hour,a,b,xt24,irr_day,timing + INTEGER :: end_hour,a,b,xt24,irr_day,timing, idump REAL :: constants_irrigation,tloc,irr_start,phase INTEGER, OPTIONAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: irr_rand_field + REAL :: time_from_output + INTEGER, INTENT(IN) :: max_time_step + LOGICAL, INTENT(IN) :: adaptive_ts + LOGICAL :: reset_arrays + ! To accommodate shared physics character*256 :: errmsg integer :: errflg + i_mid = (ids+ide)/2 + j_mid = (jds+jde)/2 + !--------------------------------------------------------------------- ! check for microphysics type. We need a clean way to ! specify these things! @@ -823,6 +840,31 @@ SUBROUTINE microphysics_driver( & CALL wrf_message ( wrf_err_message ) ENDIF + IF (PRESENT(history_interval_m) .and. PRESENT(history_interval_s)) THEN + idump = (history_interval_m * 60. + history_interval_s) / dt + time_from_output = mod( xtime, REAL(history_interval_m*60.+history_interval_s) ) + ELSE + idump = (history_interval * 60.) / dt + time_from_output = mod( xtime, REAL(history_interval) ) + END IF + + reset_arrays = .false. + IF ( adaptive_ts .EQV. .TRUE. ) THEN +! if timestep is adaptive, use new resetting method; +! also, we are multiplying max_time_step with 1.05; because of rounding error, +! the time_from_output can become slightly larger than max_time_step when +! adaptive time step is maximized, in which case the if condition below fails to detect +! that we just wrote an output + IF ( ( time_from_output .GT. 0 ) .AND. ( time_from_output .LE. ( ( max_time_step * 1.05 ) / 60. ) ) ) THEN + reset_arrays = .TRUE. + ENDIF + ELSE +! if timestep is fixed, use original resetting method + IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN + reset_arrays = .TRUE. + ENDIF + ENDIF + #ifdef XEON_OPTIMIZED_WSM5 ! the OpenMP loops are inside the scheme when running on MIC IF ( mp_physics .EQ. WSM5SCHEME ) THEN @@ -2902,6 +2944,38 @@ SUBROUTINE microphysics_driver( & END SELECT micro_select + IF (mp_physics == THOMPSON .or. & + mp_physics == THOMPSONGH .or. & + mp_physics == THOMPSONAERO .or. & + mp_physics == MORR_TWO_MOMENT .or. & + mp_physics == MORR_TM_AERO .or. & + mp_physics == MILBRANDT2MOM .or. & + mp_physics == NSSL_2MOM .or. & + mp_physics == LINSCHEME .or. & + mp_physics == P3_1CATEGORY .or. & + mp_physics == P3_1CATEGORY_NC .or. & + mp_physics == P3_2CATEGORY .or. & + mp_physics == P3_1CAT_3MOM .or. & + mp_physics == WSM5SCHEME .or. & + mp_physics == WSM6SCHEME .or. & + mp_physics == WSM7SCHEME .or. & + mp_physics == WDM5SCHEME .or. & + mp_physics == WDM6SCHEME .or. & + mp_physics == WDM7SCHEME .or. & + mp_physics == FER_MP_HIRES .or. & + mp_physics == FER_MP_HIRES_ADVECT) THEN + + IF (reset_arrays .eqv. .true.) THEN + max_prate = 0. + max_prate_1min = 0. + max_prate_5min = 0. + max_prate_10min = 0. + max_compref = -35. + END IF + call calc_prate_max(itimestep,xtime,i_mid,j_mid,its,ims,ite,ime,jts,jms,jte,jme,kts,kms,kte,kme,rainncv,rainnc,dt, & + max_prate,max_prate_1min,max_prate_5min,max_prate_10min,refl_10cm,max_compref) + END IF + ENDDO !$OMP END PARALLEL DO @@ -3001,4 +3075,155 @@ subroutine Remove_multi_perturb_mp_perturbations (perts_qvapor, perts_qcloud, pe end subroutine Remove_multi_perturb_mp_perturbations + subroutine calc_prate_max(timestep,sim_min,cent_x,cent_y,its,ims,ite,ime,jts,jms,jte,jme,kts,kms,kte,kme,rainncv,rainnc,time_step, & + max_prate,max_prate_1min,max_prate_5min,max_prate_10min,reflectivity_3D,max_compref) + + ! Subroutine inputs + INTEGER, INTENT(IN) :: its, ims, ite, ime, jts, jms, jte, jme, kts, kms, kte, kme, cent_x, cent_y, timestep + REAL, INTENT(IN) :: time_step + REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: max_prate, & + max_prate_1min, & + max_prate_5min, & + max_prate_10min, & + max_compref + REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: rainncv, rainnc + REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: reflectivity_3D + ! REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: rolling_precip_1min,rolling_precip_5min,rolling_precip_10min + REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: rolling_precip + + !local variables and arrays: + character(len=16) :: desc + integer:: i,j,k,z + integer :: n_sub_windows, max_steps_per_sub + integer, dimension(5) :: steps_per_sub, sub_window_len + logical, save :: init_rolling_precip_array = .false. + REAL, DIMENSION(5) :: n_dt_sub, sub_window_time + logical :: isopen + logical, save :: un1_open = .false. , un2_open = .false., un3_open = .false. + integer, save :: un1,un2,un3 + + n_sub_windows = 3 + sub_window_len(:n_sub_windows) = 60*(/1,5,10/) ! min -> sec + do z = 1,n_sub_windows + n_dt_sub(z) = sub_window_len(z)/time_step + ! Account for time steps that do not divide evenly into 1, 5, or 10 minutes + if (abs(fraction(n_dt_sub(z))) .gt. 1e-6) then + steps_per_sub(z) = nint(n_dt_sub(z)) + 1 + sub_window_time(z) = time_step*(steps_per_sub(z)-1) + else + steps_per_sub(z) = sub_window_len(z)/time_step + 1 + sub_window_time(z) = sub_window_len(z) + end if + if (steps_per_sub(z) == 1) then + write(desc,'(I3)') int(sub_window_len(z)/60.) + write(6,*) 'WARNING: model time step is coarser than '//trim(adjustl(desc))//' minute(s) and thus the ave_prate_[1-/5-/10-]min array will be either all zeros or meaningless.' + init_rolling_precip_array = .true. + end if + end do + max_steps_per_sub = maxval(steps_per_sub) + ! if all is okay, allocate rolling precip arrays just once and never again + if (.not. init_rolling_precip_array) then + allocate(rolling_precip(max_steps_per_sub,its:ite,jts:jte)) + rolling_precip = 0. + init_rolling_precip_array = .true. + end if + + DO i = its,ite + DO j = jts,jte + + max_prate(i,j) = max(rainncv(i,j)/time_step,max_prate(i,j)) + max_compref(i,j) = max(max_compref(i,j),MAXVAL(reflectivity_3D(i,:,j))) + + ! update rolling precipitation totals for sub-history-interval average-precipitation-rate calculation + IF (allocated(rolling_precip)) THEN + DO k = max_steps_per_sub-1,1,-1 + ! IF (i == cent_x .and. j == cent_y) PRINT '(A,I3,A,I3,A,F12.8,A,F12.8)', "rolling_precip(",k+1,") will be replaced with rolling_precip(",k,"), i.e, ", & + ! rolling_precip(k+1,cent_x,cent_y), " from ",rolling_precip(k,cent_x,cent_y) + rolling_precip(k+1,i,j) = rolling_precip(k,i,j) + END DO + ! Roll precip values + rolling_precip(1,i,j) = rainnc(i,j) + ! IF (i == cent_x .and. j == cent_y) PRINT *, "----------------------------------------------" + ! IF (i == cent_x .and. j == cent_y) PRINT '(A,F12.8,3X,F12.8)', "rolling_precip(1) is now ", rolling_precip(1,cent_x,cent_y), rainnc(cent_x,cent_y) + ! IF (i == cent_x .and. j == cent_y) PRINT '(A)', "----------------------------------------------" + ! Compute precipitation rates and assign to max values + max_prate_1min(i,j) = max((rolling_precip(1,i,j) - rolling_precip(steps_per_sub(1),i,j))/sub_window_time(1),max_prate_1min(i,j)) + max_prate_5min(i,j) = max((rolling_precip(1,i,j) - rolling_precip(steps_per_sub(2),i,j))/sub_window_time(2),max_prate_5min(i,j)) + max_prate_10min(i,j) = max((rolling_precip(1,i,j) - rolling_precip(steps_per_sub(3),i,j))/sub_window_time(3),max_prate_10min(i,j)) + END IF + + END DO + END DO + + IF (.not. un1_open) THEN + DO z = 50,1000 + INQUIRE(UNIT=z,OPENED=isopen) + IF (.not. isopen) THEN + un1 = z + EXIT + END IF + END DO + OPEN(unit=un1,file="precip_1min.txt",action='write',status='replace') + un1_open = .true. + END IF + IF (.not. un2_open) THEN + DO z = un1+1,1000 + INQUIRE(UNIT=z,OPENED=isopen) + IF (.not. isopen) THEN + un2 = z + EXIT + END IF + END DO + OPEN(unit=un2,file="precip_5min.txt",action='write',status='replace') + un2_open = .true. + END IF + IF (.not. un3_open) THEN + DO z = un2+1,1000 + INQUIRE(UNIT=z,OPENED=isopen) + IF (.not. isopen) THEN + un3 = z + EXIT + END IF + END DO + OPEN(unit=un3,file="precip_10min.txt",action='write',status='replace') + un3_open = .true. + END IF + + IF ( cent_x >= its .and. cent_x <= ite .and. cent_y >= jts .and. cent_y <= jte) THEN + inquire(unit=un1,opened=isopen) + IF (isopen) THEN + IF (timestep == 1) WRITE(un1,*) "grid point for data is ",cent_x,cent_y + WRITE(un1,'(A,F8.3,7(A,ES13.6))') & + "1-min rolling precip at time ",sim_min," [min]: newest - oldest = ", & + rolling_precip(1,cent_x,cent_y)," - ",rolling_precip(steps_per_sub(1),cent_x,cent_y), & + ", diff = ",rolling_precip(1,cent_x,cent_y) - rolling_precip(steps_per_sub(1),cent_x,cent_y), " [mm], & + &rate = ", 3600*(rolling_precip(1,cent_x,cent_y) - rolling_precip(steps_per_sub(1),cent_x,cent_y))/sub_window_time(1), " [mm/hr]), & + &(max: ", 3600*max_prate_1min(cent_x,cent_y), " [mm/hr]), & + &inst. value & inst. max: ",3600*rainncv(cent_x,cent_y)/time_step," / ",3600.*max_prate(cent_x,cent_y) + END IF + inquire(unit=un2,opened=isopen) + IF (isopen) THEN + IF (timestep == 1) WRITE(un2,*) "grid point for data is ",cent_x,cent_y + WRITE(un2,'(A,F8.3,7(A,ES13.6))') & + "5-min rolling precip at time ",sim_min," [min]: newest - oldest = ", & + rolling_precip(1,cent_x,cent_y)," - ",rolling_precip(steps_per_sub(1),cent_x,cent_y), & + ", diff = ",rolling_precip(1,cent_x,cent_y) - rolling_precip(steps_per_sub(2),cent_x,cent_y), " [mm], & + &rate = ", 3600*(rolling_precip(1,cent_x,cent_y) - rolling_precip(steps_per_sub(2),cent_x,cent_y))/sub_window_time(2), " [mm/hr]), & + &(max: ", 3600*max_prate_5min(cent_x,cent_y), " [mm/hr]), & + &inst. value & inst. max: ",3600*rainncv(cent_x,cent_y)/time_step," / ",3600.*max_prate(cent_x,cent_y) + END IF + inquire(unit=un3,opened=isopen) + IF (isopen) THEN + IF (timestep == 1) WRITE(un3,*) "grid point for data is ",cent_x,cent_y + WRITE(un3,'(A,F8.3,7(A,ES13.6))') & + "10-min rolling precip at time ",sim_min," [min]: newest - oldest = ", & + rolling_precip(1,cent_x,cent_y)," - ",rolling_precip(steps_per_sub(1),cent_x,cent_y), & + ", diff = ",rolling_precip(1,cent_x,cent_y) - rolling_precip(steps_per_sub(3),cent_x,cent_y), " [mm], & + &rate = ", 3600*(rolling_precip(1,cent_x,cent_y) - rolling_precip(steps_per_sub(3),cent_x,cent_y))/sub_window_time(3), " [mm/hr]), & + &(max: ", 3600*max_prate_10min(cent_x,cent_y), " [mm/hr]), & + &inst. value & inst. max: ",3600*rainncv(cent_x,cent_y)/time_step," / ",3600.*max_prate(cent_x,cent_y) + END IF + END IF + end subroutine calc_prate_max + END MODULE module_microphysics_driver diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index 551d4e9132..227946c083 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -2098,12 +2098,16 @@ END FUNCTION bep_bem_ngr_u END IF !----------------------------------------------------------------------- -! For nwp_diagnostics = 1, history_interval must be used. +! For nwp_diagnostics = 1, a non-zero history_interval must be used. !----------------------------------------------------------------------- IF ( ( model_config_rec%nwp_diagnostics .NE. 0 ) .AND. & - ( model_config_rec%history_interval(1) .EQ. 0 ) ) THEN - wrf_err_message = '--- ERROR: nwp_diagnostics requires the use of "history_interval" namelist.' + model_config_rec%history_interval(1) .EQ. 0 .AND. & + model_config_rec%history_interval_d(1) .EQ. 0 .AND. & + model_config_rec%history_interval_h(1) .EQ. 0 .AND. & + model_config_rec%history_interval_m(1) .EQ. 0 .AND. & + model_config_rec%history_interval_s(1) .EQ. 0 ) THEN + wrf_err_message = '--- ERROR: nwp_diagnostics requires a non-zero total history interval.' CALL wrf_message ( wrf_err_message ) wrf_err_message = '--- Replace interval variable with "history_interval".' CALL wrf_debug ( 0, TRIM( wrf_err_message ) )