diff --git a/README b/README index 62766946e1..1308214e47 100644 --- a/README +++ b/README @@ -1,4 +1,4 @@ -WRF Model Version 4.4 +WRF Model Version 4.4.1 https://www2.mmm.ucar.edu/wrf/users/ diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 4c92f156b3..da9904a546 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -3145,6 +3145,7 @@ package gdscheme cu_physics==93 - - package sasscheme cu_physics==94 - - package osasscheme cu_physics==95 - - package kfscheme cu_physics==99 - state:w0avg +package risticscheme cu_physics==147 - - package g3tave cu_diag==1 - state:GD_CLOUD,GD_CLOUD2,GD_CLDFR,GD_CLOUD_A,GD_CLOUD2_A,kbcon_deep,ktop_deep,k22_deep diff --git a/chem/depend.chem b/chem/depend.chem index a7648798c7..01ea3f5121 100644 --- a/chem/depend.chem +++ b/chem/depend.chem @@ -229,7 +229,7 @@ module_mosaic_sect_intr.o: module_mosaic_coag1d.o module_mosaic_coag3d.o module_ module_mosaic_aerdynam_intr.o: module_mosaic_sect_intr.o module_mosaic_aerchem_intr.o -module_mosaic_addemiss.o: module_data_mosaic_asect.o module_data_sorgam.o +module_mosaic_addemiss.o: module_data_mosaic_asect.o module_data_sorgam.o module_gocart_dust.o module_dep_simple.o: module_data_sorgam.o module_aerosols_soa_vbs.o diff --git a/dyn_em/module_stoch.F b/dyn_em/module_stoch.F index 0f925924a9..8f1a6f13c2 100644 --- a/dyn_em/module_stoch.F +++ b/dyn_em/module_stoch.F @@ -108,7 +108,6 @@ SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & ipsy, ipey, jpsy, jpey, kpsy, kpey ) - USE module_dm USE module_configure USE module_domain, ONLY : domain USE module_wrf_error diff --git a/inc/version_decl b/inc/version_decl index 8cc4643cbc..25d2e88587 100644 --- a/inc/version_decl +++ b/inc/version_decl @@ -1 +1 @@ - CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.4' + CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.4.1' diff --git a/phys/Makefile b/phys/Makefile index 02d81793aa..f6c0b28bbf 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -80,6 +80,7 @@ MODULES = \ module_cu_scalesas.o \ module_cu_osas.o \ module_cu_kfcup.o \ + module_cu_ristic.o \ module_madwrf.o \ module_mp_radar.o \ module_mp_kessler.o \ diff --git a/phys/module_cu_ristic.F b/phys/module_cu_ristic.F new file mode 100644 index 0000000000..c7add716de --- /dev/null +++ b/phys/module_cu_ristic.F @@ -0,0 +1,583 @@ +!----------------------------------------------------------------------- +! +!WRF:MODEL_LAYER:PHYSICS +! +! To improve cloud and precipitation forecast we developed new convective +! scheme and we implemented it in WRF model. Convective clouds have always +! been a great challenge for meteorologists, among other things, due to +! the inability to describe processes of cloud formation, development and +! dissipation in a satisfactory manner. Applying parameterization in to +! the models has lead to simpler form of equations that could be used in +! practice, and thus different types of convective schemes in numerical +! weather prediction models appeared. Proposed convective scheme is based +! on basic elements that affect convection such as convective available +! potential energy (CAPE), vertical velocity at the base of the cloud, the +! amount of ice in the cloud and important assumptions. The scheme is +! conceived as a wet vertical turbulent diffusion and a logical +! continuation of dry vertical planetary boundary layer (PBL) turbulent +! diffusion. The scheme determines the vertical levels in the model where +! the convective cloud begins and ends. Integrated in the model this +! scheme showed good results in practice. +! A complete description is now found in Ristic I., Kordic I., April 2022: +! Convective velocity scale and its application in convective parametrization +! +! Author: Ivan Ristic, WEATHER2, ivanr@weather2.rs +! Last modified: 23 December 2021 +! +!----------------------------------------------------------------------- +! + MODULE MODULE_CU_RISTIC +! +!----------------------------------------------------------------------- +! + USE MODULE_MODEL_CONSTANTS +! +!----------------------------------------------------------------------- +! +CONTAINS +! +!----------------------------------------------------------------------- + SUBROUTINE RISTICDRV( & + & IDS,IDE,JDS,JDE,KDS,KDE & + & ,IMS,IME,JMS,JME,KMS,KME & + & ,ITS,ITE,JTS,JTE,KTS,KTE & + & ,USTAR & + & ,W,QC & + & ,DT,ITIMESTEP,STEPCU & + & ,RAINCV,PRATEC,CUTOP,CUBOT & + & ,TH,T,QV & + & ,PMID,PI,RHO,DZ8W & + & ,CP,R,ELWV,ELIV,G & + & ,CU_ACT_FLAG & + & ,RTHCUTEN,RQVCUTEN & + & ) +!----------------------------------------------------------------------- + IMPLICIT NONE +!----------------------------------------------------------------------- + INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE & + & ,IMS,IME,JMS,JME,KMS,KME & + & ,ITS,ITE,JTS,JTE,KTS,KTE +! + REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: USTAR + REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: W,QC +! + INTEGER,INTENT(IN) :: ITIMESTEP,STEPCU +! + REAL,INTENT(IN) :: CP,DT,ELIV,ELWV,G,R +! + REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ8W & + & ,PI & + & ,PMID,QV & + & ,RHO,T,TH +! + REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) & + & ,OPTIONAL & + & ,INTENT(INOUT) :: RQVCUTEN,RTHCUTEN +! + REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: & + PRATEC,RAINCV +! + REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CUBOT,CUTOP +! + LOGICAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CU_ACT_FLAG +! +!----------------------------------------------------------------------- +!*** +!*** LOCAL VARIABLES +!*** +!----------------------------------------------------------------------- + INTEGER :: LBOT,LTOP +! + REAL :: DTCNVC,PCPCOL +! + REAL,DIMENSION(KTS:KTE) :: DPCOL,DQDT,DTDT,PCOL,QCOL,TCOL + REAL,DIMENSION(KTS:KTE) :: WCOL,CWMCOL +! + INTEGER :: I,J,K,KFLIP + +!*** Begin debugging convection + REAL :: DELQ,DELT,PLYR + INTEGER :: IMD,JMD + LOGICAL :: PRINT_DIAG +!*** End debugging convection +! +!----------------------------------------------------------------------- +!*********************************************************************** +!----------------------------------------------------------------------- +! +!*** PREPARE TO CALL RISTIC CONVECTION SCHEME +! +!----------------------------------------------------------------------- + +!*** Begin debugging convection + IMD=(IMS+IME)/2 + JMD=(JMS+JME)/2 + PRINT_DIAG=.FALSE. +!*** End debugging convection + +! + DO J=JTS,JTE + DO I=ITS,ITE + CU_ACT_FLAG(I,J)=.TRUE. + ENDDO + ENDDO +! + DTCNVC=DT*STEPCU +! + DO J=JTS,JTE + DO I=ITS,ITE +! + DO K=KTS,KTE + DQDT(K)=0. + DTDT(K)=0. + ENDDO +! + PCPCOL=0. + RAINCV(I,J)=0. + PRATEC(I,J)=0. +! +!*** FILL 1-D VERTICAL ARRAYS +!*** AND FLIP DIRECTION SINCE RISTIC SCHEME +!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP +! + DO K=KTS,KTE + KFLIP=KTE+1-K +! +!*** CONVERT FROM MIXING RATIO TO SPECIFIC HUMIDITY +! + QCOL(K)=MAX(EPSQ,QV(I,KFLIP,J)/(1.+QV(I,KFLIP,J))) + TCOL(K)=T(I,KFLIP,J) + PCOL(K)=PMID(I,KFLIP,J) + DPCOL(K)=RHO(I,KFLIP,J)*G*DZ8W(I,KFLIP,J) + WCOL(K)=0.5*(W(I,KFLIP,J)+W(I,KFLIP+1,J))*DT/200. + CWMCOL(K)=QC(I,KFLIP,J) + ENDDO +!----------------------------------------------------------------------- +!*** +!*** CALL CONVECTION +!*** +!----------------------------------------------------------------------- +!*** Begin debugging convection +! PRINT_DIAG=.FALSE. +! IF(I==IMD.AND.J==JMD)PRINT_DIAG=.TRUE. +!*** End debugging convection +!----------------------------------------------------------------------- + CALL RISTIC(ITIMESTEP,I,J,DTCNVC & + & ,DPCOL,PCOL,QCOL,TCOL & + & ,USTAR(I,J) & + & ,WCOL,CWMCOL & + & ,DQDT,DTDT,PCPCOL,LBOT,LTOP & + & ,CP,R,ELWV,ELIV,G & + & ,PRINT_DIAG & + & ,IDS,IDE,JDS,JDE,KDS,KDE & + & ,IMS,IME,JMS,JME,KMS,KME & + & ,ITS,ITE,JTS,JTE,KTS,KTE) +!----------------------------------------------------------------------- +! +!*** COMPUTE HEATING AND MOISTENING TENDENCIES +! + IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN )) THEN + DO K=KTS,KTE + KFLIP=KTE+1-K + RTHCUTEN(I,K,J)=DTDT(KFLIP)/PI(I,K,J) +! +!*** CONVERT FROM SPECIFIC HUMIDTY BACK TO MIXING RATIO +! + RQVCUTEN(I,K,J)=DQDT(KFLIP)/(1.-QCOL(KFLIP))**2 + ENDDO + ENDIF +! +!*** ALL UNITS IN RISTIC SCHEME ARE MKS, THUS CONVERT PRECIP FROM METERS +!*** TO MILLIMETERS PER STEP FOR OUTPUT. +! + RAINCV(I,J)=PCPCOL*1.E3/STEPCU + PRATEC(I,J)=PCPCOL*1.E3/(STEPCU * DT) +! +!*** CONVECTIVE CLOUD TOP AND BOTTOM FROM THIS CALL +! + CUTOP(I,J)=REAL(KTE+1-LTOP) + CUBOT(I,J)=REAL(KTE+1-LBOT) +! +!----------------------------------------------------------------------- +! + ENDDO + ENDDO +! + END SUBROUTINE RISTICDRV +!----------------------------------------------------------------------- +!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +!----------------------------------------------------------------------- + SUBROUTINE RISTIC & +!----------------------------------------------------------------------- + & (ITIMESTEP,I,J,DTCNVC & + & ,DPRS,PRSMID,Q,T & + & ,USTAR & + & ,W,CWM & + & ,DQDT,DTDT,PCPCOL,LBOT,LTOP & + & ,CP,R,ELWV,ELIV,G & + & ,PRINT_DIAG & + & ,IDS,IDE,JDS,JDE,KDS,KDE & + & ,IMS,IME,JMS,JME,KMS,KME & + & ,ITS,ITE,JTS,JTE,KTS,KTE) +!----------------------------------------------------------------------- + IMPLICIT NONE +!----------------------------------------------------------------------- + INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE & + ,IMS,IME,JMS,JME,KMS,KME & + ,ITS,ITE,JTS,JTE,KTS,KTE & + ,I,J,ITIMESTEP +! + INTEGER,INTENT(OUT) :: LBOT,LTOP +! + REAL,INTENT(IN) :: CP,DTCNVC,ELIV,ELWV,G,R +! + REAL,DIMENSION(KTS:KTE),INTENT(IN) :: DPRS,PRSMID,Q,T + REAL,INTENT(IN) :: USTAR + REAL,DIMENSION(KTS:KTE),INTENT(IN) :: W,CWM +! + REAL,INTENT(INOUT) :: PCPCOL +! + REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: DQDT,DTDT +! +!----------------------------------------------------------------------- +!*** DEFINE LOCAL VARIABLES +!----------------------------------------------------------------------- +! + REAL,DIMENSION(KTS:KTE+1) :: TH,QH,APEL,TP,QP,FLUX,TN,QN,C2 +! + REAL,DIMENSION(KTS:KTE) :: CPCP,DZ,DIFT,DIFQ +! +!*** Begin debugging convection + LOGICAL :: PRINT_DIAG +!*** End debugging convection +! +!----------------------------------------------------------------------- +!*** +!*** LOCAL SCALARS +!*** + REAL :: CAPA,RCP,DTPH,RDTPH,TWODT,RTWODT & + & ,CAPE,WLCL,FVIJE,FCLOUD,FLUXB,CFAC,RHOL,TTHBT,QBT & + & ,ZLO,TKL,QKL,CWMKL,PKL,APE,TBT,WKL,THICK,RTHICK & + & ,RHO,ZINT,ZMID,DZKL & + & ,TT,QQ,PP,QC,QW,TC,TG1,QWAT,TVC,TVE,FLUNI,ZLCL,FRC1 & + & ,FPCP,C0,CPDR,ZTOP,FACW,FOUT,FFUP & + & ,RP,TGS,FO,TGUESS,F1,DTG,RROW,PRECRL,SUMQ,SUMT +! + INTEGER :: L,LMHK,LMHP,LMHM,LBTK,LTPK,ITCNT,NSTEP,ITER +! + REAL,PARAMETER :: ROW=1.E3 +!----------------------------------------------------------------------- + PCPCOL=0. + LBOT=0 + LTOP=KTE +!----------------------------------------------------------------------- +!-----PREPARATORY CALCULATIONS------------------------------------------ +!----------------------------------------------------------------------- + CAPA=R/CP + RCP=1./CP + RROW=1./ROW + DTPH=DTCNVC + RDTPH=1./DTPH + TWODT=DTPH + RTWODT=1./TWODT +!----------------------------------------------------------------------- +!-----START OF CONVECTION----------------------------------------------- +!----------------------------------------------------------------------- + LMHK=KTE + LMHP=LMHK+1 + LMHM=LMHK-1 +! + LBTK=LMHK + LTPK=LMHK +! + TP=0. + QP=0. + TH=0. + QH=0. + TN=0. + QN=0. + C2=0. + TTHBT=0. + QBT=0. + FLUX=0. + CPCP=0. + CAPE=0. + WLCL=0. + FVIJE=0. + FCLOUD=0. + ZLO=0. +! + FLUXB=1000. + C2(LMHP)=0.03*USTAR + RHOL=1. +! + DO 90 L=LMHK,1,-1 +! + TKL=T(L) + QKL=Q(L) + CWMKL=CWM(L) + PKL=PRSMID(L) + APE=(1.E5/PKL)**CAPA + TBT=TTHBT/APE + WKL=W(L) +! + THICK=DPRS(L)/G + RTHICK=1./THICK +! + RHO=PKL/(R*TKL*(1.+0.608*QKL)) + CFAC=RHO/RHOL + RHOL=RHO + ZINT=ZLO+THICK/RHO + DZKL=ZINT-ZLO + ZMID=0.5*(ZLO+ZINT) + ZLO=ZINT +! + DZ(L)=DZKL + TH(L)=(TKL-RCP*ELWV*CWMKL)*APE + QH(L)=QKL+CWMKL + APEL(L)=APE +!----------------------------------------------------------------------- +!-----QC---------------------------------------------------------------- +!----------------------------------------------------------------------- + TT=TBT + QQ=QBT + PP=PKL + CALL FSLOPE + QC=QW + TC=TG1 + QWAT=QBT-QC + IF(QWAT.LE.0.)THEN + QC=QBT + TC=TBT + QWAT=0. + ENDIF + TVC=TC*(1.+0.608*QC)/(1.+QWAT) + TVE=TKL*(1.+0.608*QKL)/(1.+CWMKL) + CAPE=CAPE+(TVC-TVE)/TVE*G*DZKL + IF(TVC.GT.TVE)FVIJE=1. + FLUNI=0.5+SIGN(0.5,T(L)-T(L-1)) +!----------------------------------------------------------------------- +!-----CLOUD BOTTOM------------------------------------------------------ +!----------------------------------------------------------------------- + IF(FCLOUD.EQ.0.)THEN + IF(2*L.LT.LMHK)THEN + GOTO 100 + ELSEIF(QWAT.LE.0..OR.CAPE.LE.0..AND.WLCL.LE.0.)THEN + FLUX=0. + FVIJE=0. + FLUXB=1000. + ZLCL=ZMID + LBTK=L + CAPE=C2(LMHP)*C2(LMHP)+4.64*AMAX1(0.,WKL)**(1./3.)/TVE*G*DZKL + WLCL=WKL + ELSE + IF(CAPE.LE.0.)GOTO 100 + FCLOUD=1. + ENDIF + ENDIF +!----------------------------------------------------------------------- +!-----PRECIPITATION----------------------------------------------------- +!----------------------------------------------------------------------- + IF(TC.GE.268.16)THEN + FRC1=0. + ELSEIF(TC.LE.248.16)THEN + FRC1=1. + ELSE + FRC1=(268.16-TC)/(268.16-248.16) + ENDIF + IF(FCLOUD*CAPE.LE.0.)THEN + IF(FCLOUD*0.35.GT.FLUNI*USTAR)THEN + FLUX(L+1)=0. + CPCP(L+1)=CPCP(L+1)/FPCP + ENDIF + FPCP=FCLOUD + ELSE + C0=0.03*(0.1+0.9*FRC1)/SQRT(CAPE) + FPCP=1.-EXP(-C0*DZKL) + ENDIF + CPDR=QWAT*FPCP + QBT=QBT-CPDR + TTHBT=TTHBT+RCP*ELWV*CPDR*APE + CPCP(L)=FLUX(L+1)*CPDR*FRC1 +!----------------------------------------------------------------------- +!-----CLOUD TOP--------------------------------------------------------- +!----------------------------------------------------------------------- + IF(FCLOUD.EQ.1.)THEN + IF(CAPE.LE.0..OR.QWAT.LE.0.)THEN + FLUX(L)=0. + LTPK=L + ZTOP=ZMID + GOTO 95 + ENDIF + ELSEIF(CAPE.LE.0.)THEN + GOTO 100 + ENDIF +!----------------------------------------------------------------------- +!-----FLUX-------------------------------------------------------------- +!----------------------------------------------------------------------- + C2(L)=(C2(LMHP)*CAPE)**(1./3.) + FLUX(L)=RHO*C2(L)*TWODT + FLUXB=AMIN1(FLUXB,DZKL/SQRT(CAPE)*RTWODT) +!----------------------------------------------------------------------- +!-----UPDRAFT ENTRAINMENT----------------------------------------------- +!----------------------------------------------------------------------- + FACW=(FLUX(L)-FLUX(L+1))/FLUX(L) + FOUT=AMIN1(1.,0.5*(1.-FACW)*(3.*CFAC-1.)) + FFUP=1.-FOUT + TTHBT=TTHBT*FOUT+TH(L)*FFUP + QBT=QBT*FOUT+QH(L)*FFUP + TP(L)=TTHBT + QP(L)=QBT +!----------------------------------------------------------------------- + 90 CONTINUE +!----------------------------------------------------------------------- + 95 CONTINUE +!----------------------------------------------------------------------- +!-----NO CONVECTION----------------------------------------------------- +!----------------------------------------------------------------------- + CPCP=CPCP*FVIJE +!----------------------------------------------------------------------- +!-----PRECIPITATION----------------------------------------------------- +!----------------------------------------------------------------------- + DIFT=0. + DIFQ=0. + PRECRL=0. + NSTEP=NINT(1./FLUXB+0.5) + FLUXB=1./FLOAT(NSTEP) + DO ITER=1,NSTEP + DO L=LTPK,LMHK + THICK=DPRS(L)/G + RTHICK=1./THICK + TN(L)=TH(L)+FLUXB*RTHICK*(RCP*ELWV*CPCP(L)*APEL(L) & + & +FLUX(L+1)*(TP(L+1)-TH(L )) & + & -FLUX(L )*(TP(L )-TH(L-1))) + QN(L)=QH(L)+FLUXB*RTHICK*(-CPCP(L) & + & +FLUX(L+1)*(QP(L+1)-QH(L )) & + & -FLUX(L )*(QP(L )-QH(L-1))) + PRECRL=PRECRL+FLUXB*CPCP(L) + ENDDO + TH=TN + QH=QN + ENDDO +!----------------------------------------------------------------------- +!-------------------THE PRECIPITATION ON SFC---------------------------- +!----------------------------------------------------------------------- + SUMQ=0. + SUMT=0. + DO L=LTPK,LMHK + DIFT(L)=TH(L)/APEL(L)+CWM(L)*RCP*ELWV-T(L) + DIFQ(L)=QH(L)-CWM(L)-Q(L) + THICK=DPRS(L)/G + SUMQ=SUMQ+DIFQ(L)*THICK + SUMT=SUMT+DIFT(L)*THICK*APEL(L)*CP/ELWV + ENDDO + IF(PRINT_DIAG)write(6,*) PRECRL,SUMT,SUMQ,FLUXB +! +!--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------ +! + LTOP=LTPK + LBOT=LBTK + PCPCOL=PRECRL*RROW + DQDT=DIFQ*RTWODT + DTDT=DIFT*RTWODT +!----------------------------------------------------------------------- + 100 CONTINUE + RETURN + CONTAINS + SUBROUTINE FSLOPE + IF(L.EQ.LMHK)THEN + TG1=TT + QW=QQ + GOTO 25 + ENDIF + RP=TT+RCP*ELWV*QQ + TGS=TT + QW=QSAT(PP,TGS,0.,0.) + FO=TGS+RCP*ELWV*QW-RP + TG1=AMAX1(TGS/2.,TGS-.5*FO) + TGUESS=TGS + ITCNT=0 + 10 QW=QSAT(PP,TG1,0.,0.) + F1=TG1+RCP*ELWV*QW-RP + IF(ABS(F1).LT..001.OR.ABS(F1-FO).LT.1.E-10.OR.ITCNT.GT.30)GOTO 25 + ITCNT=ITCNT+1 + DTG=F1*(TG1-TGUESS)/(F1-FO) + TGUESS=TG1 + FO=F1 + TG1=TG1-DTG + GOTO 10 + 25 CONTINUE + END SUBROUTINE FSLOPE +!----------------------------------------------------------------------- + END SUBROUTINE RISTIC +!----------------------------------------------------------------------- + REAL FUNCTION QSAT(PP,TT,FIW,FCI) +!----------------------------------------------------------------------- + REAL, PARAMETER :: PQ0=379.90516,A2=17.2693882,A3=273.16,A4=35.86 +!----------------------------------------------------------------------- +!-----AI, BI------------------------------------------------------------ +!----------------------------------------------------------------------- + TMT0=TT-A3 + + IF(TMT0.LT.-20.)THEN + AI=0.007225 + BI=0.9674 + ELSEIF(TMT0.LT.0.)THEN + AI=0.008855 + BI=1. + ELSE + AI=0. + BI=1. + IF(FIW.EQ.1.)TMT0=0. + ENDIF + + AI=FIW*AI*FCI + BI=1.-FIW+FIW*BI*FCI + + QSAT=PQ0/PP*EXP(A2*TMT0/(TT-A4))*(BI+AI*TMT0) +!----------------------------------------------------------------------- + RETURN + END +!----------------------------------------------------------------------- + SUBROUTINE RISTICINIT(RTHCUTEN,RQVCUTEN & + & ,RESTART & + & ,ALLOWED_TO_READ & + & ,IDS,IDE,JDS,JDE,KDS,KDE & + & ,IMS,IME,JMS,JME,KMS,KME & + & ,ITS,ITE,JTS,JTE,KTS,KTE) +!----------------------------------------------------------------------- + IMPLICIT NONE +!----------------------------------------------------------------------- + LOGICAL,INTENT(IN) :: RESTART,ALLOWED_TO_READ + INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE & + & ,IMS,IME,JMS,JME,KMS,KME & + & ,ITS,ITE,JTS,JTE,KTS,KTE +! + REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: & + & RTHCUTEN & + & ,RQVCUTEN +! + INTEGER :: I,J,K,ITF,JTF,KTF +!----------------------------------------------------------------------- + + JTF=MIN0(JTE,JDE-1) + KTF=MIN0(KTE,KDE-1) + ITF=MIN0(ITE,IDE-1) +! + IF(.NOT.RESTART)THEN + DO J=JTS,JTF + DO K=KTS,KTF + DO I=ITS,ITF + RTHCUTEN(I,K,J)=0. + RQVCUTEN(I,K,J)=0. + ENDDO + ENDDO + ENDDO + ENDIF +!----------------------------------------------------------------------- + END SUBROUTINE RISTICINIT +!----------------------------------------------------------------------- +! + END MODULE MODULE_CU_RISTIC +! +!----------------------------------------------------------------------- diff --git a/phys/module_cumulus_driver.F b/phys/module_cumulus_driver.F index c6d3b1f361..9a526eceec 100644 --- a/phys/module_cumulus_driver.F +++ b/phys/module_cumulus_driver.F @@ -150,6 +150,7 @@ SUBROUTINE cumulus_driver(grid & ,OSASSCHEME & ,SCALESASSCHEME & ! scale-sware sas ,KSASSCHEME, NSASSCHEME & + ,RISTICSCHEME & #if (EM_CORE == 1) ,MSKFSCHEME & ,CAMMGMPSCHEME & @@ -169,11 +170,11 @@ SUBROUTINE cumulus_driver(grid & USE module_state_description, ONLY: DUCUSCHEME #endif - ! *** add new modules of schemes here USE module_cu_kf , ONLY : kfcps USE module_cu_bmj , ONLY : bmjdrv + USE module_cu_ristic , ONLY : risticdrv #ifdef DM_PARALLEL USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks #if (EM_CORE == 1) @@ -402,7 +403,6 @@ SUBROUTINE cumulus_driver(grid & !BSINGH - ENDS #endif - INTEGER, INTENT(IN ) :: cu_physics INTEGER, INTENT(IN ) :: STEPCU LOGICAL, INTENT(IN ) :: warm_rain @@ -1008,6 +1008,25 @@ SUBROUTINE cumulus_driver(grid & ,RTHCUTEN=rthcuten ,RQVCUTEN=rqvcuten & ) + CASE (RISTICSCHEME) + CALL wrf_debug(100,'in ristic_cps') + CALL RISTICDRV( & + TH=th,T=T ,RAINCV=raincv, PRATEC=tmppratec & + ,RHO=rho & + ,DT=dt ,ITIMESTEP=itimestep ,STEPCU=stepcu & + ,CUTOP=htop, CUBOT=hbot & + ,DZ8W=dz8w, PMID=p, PI=pi & + ,CP=cp ,R=r_d ,ELWV=xlv ,ELIV=xls ,G=g & + ,CU_ACT_FLAG=cu_act_flag & + ,QV=qv_curr & + ,ustar=ust & + ,W=w,QC=qc_curr & + ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde & + ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme & + ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte & + ,RTHCUTEN=rthcuten ,RQVCUTEN=rqvcuten & + ) + CASE (KFETASCHEME) CALL wrf_debug(100,'in kf_eta_cps') CALL KF_ETA_CPS( & diff --git a/phys/module_physics_addtendc.F b/phys/module_physics_addtendc.F index 06136a1226..eaa1f5ef0a 100644 --- a/phys/module_physics_addtendc.F +++ b/phys/module_physics_addtendc.F @@ -1124,7 +1124,7 @@ SUBROUTINE phy_cu_ten(config_flags,rk_step,n_moist,n_scalar, & ENDIF - CASE (BMJSCHEME) + CASE (BMJSCHEME, RISTICSCHEME) CALL add_a2a(rt_tendf,RTHCUTEN, & config_flags, & ids,ide, jds, jde, kds, kde, & @@ -2131,7 +2131,7 @@ SUBROUTINE advance_ppt(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, & ENDDO ENDDO - CASE (BMJSCHEME, CAMZMSCHEME) + CASE (BMJSCHEME, CAMZMSCHEME, RISTICSCHEME) DO J = j_start,j_end DO i = i_start,i_end diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 8725908b3a..1c76f98407 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -3975,6 +3975,7 @@ SUBROUTINE cu_init(DX,STEPCU,CUDT,DT,RUCUTEN,RVCUTEN,RTHCUTEN, & USE module_cu_kfeta USE module_cu_mskf USE MODULE_CU_BMJ + USE MODULE_CU_RISTIC USE module_cu_gd, ONLY : GDINIT USE module_cu_g3, ONLY : G3INIT USE module_cu_sas @@ -4087,6 +4088,14 @@ SUBROUTINE cu_init(DX,STEPCU,CUDT,DT,RUCUTEN,RVCUTEN,RTHCUTEN, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) + CASE (RISTICSCHEME) + CALL risticinit(RTHCUTEN,RQVCUTEN, & + restart, & + allowed_to_read, & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + its, ite, jts, jte, kts, kte ) + CASE (KFETASCHEME) CALL kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, & RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, & diff --git a/phys/module_ra_clWRF_support.F b/phys/module_ra_clWRF_support.F index dd0fac51a0..9eb3d1b78f 100755 --- a/phys/module_ra_clWRF_support.F +++ b/phys/module_ra_clWRF_support.F @@ -168,7 +168,13 @@ SUBROUTINE read_CAMgases(yr, julian, READtrFILE, model, co2vmr, n2ovmr, ch4vmr, CALL wrf_message( message) ENDIF - IF (istatus /= -1) THEN + IF (istatus == -1) THEN + WRITE(message,*) "Normal ending of CAMtr_volume_mixing_ratio file" + CALL wrf_message( message) + ELSE IF (istatus == -4001) THEN ! Cray CCE compiler throws -4001 rather than -1 + WRITE(message,*) "Normal ending of CAMtr_volume_mixing_ratio file" + CALL wrf_message( message) + ELSE ! if not -1 or -4001, then abort WRITE(message,*) " Not normal ending of CAMtr_volume_mixing_ratio file" call wrf_error_fatal( message) END IF diff --git a/phys/module_ra_rrtmg_lw.F b/phys/module_ra_rrtmg_lw.F index b75f36d3af..e31d8c3ab4 100644 --- a/phys/module_ra_rrtmg_lw.F +++ b/phys/module_ra_rrtmg_lw.F @@ -11701,9 +11701,8 @@ SUBROUTINE RRTMG_LWRAD( & ! Ozone REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & - OPTIONAL , & - INTENT(IN ) :: O33D - INTEGER, OPTIONAL, INTENT(IN ) :: o3input + INTENT(INOUT) :: O33D + INTEGER, INTENT(IN ) :: o3input real, parameter :: thresh=1.e-9 real slope @@ -12011,7 +12010,7 @@ SUBROUTINE RRTMG_LWRAD( & QV1D(K)=max(0.,QV1D(K)) ENDDO - IF (PRESENT(O33D)) THEN + IF (o3input.eq.2) THEN DO K=kts,kte O31D(K)=O33D(I,K,J) ENDDO @@ -12401,27 +12400,31 @@ SUBROUTINE RRTMG_LWRAD( & call inirad (o3mmr,plev,kts,nlay-1) ! Steven Cavallo: Changed to nlayers from kte+1 - if(present(o33d)) then + if(o3input.eq.2) then do k = kts, nlayers o3vmr(ncol,k) = o3mmr(k) * amdo - IF ( PRESENT( O33D ) ) THEN - if(o3input .eq. 2)then - if(k.le.kte)then - o3vmr(ncol,k) = o31d(k) - else + if(k.le.kte)then + o3vmr(ncol,k) = o31d(k) + else ! apply shifted climatology profile above model top - o3vmr(ncol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo - if(o3vmr(ncol,k) .le. 0.)o3vmr(ncol,k) = o3mmr(k)*amdo - endif + o3vmr(ncol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo + if(o3vmr(ncol,k) .le. 0.)o3vmr(ncol,k) = o3mmr(k)*amdo endif - ENDIF enddo else do k = kts, nlayers o3vmr(ncol,k) = o3mmr(k) * amdo + o31d(k) = o3vmr(ncol,k) enddo endif +! output o3 for o3input=0 + IF (o3input.ne.2) THEN + DO K=kts,kte + O33D(I,K,J)=O31D(K) + ENDDO + ENDIF + ! Set surface emissivity in each RRTMG longwave band do nb = 1, nbndlw emis(ncol, nb) = emiss(i,j) diff --git a/phys/module_ra_rrtmg_lwf.F b/phys/module_ra_rrtmg_lwf.F index ae4b1dea7c..c691bec1f0 100644 --- a/phys/module_ra_rrtmg_lwf.F +++ b/phys/module_ra_rrtmg_lwf.F @@ -15375,8 +15375,8 @@ SUBROUTINE RRTMG_LWRAD_FAST( & ! Ozone REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & OPTIONAL , & - INTENT(IN ) :: O33D - INTEGER, OPTIONAL, INTENT(IN ) :: o3input + INTENT(INOUT) :: O33D + INTEGER, INTENT(IN ) :: o3input real, parameter :: thresh=1.e-9 real slope @@ -15674,7 +15674,7 @@ SUBROUTINE RRTMG_LWRAD_FAST( & QV1D(K)=max(0.,QV1D(K)) ENDDO - IF (PRESENT(O33D)) THEN + IF (o3input.eq.2) THEN DO K=kts,kte O31D(K)=O33D(I,K,J) ENDDO @@ -16015,20 +16015,16 @@ SUBROUTINE RRTMG_LWRAD_FAST( & call inirad (o3mmr,plev(icol,:),kts,nlay-1) ! Steven Cavallo: Changed to nlayers from kte+1 - if(present(o33d)) then + if(o3input.eq.2) then do k = kts, nlayers o3vmr(icol,k) = o3mmr(k) * amdo - IF ( PRESENT( O33D ) ) THEN - if(o3input .eq. 2)then - if(k.le.kte)then - o3vmr(icol,k) = o31d(k) - else + if(k.le.kte)then + o3vmr(icol,k) = o31d(k) + else ! apply shifted climatology profile above model top - o3vmr(icol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo - if(o3vmr(icol,k) .le. 0.)o3vmr(icol,k) = o3mmr(k)*amdo - endif + o3vmr(icol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo + if(o3vmr(icol,k) .le. 0.)o3vmr(icol,k) = o3mmr(k)*amdo endif - ENDIF enddo else do k = kts, nlayers diff --git a/phys/module_ra_rrtmg_sw.F b/phys/module_ra_rrtmg_sw.F index 1b0c2a8b88..c0eb328a4d 100644 --- a/phys/module_ra_rrtmg_sw.F +++ b/phys/module_ra_rrtmg_sw.F @@ -10214,9 +10214,8 @@ SUBROUTINE RRTMG_SWRAD( & ! Ozone REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & - OPTIONAL , & - INTENT(IN ) :: O33D - INTEGER, OPTIONAL, INTENT(IN ) :: o3input + INTENT(INOUT) :: O33D + INTEGER, INTENT(IN ) :: o3input ! EC aerosol: no_src = naerec = 6 INTEGER, INTENT(IN ) :: no_src REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:no_src ) , & @@ -10598,7 +10597,7 @@ SUBROUTINE RRTMG_SWRAD( & QV1D(K)=max(0.,QV1D(K)) ENDDO - IF (PRESENT(O33D)) THEN + IF (o3input.eq.2) THEN DO K=kts,kte O31D(K)=O33D(I,K,J) ENDDO @@ -10919,24 +10918,21 @@ SUBROUTINE RRTMG_SWRAD( & ! Get ozone profile including amount in extra layer above model top call inirad (o3mmr,plev,kts,kte) - if(present(o33d)) then + if(o3input.eq.2) then do k = kts, kte+1 o3vmr(ncol,k) = o3mmr(k) * amdo - IF ( PRESENT( O33D ) ) THEN - if(o3input .eq. 2)then - if(k.le.kte)then - o3vmr(ncol,k) = o31d(k) - else + if(k.le.kte)then + o3vmr(ncol,k) = o31d(k) + else ! apply shifted climatology profile above model top - o3vmr(ncol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo - if(o3vmr(ncol,k) .le. 0.)o3vmr(ncol,k) = o3mmr(k)*amdo - endif + o3vmr(ncol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo + if(o3vmr(ncol,k) .le. 0.)o3vmr(ncol,k) = o3mmr(k)*amdo endif - ENDIF enddo else do k = kts, kte+1 o3vmr(ncol,k) = o3mmr(k) * amdo + o31d(k) = o3vmr(ncol,k) enddo endif diff --git a/phys/module_ra_rrtmg_swf.F b/phys/module_ra_rrtmg_swf.F index ad4f454a7c..9ba21d34ec 100644 --- a/phys/module_ra_rrtmg_swf.F +++ b/phys/module_ra_rrtmg_swf.F @@ -11386,8 +11386,7 @@ SUBROUTINE RRTMG_SWRAD_FAST( & INTEGER, INTENT(IN ), OPTIONAL :: progn ! Ozone REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & - OPTIONAL , & - INTENT(IN ) :: O33D + INTENT(INOUT) :: O33D INTEGER, OPTIONAL, INTENT(IN ) :: o3input ! EC aerosol: no_src = naerec = 6 INTEGER, INTENT(IN ) :: no_src @@ -11731,7 +11730,7 @@ SUBROUTINE RRTMG_SWRAD_FAST( & QV1D(K)=max(0.,QV1D(K)) ENDDO - IF (PRESENT(O33D)) THEN + IF (o3input.eq.2) THEN DO K=kts,kte O31D(K)=O33D(I,K,J) ENDDO @@ -11992,20 +11991,16 @@ SUBROUTINE RRTMG_SWRAD_FAST( & ! call inirad (o3mmr,plev,kts,kte) call inirad (o3mmr,plev(icol,:),kts,kte) - if(present(o33d)) then + if(o3input.eq.2) then do k = kts, kte+1 o3vmr(icol,k) = o3mmr(k) * amdo - IF ( PRESENT( O33D ) ) THEN - if(o3input .eq. 2)then - if(k.le.kte)then - o3vmr(icol,k) = o31d(k) - else + if(k.le.kte)then + o3vmr(icol,k) = o31d(k) + else ! apply shifted climatology profile above model top - o3vmr(icol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo - if(o3vmr(icol,k) .le. 0.)o3vmr(icol,k) = o3mmr(k)*amdo - endif + o3vmr(icol,k) = o31d(kte) - o3mmr(kte)*amdo + o3mmr(k)*amdo + if(o3vmr(icol,k) .le. 0.)o3vmr(icol,k) = o3mmr(k)*amdo endif - ENDIF enddo else do k = kts, kte+1 diff --git a/run/README.namelist b/run/README.namelist index 1609978d9b..082e530135 100644 --- a/run/README.namelist +++ b/run/README.namelist @@ -406,13 +406,13 @@ Namelist variables for controlling the adaptive time step option: max(vert cfl, horiz cfl) <= target_cfl, then the time will increase by max_step_increase_pct. Use something large for nests (51% suggested) - starting_time_step(max_dom) = -1,-1 ! flag = -1 implies use 6 * dx (defined in start_em), + starting_time_step(max_dom) = -1,-1 ! flag = -1 implies use 4*dx (defined in start_em), starting_time_step = 100 means the starting time step for the coarse grid is 100 s - max_time_step(max_dom) = -1,-1 ! flag = -1 implies max time step is 3 * starting_time_step, + max_time_step(max_dom) = -1,-1 ! flag = -1 implies max time step is 8*dx, max_time_step = 100 means that the time step will not exceed 100 s - min_time_step(max_dom) = -1,-1 ! flag = -1 implies max time step is 0.5 * starting_time_step, + min_time_step(max_dom) = -1,-1 ! flag = -1 implies max time step is 3*dx, min_time_step = 100 means that the time step will not be less than 100 s adaptation_domain = 1 ! default, all fine grid domains adaptive dt driven by coarse-grid diff --git a/run/VEGPARM.TBL b/run/VEGPARM.TBL index 930f5d0267..774bd502d2 100644 --- a/run/VEGPARM.TBL +++ b/run/VEGPARM.TBL @@ -181,11 +181,27 @@ NATURAL CROP 12 LCZ_1 -24 +31 LCZ_2 -26 +32 LCZ_3 -99 +33 +LCZ_4 +34 +LCZ_5 +35 +LCZ_6 +36 +LCZ_7 +37 +LCZ_8 +38 +LCZ_9 +39 +LCZ_10 +40 +LCZ_11 +41 Vegetation Parameters USGS-RUC 28,1, 'ALBEDO Z0 LEMI PC SHDFAC IFOR RS RGL HS SNUP LAI MAXALB' diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index 0d2b367b92..c6da058d6f 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -2507,11 +2507,11 @@ END FUNCTION bep_bem_ngr_u END DO !----------------------------------------------------------------------- -! Check if qna_update=0 when aer_init_opt>0 for Thompson-MP-Aero (mp_physics=28) +! Check if qna_update=0 when aer_init_opt>1 for Thompson-MP-Aero (mp_physics=28) !----------------------------------------------------------------------- DO i = 1, model_config_rec % max_dom IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO ) THEN - IF ( model_config_rec%aer_init_opt .GT. 0 .and. model_config_rec%qna_update .EQ. 0 ) THEN + IF ( model_config_rec%aer_init_opt .GT. 1 .and. model_config_rec%qna_update .EQ. 0 ) THEN wrf_err_message = '--- ERROR: Time-varying sfc aerosol emissions not set for mp_physics=28 ' CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) wrf_err_message = '--- ERROR: Please set qna_update=1 and control through auxinput17 options ' @@ -2555,6 +2555,21 @@ END FUNCTION bep_bem_ngr_u END IF END DO +!----------------------------------------------------------------------- +! Stop the model if full_khain_lynn or mp_physics = 32 is selected +!----------------------------------------------------------------------- + DO i = 1, model_config_rec % max_dom + IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE + IF ( model_config_rec%mp_physics(i) .eq. full_khain_lynn) THEN + oops = oops + 1 + wrf_err_message = '--- ERROR: full bin spectral microphysics should not be used ' + CALL wrf_message ( wrf_err_message ) + wrf_err_message = '--- ERROR: use fast version instead (mp_physics=30)' + CALL wrf_message ( wrf_err_message ) + count_fatal_error = count_fatal_error + 1 + END IF + ENDDO ! Loop over domains + !----------------------------------------------------------------------- ! DJW Check that we're not using ndown and vertical nesting. !----------------------------------------------------------------------- diff --git a/tools/data.h b/tools/data.h index 8841372e75..081ece8616 100644 --- a/tools/data.h +++ b/tools/data.h @@ -75,7 +75,7 @@ typedef struct node_struct { char pkg_4dscalars[NAMELEN_LONG] ; /* fields used by Comm (halo, period, xpose) nodes */ - char comm_define[2*8192] ; + char comm_define[4*8192] ; /* marker */ int mark ;