Merge tag 'v4.2.2' into letmaik/wrf-4.2.2

Finalize WRFV4.2.2 by merging bug fixes from release-v4.2.2 branch onto master.
See the annotation of the v4.2.2 tag for full release notes (git show v4.2.2)

README

@@ -1,4 +1,4 @@
-WRF Model Version 4.2.1
+WRF Model Version 4.2.2
 
 http://www2.mmm.ucar.edu/wrf/users/
 

Registry/Registry.EM_COMMON

@@ -3063,7 +3063,7 @@ package   icedepth_one   seaice_thickness_opt==1     -             state:icedept
 #Time series options for text output
 package   notseries             process_time_series==0                 -             -
 package   tseries               process_time_series==1                 -             state:ts_hour,ts_u,ts_v,ts_q,ts_t,ts_psfc,ts_glw,ts_gsw,ts_hfx,ts_lh,ts_tsk,ts_tslb,ts_clw,ts_rainc,ts_rainnc,ts_u_profile,ts_v_profile,ts_gph_profile,ts_th_profile,ts_p_profile,ts_w_profile
-package   tseries_add_solar     process_time_series==2                 -             state:ts_hour,ts_u,ts_v,ts_q,ts_t,ts_psfc,ts_glw,ts_gsw,ts_hfx,ts_lh,ts_tsk,ts_tslb,ts_clw,ts_rainc,ts_rainnc,ts_u_profile,ts_v_profile,ts_gph_profile,ts_th_profile,ts_cldfrac2d,ts_wvp,ts_lwp,ts_iwp,ts_swp,ts_lwp_tot,ts_iwp_tot,ts_swp_tot,ts_re_qc,ts_re_qi,ts_re_qs,ts_re_qc_tot,ts_re_qi_tot,ts_re_qs_tot,ts_tau_qc,ts_tau_qi,ts_tau_qs,ts_tau_qc_tot,ts_tau_qi_tot,ts_tau_qs_tot,ts_cbaseht,ts_ctopht,ts_cbaseht_tot,ts_ctopht_tot,ts_clrnidx,ts_p_profile,ts_w_profile
+package   tseries_add_solar     process_time_series==2                 -             state:ts_hour,ts_u,ts_v,ts_q,ts_t,ts_psfc,ts_glw,ts_gsw,ts_hfx,ts_lh,ts_tsk,ts_tslb,ts_clw,ts_rainc,ts_rainnc,ts_u_profile,ts_v_profile,ts_gph_profile,ts_th_profile,ts_cldfrac2d,ts_wvp,ts_lwp,ts_iwp,ts_swp,ts_lwp_tot,ts_iwp_tot,ts_swp_tot,ts_re_qc,ts_re_qi,ts_re_qs,ts_re_qc_tot,ts_re_qi_tot,ts_re_qs_tot,ts_tau_qc,ts_tau_qi,ts_tau_qs,ts_tau_qc_tot,ts_tau_qi_tot,ts_tau_qs_tot,ts_cbaseht,ts_ctopht,ts_cbaseht_tot,ts_ctopht_tot,ts_clrnidx,ts_p_profile,ts_w_profile,ts_swdown,ts_swddni,ts_swddif,ts_swdownc,ts_swddnic,ts_swdown2,ts_swddni2,ts_swddif2,ts_swdownc2,ts_swddnic2
 
 # WRF-HAILCAST
 state    real   HAILCAST_DHAIL1   ij     misc        1         -      r          "HAILCAST_DHAIL1"     "WRF-HAILCAST Hail Diameter, 1st rank order"     "mm"

Registry/registry.solar_fields

@@ -61,6 +61,16 @@ state    real   ts_cbaseht_tot  ?!       misc       -         -     -      "TS_C
 state    real   ts_ctopht_tot   ?!       misc       -         -     -      "TS_CTOPHT_TOT"   "CLOUD TOP HEIGHT RES + UNRES"
 state    real   ts_clrnidx      ?!       misc       -         -     -      "TS_CLRNIDX"      "CLEARNESS INDEX"
 state    real   ts_sza          ?!       misc       -         -     -      "TS_SZA"          "SOLAR ZENITH ANGLE"
+state    real   ts_swdown       ?!       misc       -         -     -      "TS_SWDOWN"       "DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE"
+state    real   ts_swddni       ?!       misc       -         -     -      "TS_SWDDNI"       "SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE"
+state    real   ts_swddif       ?!       misc       -         -     -      "TS_SWDDIF"       "SHORTWAVE SURFACE DOWNWARD DIFFUSE IRRADIANCE"
+state    real   ts_swdownc      ?!       misc       -         -     -      "TS_SWDOWNC"      "DOWNWARD CLEAR-SKY SHORTWAVE FLUX AT GROUND SURFACE"
+state    real   ts_swddnic      ?!       misc       -         -     -      "TS_SWDDNIC"      "CLEAR-SKY SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE"
+state    real   ts_swdown2      ?!       misc       -         -     -      "TS_SWDOWN2"      "DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE FROM FARMS"
+state    real   ts_swddni2      ?!       misc       -         -     -      "TS_SWDDNI2"      "SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE FROM FARMS"
+state    real   ts_swddif2      ?!       misc       -         -     -      "TS_SWDDIF2"      "SHORTWAVE SURFACE DOWNWARD DIFFUSE IRRADIANCE FROM FARMS"
+state    real   ts_swdownc2     ?!       misc       -         -     -      "TS_SWDOWNC2"     "DOWNWARD CLEAR-SKY SHORTWAVE FLUX AT GROUND SURFACE FROM FARMS"
+state    real   ts_swddnic2     ?!       misc       -         -     -      "TS_SWDDNIC2"     "CLEAR-SKY SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE FROM FARMS"
 
 #  Package declarations
 

arch/configure.defaults

@@ -44,7 +44,7 @@ RLFLAGS		=
 CC_TOOLS        =      cc 
 
 ###########################################################
-#ARCH    Linux i486 i586 i686, gfortran compiler with gcc #serial smpar dmpar dm+sm
+#ARCH    Linux i486 i586 i686 armv7l aarch64, gfortran compiler with gcc #serial smpar dmpar dm+sm
 #
 DESCRIPTION     =       GNU ($SFC/$SCC)
 DMPARALLEL      =       # 1
@@ -1834,7 +1834,7 @@ RLFLAGS		=
 CC_TOOLS        =      $(SCC) 
 
 ###########################################################
-#ARCH    Linux x86_64 ppc64le i486 i586 i686 #serial smpar dmpar dm+sm
+#ARCH    Linux x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm
 #
 DESCRIPTION     =       INTEL ($SFC/$SCC): HSW/BDW
 DMPARALLEL      =       # 1
@@ -1878,7 +1878,7 @@ RLFLAGS		=
 CC_TOOLS        =      $(SCC) 
 
 ###########################################################
-#ARCH    Linux KNL x86_64 ppc64le i486 i586 i686 #serial smpar dmpar dm+sm
+#ARCH    Linux KNL x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm
 #
 DESCRIPTION     =       INTEL ($SFC/$SCC): KNL MIC
 DMPARALLEL      =       # 1

chem/chemics_init.F

@@ -1944,7 +1944,11 @@ subroutine chem_init (id,chem,emis_ant,scalar,dt,bioemdt,photdt,chemdt,stepbioe,
 !
    drydep_select: SELECT CASE(config_flags%gas_drydep_opt)
      CASE (WESELY)
-       CALL wrf_debug(15,'initializing dry dep (wesely)')
+        IF (numgas .eq. 0) THEN
+            CALL wrf_error_fatal("ERROR: numgas = 0, SELECTED CHEM OPT IS &
+            &NOT COMPATIBLE WITH WESELY DRY DEPOSITION")
+        ENDIF
+        CALL wrf_debug(15,'initializing dry dep (wesely)') 
 
         call dep_init( id, config_flags, numgas, mminlu_loc, &
                       its, ite, jts, jte, ide, jde )

chem/module_emissions_anthropogenics.F

@@ -241,7 +241,7 @@ subroutine add_anthropogenics(id,dtstep,dz8w,config_flags,rho_phy,alt,  &
    if (config_flags%aircraft_emiss_opt == 1 ) then
       do j=jts,jte  
         do k=kts,min(config_flags%kemit_aircraft,kte)
-          conv_rho(its:ite)=4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
+          conv_rho(its:ite)=4.828e-4/rho_phy(its:ite,k,j)*dtstep/(dz8w(its:ite,k,j)*60.)
           if( p_no >= param_first_scalar ) then
             chem(its:ite,k,j,p_no)  = chem(its:ite,k,j,p_no)  + emis_aircraft(its:ite,k,j,p_eac_no) *conv_rho(its:ite)
           endif

chem/module_mosaic_addemiss.F

@@ -601,7 +601,7 @@ subroutine mosaic_addemiss( id, dtstep, u10, v10, alt, dz8w, xland,     &
               do i = its, ite
 ! compute mass biomass burning emissions [(ug/m3)*m/s] for each species
 ! using the apportioning fractions
-                aem_so4 = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_sulf)      
+                IF ( p_ebu_sulf .gt. 1) aem_so4 = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_sulf)      
                 aem_oc  = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_oc)    
                 aem_bc  = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_bc) 
 ! Option to calculate OIN fraction of total PM for fire emissions 

dyn_em/module_big_step_utilities_em.F

@@ -1506,7 +1506,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
    jtf=MIN(jte,jde-1)
 
    IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
-   IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
+   IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jte-2
 
    DO j = j_start, jtf
 
@@ -1539,7 +1539,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
    jtf=MIN(jte,jde-1)
 
    IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
-   IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
+   IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = ite-2
 
    DO j = j_start, jtf
 
@@ -1574,7 +1574,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
    jtf=MIN(jte,jde-1)
 
    IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
-   IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
+   IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jte-3
 
    DO j = j_start, jtf
 
@@ -1664,7 +1664,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
    jtf=MIN(jte,jde-1)
 
    IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
-   IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
+   IF ( (config_flags%open_xe) .and. ite == ide ) itf = ite-3
 
    DO j = j_start, jtf
 

dyn_em/module_diffusion_em.F

@@ -2160,9 +2160,18 @@ SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv,         &
 ! tke_seed if the drag and flux are off.
 
     c_k = config_flags%c_k
-    tke_seed = tke_seed_value
-    if( (config_flags%tke_drag_coefficient .gt. epsilon) .or.  &
-        (config_flags%tke_heat_flux .gt. epsilon)  ) tke_seed = 0.
+    tke_seed = 0.
+    IF (config_flags%isfflx .eq. 0) THEN
+      IF ((config_flags%diff_opt .eq. 2) .and. (config_flags%bl_pbl_physics .eq. 0)) THEN
+        IF( (config_flags%tke_drag_coefficient .lt. epsilon) .and.  &
+            (config_flags%tke_heat_flux .lt. epsilon)  )  THEN
+          tke_seed = tke_seed_value
+        ENDIF
+      ELSE
+        !tke_drag_coefficient and tke_heat_flux are irrelevant here
+        tke_seed = tke_seed_value
+      ENDIF
+    ENDIF
 
     DO j = j_start, j_end
     DO k = kts+1, ktf-1
@@ -4282,11 +4291,13 @@ SUBROUTINE vertical_diffusion_2   ( ru_tendf, rv_tendf, rw_tendf, rt_tendf,   &
     DO i = i_start, i_end
        cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
        hfx(i,j)=heat_flux*cpm*rho(i,kts,j)         ! provided for output only
-       rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
-             -g*heat_flux*rho(i,kts,j)/dnw(kts)
        if(config_flags%use_theta_m == 1)THEN
          rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
+             -g*heat_flux*(1.+rvovrd*moist(i,kts,j,P_QV))*rho(i,kts,j)/dnw(kts) &
              -g*1.61*theta(i,kts,j)*qfx(i,j)/dnw(kts)
+       ELSE
+         rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
+             -g*heat_flux*rho(i,kts,j)/dnw(kts)
        ENDIF
     ENDDO
     ENDDO
@@ -4297,11 +4308,13 @@ SUBROUTINE vertical_diffusion_2   ( ru_tendf, rv_tendf, rw_tendf, rt_tendf,   &
 
        cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
        heat_flux = hfx(i,j)/cpm
-       rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
-            -g*heat_flux/dnw(kts)
        if(config_flags%use_theta_m == 1)THEN
          rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
+             -g*heat_flux*(1.+rvovrd*moist(i,kts,j,P_QV))/dnw(kts) &
              -g*1.61*theta(i,kts,j)*qfx(i,j)/dnw(kts)
+       ELSE
+         rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
+             -g*heat_flux/dnw(kts)
        ENDIF
 
     ENDDO
@@ -5133,7 +5146,7 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl,     &
   DO i = i_start, i_end
       delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j))
       pu1=pu(delxy,hpbl(i,j))
-      IF(sfcflg(i,j).and.sflux(i,j).GT.0.0)THEN
+      IF((sfcflg(i,j).and.sflux(i,j).GT.0.0) .and. (hpbl(i,j) .GT. 0)) THEN
 !nonlocal momentum flux based on Brown and Grant (1997)
        brint = -sm*ust(i,j)*ust(i,j)*wstar3(i,j)/(hpbl(i,j)*wscale(i,j)**4)
        gamu(i,j) = pu1 * brint*u_phy(i,1,j)/wspd(i,j)
@@ -5165,8 +5178,12 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl,     &
        deltaoh(i,j)  = max(ezfac*deltaoh(i,j),hpbl(i,j)-za(i,kpbl(i,j)-1,j)-1.)
        deltaoh(i,j)  = min(deltaoh(i,j), hpbl(i,j))
  
-       rigs(i,j)     = govrth(i,j)*dthv(i,j)*deltaoh(i,j)/(du**2.+dv**2.)
-       rigs(i,j)     = max(min(rigs(i,j), rigsmax),rimin)
+       if ((du .ne. 0) .or. (dv .ne. 0)) then
+           rigs(i,j) = govrth(i,j)*dthv(i,j)*deltaoh(i,j)/(du**2.+dv**2.)
+           rigs(i,j) = max(min(rigs(i,j), rigsmax),rimin)
+       else
+           rigs(i,j) = rigsmax
+       endif
        enlfrac2(i,j) = max(min(wm3/wstar3(i,j)/(1.+cpent/rigs(i,j)),entfmax),entfmin)
        enlfrac2(i,j) = enlfrac2(i,j)*enlfrac
      ENDIF
@@ -5182,44 +5199,41 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl,     &
   ENDDO
   ENDDO
   ENDDO
- 
-  DO j = j_start, j_end
-  DO i = i_start, i_end
-         deltaoh(i,j) = deltaoh(i,j)/hpbl(i,j)
-  ENDDO
-  ENDDO
 
   DO j = j_start, j_end
   DO i = i_start, i_end
-     delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j))
-     mlfrac        = mltop-deltaoh(i,j)
-     ezfrac        = mltop+deltaoh(i,j)
-     zfacmf(i,1,j) = min(max((zq(i,2,j)/hpbl(i,j)),zfmin),1.)
-     sfcfracn      = max(sfcfracn1,zfacmf(i,1,j))
-!
-     sflux0      = (a11+a12*sfcfracn)*sflux(i,j)
-     snlflux0    = nlfrac*sflux0
-     amf1        = snlflux0/sfcfracn
-     amf2        = -snlflux0/(mlfrac-sfcfracn)
-     bmf2        = -mlfrac*amf2
-     amf3        = snlflux0*enlfrac2(i,j)/deltaoh(i,j)
-     bmf3        = -amf3*mlfrac
-     hfxpbl(i,j) = amf3+bmf3
-     pth1        = pthnl(delxy,hpbl(i,j))
-     hfxpbl(i,j) = hfxpbl(i,j)*pth1
+     IF (pblflg(i,j)) THEN
+         deltaoh(i,j) = deltaoh(i,j)/hpbl(i,j)
+         delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j))
+         mlfrac        = mltop-deltaoh(i,j)
+         ezfrac        = mltop+deltaoh(i,j)
+         zfacmf(i,1,j) = min(max((zq(i,2,j)/hpbl(i,j)),zfmin),1.)
+         sfcfracn      = max(sfcfracn1,zfacmf(i,1,j))
+    !
+         sflux0      = (a11+a12*sfcfracn)*sflux(i,j)
+         snlflux0    = nlfrac*sflux0
+         amf1        = snlflux0/sfcfracn
+         amf2        = -snlflux0/(mlfrac-sfcfracn)
+         bmf2        = -mlfrac*amf2
+         amf3        = snlflux0*enlfrac2(i,j)/deltaoh(i,j)
+         bmf3        = -amf3*mlfrac
+         hfxpbl(i,j) = amf3+bmf3
+         pth1        = pthnl(delxy,hpbl(i,j))
+         hfxpbl(i,j) = hfxpbl(i,j)*pth1
 
-    DO k = kts, ktf
-       zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin)
-       IF(pblflg(i,j).and.k.LT.kpbl(i,j)) THEN
-         IF(zfacmf(i,k,j).LE.sfcfracn) THEN
-           nlflux(i,k,j) =  amf1*zfacmf(i,k,j)
-         ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN
-           nlflux(i,k,j) =  amf2*zfacmf(i,k,j)+bmf2
-         ENDIF
-         nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j)) 
-         nlflux(i,k,j) = nlflux(i,k,j)*pth1
-       ENDIF
-    ENDDO
+        DO k = kts, ktf
+           zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin)
+           IF(k.LT.kpbl(i,j)) THEN
+             IF(zfacmf(i,k,j).LE.sfcfracn) THEN
+               nlflux(i,k,j) =  amf1*zfacmf(i,k,j)
+             ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN
+               nlflux(i,k,j) =  amf2*zfacmf(i,k,j)+bmf2
+             ENDIF
+             nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j))
+             nlflux(i,k,j) = nlflux(i,k,j)*pth1
+           ENDIF
+        ENDDO
+    ENDIF
   ENDDO
   ENDDO
 END SUBROUTINE nonlocal_flux
@@ -5237,10 +5251,15 @@ FUNCTION pthnl(d,h)
    REAL,PARAMETER :: b1 = 2.0, b2 = 0.875
    real           :: d,h,doh,num,den
 
-   doh   = d/h
-   num   = a1*(doh)**b1 + a2*(doh)**b2+a3
-   den   = a4*(doh)**b1 + a5*(doh)**b2+a6
-   pthnl = a7*num/den + (1. - a7)
+   if (h .ne. 0) then
+      doh   = d/h
+      num   = a1*(doh)**b1 + a2*(doh)**b2+a3
+      den   = a4*(doh)**b1 + a5*(doh)**b2+a6
+      pthnl = a7*num/den + (1. - a7)
+   else
+      pthnl = 1.
+   endif
+
    pthnl = max(pthnl,pmin)
    pthnl = min(pthnl,pmax)
 
@@ -5261,10 +5280,14 @@ FUNCTION pthl(d,h)
    REAL,PARAMETER :: b1 = 2.0, b2 = 0.5
    REAL           :: d,h,doh,num,den
 
-   doh  = d/h
-   num  = a1*(doh)**b1 + a2*(doh)**b2+a3
-   den  = a4*(doh)**b1 + a5*(doh)**b2+a6
-   pthl = a7*num/den+(1. - a7)
+   if (h .ne. 0) then
+      doh  = d/h
+      num  = a1*(doh)**b1 + a2*(doh)**b2+a3
+      den  = a4*(doh)**b1 + a5*(doh)**b2+a6
+      pthl = a7*num/den+(1. - a7)
+   else
+      pthl = 1.
+   endif
    pthl = max(pthl,pmin)
    pthl = min(pthl,pmax)
 
@@ -5284,10 +5307,14 @@ FUNCTION pu(d,h)
    REAL,PARAMETER :: b1 = 2.0, b2 = 0.6666667
    REAL           :: d,h,doh,num,den
 
-   doh = d/h
-   num = a1*(doh)**b1 + a2*(doh)**b2
-   den = a3*(doh)**b1 + a4*(doh)**b2+a5
-   pu  = num/den
+   if (h .ne. 0) then
+      doh = d/h
+      num = a1*(doh)**b1 + a2*(doh)**b2
+      den = a3*(doh)**b1 + a4*(doh)**b2+a5
+      pu  = num/den
+   else
+      pu = 1.
+   endif
    pu  = max(pu,pmin)
    pu  = min(pu,pmax)
 
@@ -7717,7 +7744,7 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,&
    INTEGER :: i_start, i_end, j_start, j_end
 
    REAL :: V0_u,V0_v,ustar,beta
-   REAL :: heat_flux, moist_flux
+   REAL :: heat_flux, moist_flux, qfac
    REAL :: cpm,rdz_w,rhoavg_u,rhoavg_v,rdz_z
 ! End declarations.
 !-----------------------------------------------------------------------
@@ -8094,17 +8121,23 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,&
                   - dt*rdz_w*nlflux_rho(i,k+1,j)        &
                   + dt*rdz_w*hfx(i,j)/cpm
           ELSE
+             qfac = 1.+rvovrd*moist(i,kts,j,P_QV)
              d(k) = var_mix(i,kts,j)                    &
-                  - dt*rdz_w*nlflux_rho(i,k+1,j)        &
-                  + dt*rdz_w*(hfx(i,j)/cpm+1.61*theta(i,kts,j)*qfx(i,j))
+                  - dt*rdz_w*nlflux_rho(i,k+1,j)*qfac   &
+                  + dt*rdz_w*(qfac*hfx(i,j)/cpm+1.61*theta(i,kts,j)*qfx(i,j))
           ENDIF
 
+          IF(config_flags%use_theta_m == 1)THEN
+              qfac = 1.+rvovrd*moist(i,kts,j,P_QV)
+          ELSE
+              qfac = 1.
+          ENDIF
           DO  k = kts+1, ktf-1
             rdz_w = -g/dnw(k)/(c1h(k)*mu(i,j) + c2h(k))
             a(k)  = -rdz_w*rdz(i,k,j)*xkxavg(i,k,j)*dt
             b(k)  = 1.+rdz_w*(rdz(i,k+1,j)*xkxavg(i,k+1,j)+rdz(i,k,j)*xkxavg(i,k,j))*dt
             c(k)  = -rdz_w*rdz(i,k+1,j)*xkxavg(i,k+1,j)*dt
-            d(k)  = -rdz_w*(nlflux_rho(i,k+1,j)-nlflux_rho(i,k,j))*dt + var_mix(i,k,j)
+            d(k)  = -rdz_w*(nlflux_rho(i,k+1,j)-nlflux_rho(i,k,j))*qfac*dt + var_mix(i,k,j)
           ENDDO
 
             a(ktf) = 0.
@@ -8146,15 +8179,21 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,&
     DO  i = i_start, i_end
         DO  k = kts, ktf-1
            rdz_w = - g/dnw(k)/(c1h(k)*mu(i,j)+c2h(k))
-           beta  =   1.5*sqrt(tke(i,k,j))/l_diss(i,k,j)
+           IF (l_diss(i,k,j) .ne. 0) THEN
+               beta  =   1.5*sqrt(tke(i,k,j))/l_diss(i,k,j)
+           ELSE
+               beta  =   0.
+           ENDIF
            a(k)  = - 2.0*xkxavg(i,k,j)*dt*rdz_w*rdz(i,k,j)
            b(k)  =   1.0 + 2.0*dt*rdz_w*(rdz(i,k,j)*xkxavg(i,k,j)    &
                    + rdz(i,k+1,j)*xkxavg(i,k+1,j))                   &
                    + dt*beta
            c(k)  = - 2.0*xkxavg(i,k+1,j)*dt*rdz(i,k+1,j)*rdz_w
 
-           d(k)  =   tke(i,k,j)                                  &
-                   + 0.5*dt*tke(i,k,j)**1.5/l_diss(i,k,j)
+           d(k)  =   tke(i,k,j)
+           IF (l_diss(i,k,j) .ne. 0) THEN
+               d(k) = d(k) + 0.5*dt*tke(i,k,j)**1.5/l_diss(i,k,j)
+           ENDIF
         ENDDO
 
            a(ktf) = 0. !-1