diff --git a/README b/README index 09debb7a68..18c13efec1 100644 --- a/README +++ b/README @@ -1,4 +1,4 @@ -WRF Model Version 4.5 +WRF Model Version 4.5.1 https://www2.mmm.ucar.edu/wrf/users/ diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 50af9c465c..889a92854b 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -3275,7 +3275,7 @@ state real HAILCAST_DIAM_MEAN ij misc 1 - - state real HAILCAST_DIAM_STD ij misc 1 - - "HAILCAST_DIAM_STD" "WRF-HAILCAST Stand. Dev. Hail Diameter" "mm" state real HAILCAST_WUP_MASK ij misc 1 - r "HAILCAST_WUP_MASK" "Updraft mask, 1 if > 10m/s" "" state real HAILCAST_WDUR ij misc 1 - r "HAILCAST_WDUR" "Updraft duration" "sec" -state real haildtacttime - - - - r "haildtacttime" "HAILDTACTTIME" "HAILCAST ACTIVATION TIME in s" +state real haildtacttime ij misc 1 - r "haildtacttime" "HAILDTACTTIME" "HAILCAST ACTIVATION TIME in s" package hailcast hailcast_opt==1 - state:hailcast_diam_max,hailcast_diam_mean,hailcast_diam_std,hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5 diff --git a/arch/configure.defaults b/arch/configure.defaults index e7a9836226..aed98042a2 100644 --- a/arch/configure.defaults +++ b/arch/configure.defaults @@ -2054,6 +2054,7 @@ FCDEBUG = # -g $(FCNOOPT) # -fbacktrace -ggdb -fcheck=bounds,do,m FORMAT_FIXED = -ffixed-form FORMAT_FREE = -ffree-form -ffree-line-length-none FCSUFFIX = +FCCOMPAT = BYTESWAPIO = -fconvert=big-endian -frecord-marker=4 FCBASEOPTS_NO_G = -w $(FORMAT_FREE) $(BYTESWAPIO) $(FCCOMPAT) FCBASEOPTS = $(FCBASEOPTS_NO_G) $(FCDEBUG) diff --git a/configure b/configure index eaf4a1aa17..27b01b19ce 100755 --- a/configure +++ b/configure @@ -1043,7 +1043,7 @@ rm -f tools/rpc_test.log if [ -f tools/rpc_test.exe ] ; then rpc_type=`tools/rpc_test.exe` - if [ $rpc_type == "rpc" ]; then + if [ $rpc_type = "rpc" ]; then sed -e '/^CFLAGS_LOCAL/s/#/-DRPC_TYPES=1 &/' configure.wrf > configure.wrf.edit else sed -e '/^CFLAGS_LOCAL/s/#/-DRPC_TYPES=2 &/' configure.wrf > configure.wrf.edit diff --git a/external/io_adios2/wrf_io.F90 b/external/io_adios2/wrf_io.F90 index 99e9bfe2e0..3d5fdd6844 100644 --- a/external/io_adios2/wrf_io.F90 +++ b/external/io_adios2/wrf_io.F90 @@ -1599,9 +1599,9 @@ subroutine ext_adios2_ioinit(SysDepInfo, Status) !look for adios2 xml runtime configuration INQUIRE(FILE="adios2.xml", EXIST=file_exists) if(file_exists) then - call adios2_init(adios, 'adios2.xml', MPI_COMM_WORLD, adios2_debug_mode_on, stat) + call adios2_init(adios, 'adios2.xml', MPI_COMM_WORLD, stat) else - call adios2_init(adios, MPI_COMM_WORLD, adios2_debug_mode_on, stat) + call adios2_init(adios, MPI_COMM_WORLD, stat) endif call adios2_err(stat,Status) if(Status /= WRF_NO_ERR) then diff --git a/inc/version_decl b/inc/version_decl index 1c7a81726f..98f1ed2596 100644 --- a/inc/version_decl +++ b/inc/version_decl @@ -1 +1 @@ - CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.5' + CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.5.1' diff --git a/phys/module_bl_keps.F b/phys/module_bl_keps.F index 7f80a9fb6c..8c5be78733 100644 --- a/phys/module_bl_keps.F +++ b/phys/module_bl_keps.F @@ -1166,7 +1166,7 @@ subroutine calc_countergradient(dz,g,kms,kme,kts,kte,th0,te,d,tpe,Pra,b_t) dphidz(iz)=0.5*(dphidz1+dphidz2) enddo dphidz(kte)=0. - + dphidz(kts)=0. do iz=kts,kte b_t(iz)=b_t(iz)-dphidz(iz) enddo diff --git a/phys/module_diag_hailcast.F b/phys/module_diag_hailcast.F index 1798f919e8..098e9d9655 100644 --- a/phys/module_diag_hailcast.F +++ b/phys/module_diag_hailcast.F @@ -50,7 +50,8 @@ SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags & REAL, INTENT(IN ),OPTIONAL :: haildt REAL, INTENT(IN ),OPTIONAL :: curr_secs - REAL, INTENT(INOUT),OPTIONAL :: haildtacttime + REAL, DIMENSION( ims:ime, jms:jme ), & + INTENT(INOUT),OPTIONAL :: haildtacttime INTEGER, INTENT(IN ) :: itimestep REAL, INTENT(IN ) :: dt @@ -209,10 +210,18 @@ SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags & ( (.not. grid%nested) .or. & ( (grid%nested) .and. (abs(grid%dtbc) < 0.0001) ) ) )THEN doing_adapt_dt = .TRUE. - IF ( haildtacttime .eq. 0. ) THEN - haildtacttime = CURR_SECS + haildt - END IF - END IF + + !170519 - bug fix RAS + !Change haildtacttime to an array so can be set individually + ! by gridpoint to account for quilted processes + DO j = j_start, j_end + DO i = i_start, i_end + IF ( haildtacttime(i,j) .eq. 0. ) THEN + haildtacttime(i,j) = CURR_SECS + haildt + END IF + ENDDO + ENDDO + END IF ! Calculate STEPHAIL stephail = NINT(haildt / dt) @@ -270,12 +279,19 @@ SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags & decided = .TRUE. END IF - IF ( ( .NOT. decided ) .AND. & - ( doing_adapt_dt ) .AND. & - ( curr_secs .GE. haildtacttime ) ) THEN - run_param = .TRUE. - decided = .TRUE. - haildtacttime = curr_secs + haildt + !170519 - bug fix RAS + !Changed haildtacttime to an array so can be set individually + ! by gridpoint to account for quilted processes + IF ( ( .NOT. decided ) .AND. & + ( doing_adapt_dt ) .AND. & + ( curr_secs .GE. haildtacttime(i_start, j_start) ) ) THEN + run_param = .TRUE. + decided = .TRUE. + DO j = j_start, j_end + DO i = i_start, i_end + haildtacttime(i,j) = curr_secs + haildt + ENDDO + ENDDO END IF write ( message, * ) 'timestep to run HAILCAST? ', run_param @@ -442,6 +458,13 @@ END SUBROUTINE hailcast_diagnostic_driver !!!! 06 May 2017...................................Becky Adams-Selin AER !!!! Bug fixes for some memory issues in the adiabatic liquid !!!! water profile code. +!!!! 23 Apr 2019...................................Becky Adams-Selin AER +!!!! Added check to make sure embryo spends at least some time in +!!!! cloud before returning a positive hail diameter. +!!!! 16 Dec 2022...................................Becky Adams-Selin AER +!!!! Consolodate a number of bug fixes, including fixing behavior +!!!! of haildt when tiling is used, and non-SI units errors in +!!!! HEATBUD subroutine. !!!! !!!! See Adams-Selin and Ziegler 2016, MWR for further documentation. !!!! @@ -609,7 +632,8 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,& !the difference between the saturation mixing ratio at each level !and the cloud base water vapor mixing ratio DO k=1,nz - IF (k.LT.KBAS) THEN + RWA_new(k) = RWA(k) + IF (k.LT.KBAS) THEN RWA_adiabat(k) = 0. CYCLE ENDIF @@ -640,19 +664,10 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,& ELSE IF ((k.ge.KBAS+1).AND.(RWA_adiabat(k).ge.1.E-12)) THEN RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1) IF (RWA_new(k).LT.0) RWA_new(k) = 0. - ELSE - RWA_new(k) = RWA(k) ENDIF - ! - old code - ! - !IF (k.eq.KBAS+1) THEN - ! RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - !ELSE IF ((k.ge.KBAS+2).AND.(RWA_adiabat(k).ge.1.E-12)) THEN - ! RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1) - ! IF (RWA_new(k).LT.0) RWA_new(k) = 0. - !ENDIF ENDDO !remove the height factor from RWA_new - DO k=KBAS,nz + DO k=KBAS+1,nz !bug fix, ensure index start 1 higher RWA_new(k) = RWA_new(k) / (h1d(k)-h1d(k-1)) ENDDO @@ -696,6 +711,11 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,& CALL INTERP(VUU_pert, VUFZL, WFZLP, IFOUT, PA, nz) CALL INTERP(rho1d, DENSAFZL, WFZLP, IFOUT, PA, nz) + !Check if height of cloud base is above embryo insertion point + !RAS-20190423-! + IF (ZFZL < ZBAS-1000) GOTO 400 + + !Begin hail simulation time (seconds) sec = 0. @@ -850,9 +870,15 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,& !so let's shed any excess water now. D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 D = D_ICE - CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) + CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT,TS,DENSE) ENDIF !end check for melting call + + !RAS - 20190423 - ! + 400 CONTINUE !Check for if embryo insertion point is below cloud base + + !Require embryo to have stayed aloft for at least some time + IF (sec.LT.60) D = 0. !Check to make sure something hasn't gone horribly wrong IF (D.GT.0.254) D = 0. !just consider missing for now if > 10 in @@ -993,7 +1019,7 @@ SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) RE=(GX/0.6)**0.5 !SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON - !THE BEST NUMBER + !THE BEST NUMBER. Follows RH87 pg. 2762, but fixes an error in their (B1). IF (GX.LT.550) THEN W=LOG10(GX) Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0) @@ -1115,7 +1141,7 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !experimentally derived, and are expecting cal-C-g units. Do some conversions. DENSAC = DENSA * (1.E3) * (1.E-6) !dynamic viscosity - ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5 + ANU=1.717E-5*(393.0/(TC+120.0))*(TC/273.155)**1.5 !RAS units fix to agree with Roger and Yau ps. 102 !!! CALCULATE THE REYNOLDS NUMBER - unitless RE=D*VT*DENSAC/ANU E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv) @@ -1123,9 +1149,10 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& IF(RE.LT.6000.0)THEN AE=0.78+0.308*E ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN - AE=0.76*E + !RAS 190713 + AE=0.5*0.76*E !bug fix to match RasmussenHeymsfield1987 ELSEIF(RE.GE.20000.0) THEN - AE=(0.57+9.0E-6*RE)*E + AE = 0.5*(0.57+9.0E-6*RE)*E !bug fix to match RasmussenHeymsfield1987 ENDIF @@ -1159,6 +1186,9 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !!! CALCULATE THE TOTAL MASS CHANGE DGM=DGMW+DGMI+DGMV + !Dummy arguments in the event of no water, vapor, or ice growth + DENSELW = 900. + DENSELI = 700. !!! CALCULATE DENSITY OF NEW LAYER, DEPENDS ON FW AND ITYPE IF (ITYPE.EQ.1) THEN !DRY GROWTH !If hailstone encountered supercooled water, calculate new layer density @@ -1169,7 +1199,12 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !!! FIND THE STOKES NUMBER (rasmussen heymsfield 1985) NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50) !need hail radius in cm !!! FIND IMPACT VELOCITY (rasmussen heymsfield 1985) - W = LOG10(NS) + !W = LOG10(NS) - bug fix 20221216 - RAS ! + IF (NS.LT.0.1)THEN + W=-1. + ELSE + W = LOG10(NS) + ENDIF IF (RE.GT.200) THEN IF (NS.LT.0.1) THEN VIMP = 0.0 @@ -1213,7 +1248,7 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& AFACTOR = -DC*VIMP/TSCELSIUS IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN DENSELW = 0.30*(AFACTOR)**0.44 - ELSEIF (TSCELSIUS.GT.-5.) THEN + ELSE DENSELW = EXP(-0.03115 - 1.7030*AFACTOR + & 0.9116*AFACTOR**2. - 0.1224*AFACTOR**3.) ENDIF @@ -1309,13 +1344,14 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & REAL RV, RD, G, PI, ALF, ALV, ALS, CI, CW, AK REAL H, AH, TCC, TSC, DELRWC, DENSAC, TDIFF + REAL DCM, GMC, DGMWC, DGMIC, GM1C, DGMC REAL DMLT REAL TSCm1, TSCm2 DATA RV/461.48/,RD/287.04/,G/9.78956/ DATA PI/3.141592654/,ALF/79.7/,ALV/597.3/ DATA ALS/677.0/,CI/0.5/,CW/1./ - !Convert values to non-SI units here + !Convert values to non-SI units here to match all the constants TSC = TS - 273.155 TSCm1 = TSm1 - 273.155 TSCm2 = TSm2 - 273.155 @@ -1323,12 +1359,18 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & DELRWC = DELRW * (1.E3) * (1.E-6) DENSAC = DENSA * (1.E3) * (1.E-6) !DI still in cm2/sec + DCM = D * 100. !m to cm + GMC = GM * 1000. !kg to g + GM1C = GM1 * 1000. !kg to g + DGMWC = DGMW * 1000. !kg to g + DGMIC = DGMI * 1000. !kg to g + DGMC = DGM * 1000. !kg to g !!! CALCULATE THE CONSTANTS AK=(5.8+0.0184*TCC)*1.E-5 !thermal conductivity - cal/(cm*sec*K) !dynamic viscosity - calculated in MASSAGR - !ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5 + !ANU=1.717E-5*(393.0/(TC+120.0))*(TC/273.155)**1.5 !!! CALCULATE THE REYNOLDS NUMBER - unitless !RE=D*VT*DENSAC/ANU - calculated in MASSAGR @@ -1341,10 +1383,11 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & AH=0.78+0.308*H !AE=0.78+0.308*E ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN - AH=0.76*H - !AE=0.76*E + !RAS 190713 + AH=0.5*0.76*H !bug fix to match RasmussenHeymsfield1987 ELSEIF(RE.GE.20000.0) THEN - AH=(0.57+9.0E-6*RE)*H + !RAS 190713 + AH=0.5*(0.57+9.0E-6*RE)*H !bug fix to match RasmussenHeymsfield1987 !AE=(0.57+9.0E-6*RE)*E calculated in MASSAGR ENDIF @@ -1358,11 +1401,12 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & !TSC=TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* & ! (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & ! DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC) - TSC=0.6*(TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* & - (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & - DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)) + & + !units fix + TSC=0.6*(TSC-TSC*DGMC/GM1C+SEKDEL/(GM1C*CI)* & + (2.*PI*DCM*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & + DGMWC/SEKDEL*(ALF+CW*TCC)+DGMIC/SEKDEL*CI*TCC)) + & 0.2*TSCm1 + 0.2*TSCm2 - + TS = TSC+273.155 IF (TS.GE.273.155) THEN TS=273.155 @@ -1375,14 +1419,21 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & IF (TCC.LT.0.) THEN !Original Hailcast algorithm - FW=FW-FW*DGM/GM1+SEKDEL/(GM1*ALF)* & - (2.*PI*D*(AH*AK*TCC-AE*ALV*DI*DELRWC)+ & - DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC) + !FW=FW-FW*DGM/GM1+SEKDEL/(GM1*ALF)* & + ! (2.*PI*D*(AH*AK*TCC-AE*ALV*DI*DELRWC)+ & + ! DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC) + !units fixed + FW = FW-FW*DGMC/GMC+SEKDEL/(GMC*ALF)* & + (2.*PI*DCM*(AH*AK*TCC-AE*ALV*DI*DELRWC)+ & + DGMWC/SEKDEL*(ALF+CW*TCC)+DGMIC/SEKDEL*CI*TCC) ELSE !Calculate decrease in ice mass due to melting - DMLT = (2.*PI*D*AH*AK*TCC + 2.*PI*D*AE*ALV*DI*DELRWC + & - DGMW/SEKDEL*CW*TCC) / ALF - FW = (FW*GM + DMLT) / GM + !!Bug fix 28Jul2021 RAS - non-SI units + !More non-SI units - fixed 20220307 RAS + ! Extra SEKDEL in the calculation - 20221102 RAS + DMLT = SEKDEL / ALF * (2.*PI*DCM*AH*AK*TCC + 2.*PI*DCM*AE*ALV*DI*DELRWC + & + DGMWC/SEKDEL*CW*TCC) + FW = (FW*GM1C + DMLT + DGMWC) / GMC ENDIF IF(FW.GT.1.)FW=1. @@ -1436,13 +1487,19 @@ SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT) END SUBROUTINE BREAKUP - SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) + SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT,TS,DENSE) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! This is a spherical hail melting estimate based on the Goyer !!! et al. (1969) eqn (3). The depth of the warm layer, estimated !!! terminal velocity, and mean temperature of the warm layer are !!! used. DRB. 11/17/2003. !!! + !!! Incorporated some known variables into the subroutine + !!! to use instead of constants (VT, TS, DENSE). RAS 10/2/2019 + !!! + !!! Note - this subroutine assumes melted water is immediately + !!! shed. Possibly not the most accurate assumption. + !!! !!! INPUT: TLAYER mean sub-cloud layer temperature (K) !!! PLAYER mean sub-cloud layer pressure (Pa) !!! RLAYER mean sub-cloud layer mixing ratio (kg/kg) @@ -1453,7 +1510,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) IMPLICIT NONE REAL*8 D - REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT + REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT, TS, DENSE REAL eenv, delta, ewet, de, der, wetold, wetbulb, wetbulbk REAL tdclayer, tclayer, eps, b, hplayer REAL*8 a @@ -1502,7 +1559,8 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) dv = 0.25e-4 ! diffusivity of water vapor (m2/s) pi = 3.1415927 rv = 1004. - 287. ! gas constant for water vapor - rhoice = 917.0 ! density of ice (kg/m**3) + !rhoice = 917.0 ! density of ice (kg/m**3) + rhoice = DENSE ! we know the stone density, let's use it r = D/2. ! radius of stone (m) !Compute residence time in warm layer @@ -1511,12 +1569,14 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) !Calculate dmdt based on eqn (3) of Goyer et al. (1969) !Reynolds number...from pg 317 of Atmo Physics (Salby 1996) !Just use the density of air at 850 mb...close enough. - rho = 85000./(287.*TLAYER) - re = rho*r*VT*.01/1.7e-5 + !rho = 85000./(287.*TLAYER) + rho = player/(287.*TLAYER) !we have the layer pressure, just use it + !units fix - r is now in meters, not mm. Plus we need D, not r. RAS 20220308 + re = rho*r*2*VT/1.7e-5 !Temperature difference between environment and hailstone surface - delt = wetbulb !- 0.0 !assume stone surface is at 0C - !wetbulb is in Celsius + !We know the stone surface temperature, let's use it + delt = wetbulb - (TS - 273.155) !again, wet bulb is in C !Difference in vapor density of air stream and equil vapor !density at the sfc of the hailstone @@ -1531,7 +1591,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) !Calculate new mass growth dmdt = (-1.7*pi*r*(re**0.5)/lf)*((ka*delt)+((lv-lf)*dv*dsig)) IF (dmdt.gt.0.) dmdt = 0 - mass = dmdt*tres + mass = dmdt*tres ! < 0 !Find the new hailstone diameter massorg = 1.33333333*pi*r*r*r*rhoice diff --git a/phys/module_ltng_cpmpr92z.F b/phys/module_ltng_cpmpr92z.F index f93f55c37d..4b1117f48c 100644 --- a/phys/module_ltng_cpmpr92z.F +++ b/phys/module_ltng_cpmpr92z.F @@ -6,10 +6,9 @@ ! Price, C., and D. Rind (1992), A Simple Lightning Parameterization for Calculating ! Global Lightning Distributions, J. Geophys. Res., 97(D9), 9919-9933, doi:10.1029/92JD00719. ! -! Wong, J., M. Barth, and D. Noone (2012), Evaluating a Lightning Parameterization -! at Resolutions with Partially-Resolved Convection, GMDD, in preparation. -! -! Contact: J. Wong +! Wong, J., M. Barth, and D. Noone, 2013: Evaluating a lightning parameterization +! based on cloud-top height for mesoscale numerical model simulations, +! Geosci. Model Dev., 6, 429–443, https://doi.org/10.5194/gmd-6-429-2013. ! !********************************************************************** diff --git a/phys/module_ltng_crmpr92.F b/phys/module_ltng_crmpr92.F index eb9ba7c07a..986e1bb53a 100644 --- a/phys/module_ltng_crmpr92.F +++ b/phys/module_ltng_crmpr92.F @@ -6,15 +6,14 @@ ! Price, C., and D. Rind (1992), A Simple Lightning Parameterization for Calculating ! Global Lightning Distributions, J. Geophys. Res., 97(D9), 9919-9933, doi:10.1029/92JD00719. ! -! Wong, J., M. Barth, and D. Noone (2012), Evaluating a Lightning Parameterization -! at Resolutions with Partially-Resolved Convection, GMDD, in preparation. +! Wong, J., M. Barth, and D. Noone, 2013: Evaluating a lightning parameterization +! based on cloud-top height for mesoscale numerical model simulations, +! Geosci. Model Dev., 6, 429–443, https://doi.org/10.5194/gmd-6-429-2013. ! ! Unlike previous implementation, this version will produce slightly inconsistent ! IC and CG grid-flash rates against NO emission after production via calling ! lightning_nox_decaria. ! -! Contact: J. Wong -! !********************************************************************** MODULE module_ltng_crmpr92 diff --git a/phys/module_ltng_iccg.F b/phys/module_ltng_iccg.F index d038f4a0e2..e4397419af 100644 --- a/phys/module_ltng_iccg.F +++ b/phys/module_ltng_iccg.F @@ -8,8 +8,6 @@ ! ! See comments preceeding each method for details ! -! Contact: J. Wong -! !********************************************************************** diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index 7d9c6cbefb..d283b2f7eb 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -837,7 +837,7 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, ng, & !..Create bins of cloud ice (from min diameter up to 5x min snow size). xDx(1) = D0i*1.0d0 - xDx(nbi+1) = 5.0d0*D0s + xDx(nbi+1) = 2.0d0*D0s do n = 2, nbi xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) & *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1))) @@ -2542,7 +2542,7 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & png_rcg(k) = MIN(DBLE(ng(k)*odts), png_rcg(k)) pbg_rcg(k) = prg_rcg(k)/rho_g(idx_bg(k)) !..Put in explicit drop break-up due to collisions. - pnr_rcg(k) = -5.*tnr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r) ! RAIN2M + pnr_rcg(k) = -1.5*tnr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r) ! RAIN2M endif endif endif diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 7b1487ed56..e26df70a7d 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -5528,8 +5528,14 @@ subroutine aerosol_in_cu(aeromcu,alevsiz,no_months,no_src_types,XLAT,XLONG,aerop READ(UNIT=iunit) lon_aer READ(UNIT=iunit) aeropin READ(UNIT=iunit) aeroin + ELSE + CALL wrf_error_fatal('Unable to open file "CESM_RCP4.5_Aerosol_Data.dat" (required for aercu_opt = 1 or 2). ' // & + 'Please refer to the run/README.namelist file for more information.') END IF CLOSE(iunit) + ELSE + CALL wrf_error_fatal('Missing file "CESM_RCP4.5_Aerosol_Data.dat" (required for aercu_opt = 1 or 2). ' // & + 'Please refer to the run/README.namelist file for more information.') END IF !--------------------------------------------------------------------------- !-- latitudinally interpolate ozone data (and extend longitudinally) diff --git a/phys/module_sf_bem.F b/phys/module_sf_bem.F index 2a632fa67b..e693e7c6d2 100644 --- a/phys/module_sf_bem.F +++ b/phys/module_sf_bem.F @@ -1,4 +1,11 @@ MODULE module_sf_bem +#if defined(mpas) +use mpas_atmphys_utilities, only: physics_error_fatal +#define FATAL_ERROR(M) call physics_error_fatal( M ) +#else +use module_wrf_error +#define FATAL_ERROR(M) call wrf_error_fatal( M ) +#endif ! ----------------------------------------------------------------------- ! Variables and constants used in the BEM module ! ----------------------------------------------------------------------- @@ -2459,7 +2466,7 @@ subroutine gaussjbem(a,n,b,np) icol=k endif elseif(ipiv(k).gt.1)then - CALL wrf_error_fatal('singular matrix in gaussjbem') + FATAL_ERROR('singular matrix in gaussjbem') endif enddo endif @@ -2480,7 +2487,7 @@ subroutine gaussjbem(a,n,b,np) endif - if(a(icol,icol).eq.0) CALL wrf_error_fatal('singular matrix in gaussjbem') + if(a(icol,icol).eq.0) FATAL_ERROR('singular matrix in gaussjbem') pivinv=1./a(icol,icol) a(icol,icol)=1 diff --git a/phys/module_sf_bep_bem.F b/phys/module_sf_bep_bem.F index a9bcf50d99..f249700a08 100644 --- a/phys/module_sf_bep_bem.F +++ b/phys/module_sf_bep_bem.F @@ -3722,7 +3722,7 @@ subroutine gaussj(a,n,b,np) icol=k endif elseif(ipiv(k).gt.1)then - CALL wrf_error_fatal('singular matrix in gaussj') + FATAL_ERROR('singular matrix in gaussj') endif enddo endif @@ -3743,7 +3743,7 @@ subroutine gaussj(a,n,b,np) endif - if(a(icol,icol).eq.0) CALL wrf_error_fatal('singular matrix in gaussj') + if(a(icol,icol).eq.0) FATAL_ERROR('singular matrix in gaussj') pivinv=1./a(icol,icol) a(icol,icol)=1 diff --git a/phys/module_sf_sfclayrev.F b/phys/module_sf_sfclayrev.F index 34ff70b7a4..9f65730122 100644 --- a/phys/module_sf_sfclayrev.F +++ b/phys/module_sf_sfclayrev.F @@ -1212,7 +1212,9 @@ function zolri(ri,z,z0) ! fx1=zolri2(x1,ri,z,z0) fx2=zolri2(x2,ri,z,z0) + iter = 0 Do While (abs(x1 - x2) > 0.01) + if (iter .eq. 10) return ! check added for potential divide by zero (2019/11) if(fx1.eq.fx2)return if(abs(fx2).lt.abs(fx1))then @@ -1225,6 +1227,7 @@ function zolri(ri,z,z0) zolri=x2 endif ! + iter = iter + 1 enddo ! diff --git a/phys/module_sf_urban.F b/phys/module_sf_urban.F index 6a7eec58cc..f7c04091a4 100644 --- a/phys/module_sf_urban.F +++ b/phys/module_sf_urban.F @@ -2161,9 +2161,9 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & ALLOCATE( HSESF_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) FATAL_ERROR('Error allocating HSESF_TBL in urban_param_init') ALLOCATE( PV_FRAC_ROOF_TBL(ICATE), stat=allocate_status ) - if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating PV_FRAC_ROOF_TBL in urban_param_init') + if(allocate_status /= 0) FATAL_ERROR('Error allocating PV_FRAC_ROOF_TBL in urban_param_init') ALLOCATE( GR_FRAC_ROOF_TBL(ICATE), stat=allocate_status ) - if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating GR_FRAC_ROOF_TBL in urban_param_init') + if(allocate_status /= 0) FATAL_ERROR('Error allocating GR_FRAC_ROOF_TBL in urban_param_init') endif numdir_tbl = 0 @@ -2889,7 +2889,7 @@ SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, IF (CHECK.EQ.0)THEN IF(IVGTYP(I,J).EQ.1)THEN write(mesg,*) 'Sample of Urban settings' - call wrf_message(mesg) + WRITE_MESSAGE(mesg) write(mesg,*) 'TSURFACE0_URB',TSURFACE0_URB(I,J) WRITE_MESSAGE(mesg) write(mesg,*) 'TDEEP0_URB', TDEEP0_URB(I,J) diff --git a/phys/noahmp b/phys/noahmp index 3be0b2860d..981d4f859c 160000 --- a/phys/noahmp +++ b/phys/noahmp @@ -1 +1 @@ -Subproject commit 3be0b2860dab167006a0b3c4822e234ca253c3df +Subproject commit 981d4f859ce6c64213d38a783654c05b47b3485e diff --git a/run/README.namelist b/run/README.namelist index 6e7ebc1a3d..4efccbe253 100644 --- a/run/README.namelist +++ b/run/README.namelist @@ -902,6 +902,11 @@ Namelist variables for controlling the adaptive time step option: = 0: no aerosol interaction = 1: aerosol interaction with MSKF only = 2: aerosol interaction with both MSKF and Morrison + ! aercu_opt = 1 and 2 both require the "CESM_RCP4.5_Aerosol_Data.dat" static runtime data + file. This can be downloaded from: + http://www2.mmm.ucar.edu/wrf/src/wrf_files/CESM_RCP4.5_Aerosol_Data.tar.gz + Once unpacked link/copy either of the two files to CESM_RCP4.5_Aerosol_Data.dat in + your run directory prior to running wrf.exe aercu_fct factor to multiply with aerosol amount = 1. ! default value no_src_types_cu = 1 ! number of aerosol species in global aerosol data: 10 for CESM input, diff --git a/run/URBPARM_LCZ.TBL b/run/URBPARM_LCZ.TBL index ad401d7977..80e6809c17 100644 --- a/run/URBPARM_LCZ.TBL +++ b/run/URBPARM_LCZ.TBL @@ -543,7 +543,11 @@ BUILDING HEIGHTS: 1 # height Percentage # [m] [%] - 21.0 100.0 + 40.0 17. + 45.0 21. + 50.0 24. + 55.0 21. + 60.0 17. END BUILDING HEIGHTS BUILDING HEIGHTS: 2 @@ -551,10 +555,8 @@ BUILDING HEIGHTS: 2 # height Percentage # [m] [%] - 9.0 25.0 - 15.0 45.0 - 18.0 20.0 - 21.0 10.0 + 15.0 50. + 20.0 50. END BUILDING HEIGHTS BUILDING HEIGHTS: 3 @@ -562,8 +564,8 @@ BUILDING HEIGHTS: 3 # height Percentage # [m] [%] - 6.0 55.0 - 9.0 45.0 + 5.0 72. + 10.0 28. END BUILDING HEIGHTS BUILDING HEIGHTS: 4 @@ -571,7 +573,11 @@ BUILDING HEIGHTS: 4 # height Percentage # [m] [%] - 24.0 100.0 + 40.0 17. + 45.0 21. + 50.0 24. + 55.0 21. + 60.0 17. END BUILDING HEIGHTS BUILDING HEIGHTS: 5 @@ -579,10 +585,8 @@ BUILDING HEIGHTS: 5 # height Percentage # [m] [%] - 9.0 10.0 - 15.0 25.0 - 21.0 40.0 - 24.0 25.0 + 15.0 50. + 20.0 50. END BUILDING HEIGHTS BUILDING HEIGHTS: 6 @@ -590,9 +594,8 @@ BUILDING HEIGHTS: 6 # height Percentage # [m] [%] - 6.0 30.0 - 9.0 40.0 - 15.0 30.0 + 5.0 72. + 10.0 28. END BUILDING HEIGHTS BUILDING HEIGHTS: 7 @@ -600,11 +603,7 @@ BUILDING HEIGHTS: 7 # height Percentage # [m] [%] - 5.0 100.0 - 9.0 0.0 - 15.0 0.0 - 21.0 0.0 - 24.0 0.0 + 5.0 100. END BUILDING HEIGHTS BUILDING HEIGHTS: 8 @@ -612,8 +611,8 @@ BUILDING HEIGHTS: 8 # height Percentage # [m] [%] - 6.0 35.0 - 9.0 65.0 + 5.0 72. + 10.0 28. END BUILDING HEIGHTS BUILDING HEIGHTS: 9 @@ -621,8 +620,8 @@ BUILDING HEIGHTS: 9 # height Percentage # [m] [%] - 6.0 75.0 - 9.0 25.0 + 5.0 72. + 10.0 28. END BUILDING HEIGHTS BUILDING HEIGHTS: 10 @@ -630,9 +629,9 @@ BUILDING HEIGHTS: 10 # height Percentage # [m] [%] - 3.0 10.0 - 9.0 60.0 - 15.0 30.0 + 5.0 16. + 10.0 68. + 15.0 16. END BUILDING HEIGHTS BUILDING HEIGHTS: 11 @@ -640,8 +639,7 @@ BUILDING HEIGHTS: 11 # height Percentage # [m] [%] - 6.0 100.0 - + 5.0 100.0 END BUILDING HEIGHTS diff --git a/test/em_real/namelist.input.default b/test/em_real/namelist.input.default deleted file mode 100644 index 31597f5ddc..0000000000 --- a/test/em_real/namelist.input.default +++ /dev/null @@ -1,102 +0,0 @@ - &time_control - run_days = 0, ! Number of simulation days. DO NOT ADD ADDITIONAL COLUMNS. - run_hours = 36, ! Number of simulation hours. DO NOT ADD ADDITIONAL COLUMNS. - run_minutes = 0, ! Number of simulation minutes. DO NOT ADD ADDITIONAL COLUMNS. - run_seconds = 0, ! Number of simulation seconds. DO NOT ADD ADDITIONAL COLUMNS. - ! *NOTE: real.exe does not read run_days/hours/minutes/seconds. Only read by wrf.exe - ! real.exe only reads start_* and end_* parameters. - start_year = 2019, 2019, 2019, ! Starting year of the simulation. - start_month = 09, 09, 09, ! Starting month of the simulation. - start_day = 04, 04, 24, ! Starting day of the simulation. - start_hour = 12, 12, 12, ! Starting hour of the simulation. - end_year = 2019, 2019, 2019, ! Ending year of the simulation. - end_month = 09, 09, 09, ! Ending month of the simulation. - end_day = 06, 06, 06, ! Ending day of the simulation. - end_hour = 00, 00, 00, ! Ending hour of the simulation. - interval_seconds = 21600 ! Time interval (s) between input data. DO NOT ADD ADDITIONAL COLUMNS. - input_from_file = .true.,.true.,.true., ! Logical - whether input files (e.g., met_em*) will be used. - history_interval = 180, 60, 60, ! Interval between history output (mins). *Note: to change the units of history_interval, add _s (secs), _h (hours) - frames_per_outfile = 1000, 1000, 1000, ! Number of history output times printed to each wrf output (wrfout) file. - restart = .false., ! Logical - whether the run is a "restart" simulation. DO NOT ADD ADDITIONAL COLUMNS. - restart_interval = 7200, ! Interval between restart (wrfrst) files (min). DO NOT ADD ADDITIONAL COLUMNS. - io_form_history = 2 ! File format for history (wrfout) files. See Users' Guide for options. DO NOT ADD ADDITIONAL COLUMNS. - io_form_restart = 2 ! File format for restart (wrfrst) files. See Users' Guide for options. DO NOT ADD ADDITIONAL COLUMNS. - io_form_input = 2 ! File format for the input (wrfinput) files. See Users' Guide for options. DO NOT ADD ADDITIONAL COLUMNS. - io_form_boundary = 2 ! File format for the boundary (wrfbdy) file. See Users' Guide for options. DO NOT ADD ADDITIONAL COLUMNS. - / - - &domains - time_step = 90, ! Time step for model integration (s). Must use 6xDX (or less). DO NOT ADD ADDITIONAL COLUMNS. - time_step_fract_num = 0, ! Numerator if using fractional time step. DO NOT ADD ADDITIONAL COLUMNS. - time_step_fract_den = 1, ! Denomenator if using fractional time step. DO NOT ADD ADDITIONAL COLUMNS. - max_dom = 1, ! Total number of domains for the simulation. DO NOT ADD ADDITIONAL COLUMNS. - e_we = 150, 220, 200, ! Number of grid spaces in the x direction (west-east). - e_sn = 130, 214, 210, ! Number of grid spaces in the y direction (south-north). - e_vert = 45, 45, 45, ! Number of full vertical levels. must be the same for all domains. - dzstretch_s = 1.1 ! Surface stretch factor when auto_levels_opt = 2 (which is default). DO NOT ADD ADDITIONAL COLUMNS. - p_top_requested = 5000, ! Pressure top (in Pa) to use for the simulation. The input data must go up to this value. DO NOT ADD ADDITIONAL COLUMNS. - num_metgrid_levels = 34, ! Number of vertical levels in the WPS output. DO NOT ADD ADDITIONAL COLUMNS. - num_metgrid_soil_levels = 4, ! Number of soil levels in the WPS output. DO NOT ADD ADDITIONAL COLUMNS. - dx = 15000, ! Length (m) of each grid in the x (west-east) direction (resolution). DO NOT ADD ADDITIONAL COLUMNS. - dy = 15000, ! Length (m) of each grid in the y (south-north) direction (resolution). DO NOT ADD ADDITIONAL COLUMNS. - grid_id = 1, 2, 3, ! Identifier number of each domain. - parent_id = 0, 1, 2, ! The ID of the parent of each domain. - i_parent_start = 1, 53, 30, ! Starting location of the lower left corner i-indice within the parent domain. - j_parent_start = 1, 25, 30, ! Starting location of the lower left corner j-indice within the parent domain. - parent_grid_ratio = 1, 3, 3, ! Ratio of the parent to nest grid size (resolution). - parent_time_step_ratio = 1, 3, 3, ! Time_step ratio of the parent to nest. - feedback = 1, ! Switch for feedback from the nest to it's parent at corresponding grid points (1=on;0=off). DO NOT ADD ADDITIONAL COLUMNS. - smooth_option = 0 ! Smoothing option for the boundary between parent and nest domains (when feedback=1). DO NOT ADD ADDITIONAL COLUMNS. - / - - &physics - physics_suite = 'CONUS' ! Physics suite option (CONUS or tropical). See Users' Guide for details. DO NOT ADD ADDITIONAL COLUMNS. - mp_physics = -1, -1, -1, ! Microphysics scheme. See Users' Guide for options (-1=option from suite). Every domain must be the same. - cu_physics = -1, -1, 0, ! Cumulus scheme. See Users' Guide for options (-1=option from suite). - ra_lw_physics = -1, -1, -1, ! Longwave radiation scheme. See Users' Guide for options (-1=option from suite). Every domain must be the same. - ra_sw_physics = -1, -1, -1, ! Shortwave radiation scheme. See Users' Guide for options (-1=option from suite). Every domain must be the same. - bl_pbl_physics = -1, -1, -1, ! PBL scheme. See Users' Guide for options (-1=option from suite). - sf_sfclay_physics = -1, -1, -1, ! Surface layer scheme. See Users' Guide for options (-1=option from suite). Every domain must be the same. - sf_surface_physics = -1, -1, -1, ! LSM scheme. See Users' Guide for options (-1=option from suite). Every domain must be the same. - radt = 15, 15, 15, ! Minutes between radiation physics calls (recommended 1 min per km of d01 dx). Every domain must be the same. - bldt = 0, 0, 0, ! Minutes between boundary-layer physics calls (recommended 0=every time step). Every domain must be the same. - cudt = 0, 0, 0, ! Minutes between cumulus physics calls. Set to 0 (called every time step) for all cumulus schemes except Kain-Fritsch. Every domain must be the same. - icloud = 1, ! Option to couple subgrid-scale clouds from the PBL scheme to radiation schemes (only works with RRTM/RRTMG rad schemes). See Users' Guide for options. DO NOT ADD ADDITIONAL COLUMNS. - num_land_cat = 21, ! Number of land categories from WPS input data. See Users' Guide for details. DO NOT ADD ADDITIONAL COLUMNS. - sf_urban_physics = 0, 0, 0, ! Option for activating urban canopy model (only for Noah LSM). See Users' Guide for details. Use same value for all domains. - / - - &fdda - / - - &dynamics - hybrid_opt = 2, ! Hybrid coordinate option (2=Klemp cubic form with etac; default beginning V4.0). DO NOT ADD ADDITIONAL COLUMNS. - w_damping = 0, ! Vertical velocity damping flag (0=off;1=on). DO NOT ADD ADDITIONAL COLUMNS. - diff_opt = 2, 2, 2, ! Turbulence and mixing option. See Users' Guide for options. - km_opt = 4, 4, 4, ! Eddy coefficient option. See Users' Guide for options. - diff_6th_opt = 0, 0, 0, ! 6th order numerical diffusion option. See Users' Guide for options. - diff_6th_factor = 0.12, 0.12, 0.12, ! 6th order numerical diffusion non-dimensional rate (max 1.0 = complete removal of 2dx wave in 1 timestep). - base_temp = 290. ! Base state temp (K). For real cases only. DO NOT ADD ADDITIONAL COLUMNS. - damp_opt = 3, ! Upper-level damping option. See Users' Guide for options. DO NOT ADD ADDITIONAL COLUMNS. - zdamp = 5000., 5000., 5000., ! Damping depth (m) from model top. - dampcoef = 0.2, 0.2, 0.2, ! Damping coefficient corresponding to damp_opt. See Users' Guide for details. - khdif = 0, 0, 0, ! Horizontal diffusion constant (m^2/s) - kvdif = 0, 0, 0, ! Vertical diffusion constant (m^2/s) - non_hydrostatic = .true., .true., .true., ! Set to false to run the model hydrostatically. - moist_adv_opt = 1, 1, 1, ! Advection option for moisture. See Users' Guide for details. - scalar_adv_opt = 1, 1, 1, ! Advection option for scalars. See Users' Guide for details. - gwd_opt = 1, 1, 0, ! Switch for gravity wave drag option (1=on;0=off) - / - - &bdy_control - spec_bdy_width = 5, ! Total number of rows for specified boundary value nudging for real cases only. DO NOT ADD ADDITIONAL COLUMNS. - specified = .true. ! Specified boundary condition for real cases only. DO NOT ADD ADDITIONAL COLUMNS. - / - - &grib2 - / - - &namelist_quilt - nio_tasks_per_group = 0, ! Number of I/O tasks per group for quilting. 0=no quilting. DO NOT ADD ADDITIONAL COLUMNS. - nio_groups = 1, ! Number of I/O groups for quilting. DO NOT ADD ADDITIONAL COLUMNS. - /