Fixes for WDM6 and WDM7 (#1933)

TYPE: bug fixes

KEYWORDS: WDM schemes, uninitialized value, wrong cloud autoconversion rate

SOURCE: Songyou Hong, PSL/NOAA, internal

DESCRIPTION OF CHANGES:
Problem:
1. cloud autoconversion rate is off by a factor of 600
2. if cloud water is present, cloud number concentration needs to be set at initialization.

Solution:
Both problems are fixed in this PR. 

LIST OF MODIFIED FILES: 
M    phys/module_mp_wdm6.F
M    phys/module_mp_wdm7.F

TESTS CONDUCTED: 
1. Tested in case studies.
2. The Jenkins tests are all passing.

RELEASE NOTE: This PR fixes two errors in the WDM6 and WDM7 that have been in the code since V4.0: 1. uninitialized cloud number concentration when cloud is present at the start. 2. The cloud autoconversion rate was off by a factor of 600. The effect of the first error isn't large, but the effect of the wrong autoconversion rate has resulted in doubling the surface rainfall in some tests (warm rain processes).

phys/module_mp_wdm6.F

@@ -383,6 +383,12 @@ subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz            &
 !        nccn > 100 cm-3, bae from kiaps, oct 2017
 !        bug fix hydrometeros update before condensation process of rain,
 !           bae from kiaps, oct2017
+!        bug fix in snow melting and autoconversion.
+!        ==> Revisions enchance melting and autoconverison.
+!           hong from noaa, july 2023 
+!        remedy for tiny nc in the presence of qc (nc = xncr, 300/cm**3)
+!        ==> prevents the crash in computing autoconversion
+!           hong from noaa, july 2023
 !
 !  history log :
 !
@@ -586,7 +592,6 @@ subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz            &
        dend(i,k) = den(i,k)
      enddo
    enddo
-
 !
 !  latent heat for phase changes and heat capacity. neglect the
 !  changes during microphysical process calculation
@@ -904,8 +909,7 @@ subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz            &
               psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
                          *n0sfac(i,k)*(precs1*rslope2(i,k,2)                 &
                          +precs2*work2(i,k)*coeres)/den(i,k)                  
-              psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2)     &
-                         /mstep(i)),0.)
+              psmlt(i,k) = min(max(psmlt(i,k)*dtcld,-qrs(i,k,2)),0.)
 !-------------------------------------------------------------------
 ! nsmlt: melting of snow [LH A27]
 !       (T>T0: ->NR)
@@ -928,8 +932,7 @@ subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz            &
               pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1     &
                            *rslope2(i,k,3) + precg2*work2(i,k)*coeres)       &
                            /den(i,k)                                          
-              pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i),                &
-                          -qrs(i,k,3)/mstep(i)),0.)
+              pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld,-qrs(i,k,3)),0.)
 !-------------------------------------------------------------------
 ! ngmlt: melting of graupel [LH A28]
 !       (T>T0: ->NR)
@@ -1165,7 +1168,7 @@ subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz            &
          lencon  = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k)         &
                   *rslopec2(i,k)-0.4)
          lenconcr = max(1.2*lencon, qcrmin)
-         if(qci(i,k,1).gt.qcr(i,k)) then
+         if(qci(i,k,1).gt.qcr(i,k).and.ncr(i,k,2).gt.ncmin) then
            praut(i,k) = qck1*qci(i,k,1)**(7./3.)*ncr(i,k,2)**(-1./3.)
            praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
 !----------------------------------------------------------------------
@@ -2109,7 +2112,7 @@ SUBROUTINE wdm6init(den0,denr,dens,cl,cpv, ccn0, hail_opt, allowed_to_read) ! RA
 !
    qc0  = 4./3.*pi*denr*r0**3.*xncr0/den0
    qc1  = 4./3.*pi*denr*r0**3.*xncr1/den0
-   qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
+   qck1 = .104*9.8*peaut/(denr)**(1./3.)/xmyu*den0**(4./3.) ! 4706.08203
    pidnc = pi*denr/6.
 !
    bvtr1 = 1.+bvtr

phys/module_mp_wdm7.F

@@ -380,6 +380,16 @@ SUBROUTINE wdm72D(t, q, qci, qrs, ncr, den, p, delz             &
                                                                   )
 !-------------------------------------------------------------------
   IMPLICIT NONE
+!
+!        The code was first implemented in WRF in v4.1
+!        Further Development:
+!        bug fix in snow melting and autoconversion.
+!        ==> Revisions enchance melting and autoconverison.
+!           hong from noaa, july2023
+!        remedy for tiny nc in the presence of qc (nc = xncr, 300/cm**3)
+!        ==> prevents the crashed in computing autoconversion
+!           hong from noaa, july2023
+!
 !-------------------------------------------------------------------
   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
                                       ims,ime, jms,jme, kms,kme , &
@@ -916,8 +926,7 @@ SUBROUTINE wdm72D(t, q, qci, qrs, ncr, den, p, delz             &
               psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
                          *n0sfac(i,k)*(precs1*rslope2(i,k,2)                 &
                          +precs2*work2(i,k)*coeres)/den(i,k)                  
-              psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2)     &
-                         /mstep(i)),0.)
+              psmlt(i,k) = min(max(psmlt(i,k)*dtcld,-qrs(i,k,2)),0.)
 !
 ! nsmlt: melting of snow [LH A27]
 !       (T>T0: ->NR)
@@ -939,8 +948,7 @@ SUBROUTINE wdm72D(t, q, qci, qrs, ncr, den, p, delz             &
               pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1     &
                            *rslope2(i,k,3) + precg2*work2(i,k)*coeres)       &
                            /den(i,k)                                          
-              pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i),                &
-                          -qrs(i,k,3)/mstep(i)),0.)
+              pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld,-qrs(i,k,3)),0.)
 !
 ! ngmlt: melting of graupel [LH A28]
 !       (T>T0: ->NR)
@@ -1212,7 +1220,7 @@ SUBROUTINE wdm72D(t, q, qci, qrs, ncr, den, p, delz             &
           lencon  = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k)        &
                    *rslopec2(i,k)-0.4)
           lenconcr = max(1.2*lencon, qcrmin)
-          if(qci(i,k,1).gt.qcr(i,k)) then
+          if(qci(i,k,1).gt.qcr(i,k).and.ncr(i,k,2).gt.ncmin) then
             praut(i,k) = qck1*qci(i,k,1)**(7./3.)*ncr(i,k,2)**(-1./3.)
             praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
 !
@@ -2381,7 +2389,7 @@ SUBROUTINE wdm7init(den0,denr,dens,cl,cpv, ccn0, allowed_to_read) ! RAS
 !
    qc0  = 4./3.*pi*denr*r0**3.*xncr0/den0
    qc1  = 4./3.*pi*denr*r0**3.*xncr1/den0
-   qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
+   qck1 = .104*9.8*peaut/(denr)**(1./3.)/xmyu*den0**(4./3.) ! 4706.08203
    pidnc = pi*denr/6.
 !
    bvtr1 = 1.+bvtr