<p><b>laura@ucar.edu</b> 2010-07-23 14:55:58 -0600 (Fri, 23 Jul 2010)</p><p>Kessler cloud microphysics scheme<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_physics/src/core_physics/module_mp_kessler.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_mp_kessler.F (rev 0)
+++ branches/atmos_physics/src/core_physics/module_mp_kessler.F 2010-07-23 20:55:58 UTC (rev 400)
@@ -0,0 +1,244 @@
+!WRF:MODEL_LAYER:PHYSICS
+!
+
+MODULE module_mp_kessler
+
+CONTAINS
+!----------------------------------------------------------------
+ SUBROUTINE kessler( t, qv, qc, qr, rho, pii &
+ ,dt_in, z, xlv, cp &
+ ,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater &
+ ,dz8w &
+ ,RAINNC, RAINNCV &
+ ,ids,ide, jds,jde, kds,kde & ! domain dims
+ ,ims,ime, jms,jme, kms,kme & ! memory dims
+ ,its,ite, jts,jte, kts,kte & ! tile dims
+ )
+!----------------------------------------------------------------
+ IMPLICIT NONE
+!----------------------------------------------------------------
+ ! taken from the COMMAS code - WCS 10 May 1999.
+ ! converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999.
+!----------------------------------------------------------------
+ REAL , PARAMETER :: c1 = .001
+ REAL , PARAMETER :: c2 = .001
+ REAL , PARAMETER :: c3 = 2.2
+ REAL , PARAMETER :: c4 = .875
+ REAL , PARAMETER :: fudge = 1.0
+ REAL , PARAMETER :: mxfall = 10.0
+!----------------------------------------------------------------
+ INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte
+ REAL , INTENT(IN ) :: xlv, cp
+ REAL , INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
+ REAL , INTENT(IN ) :: rhowater
+
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(INOUT) :: &
+ t , &
+ qv, &
+ qc, &
+ qr
+
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(IN ) :: &
+ rho, &
+ pii, &
+ dz8w
+
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(IN ) :: z
+
+ REAL, INTENT(IN ) :: dt_in
+
+ REAL, DIMENSION( ims:ime , jms:jme ), &
+ INTENT(INOUT) :: RAINNC, &
+ RAINNCV
+
+
+
+ ! local variables
+
+ REAL :: qrprod, ern, gam, rcgs, rcgsi
+ REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: prod
+ REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw
+ INTEGER :: i,j,k
+ INTEGER :: nfall, n, nfall_new
+ REAL :: qrr, pressure, temp, es, qvs, dz, dt
+ REAL :: f5, dtfall, rdz, product
+ REAL :: max_heating, max_condense, max_rain, maxqrp
+ REAL :: vtmax, ernmax, crmax, factorn, time_sediment
+ REAL :: qcr, factorr, ppt
+ REAL, PARAMETER :: max_cr_sedimentation = 0.75
+!----------------------------------------------------------------
+
+ INTEGER :: imax, kmax
+
+ dt = dt_in
+
+! f5 = 237.3 * 17.27 * 2.5e6 / cp
+ f5 = svp2*(svpt0-svp3)*xlv/cp
+ ernmax = 0.
+ maxqrp = -100.
+
+!------------------------------------------------------------------------------
+! parameters for the time split terminal advection
+!------------------------------------------------------------------------------
+
+ max_heating = 0.
+ max_condense = 0.
+ max_rain = 0.
+
+!-----------------------------------------------------------------------------
+! outer J loop for entire microphysics, outer i loop for sedimentation
+!-----------------------------------------------------------------------------
+
+ microphysics_outer_j_loop: DO j = jts, jte
+
+ sedimentation_outer_i_loop: DO i = its,ite
+
+! vtmax = 0.
+ crmax = 0.
+
+
+!------------------------------------------------------------------------------
+! Terminal velocity calculation and advection, set up coefficients and
+! compute stable timestep
+!------------------------------------------------------------------------------
+
+ DO k = 1, kte
+ prodk(k) = qr(i,k,j)
+ rhok(k) = rho(i,k,j)
+ qrr = amax1(0.,qr(i,k,j)*0.001*rhok(k))
+ vtden(k) = sqrt(rhok(1)/rhok(k))
+ vt(k) = 36.34*(qrr**0.1364) * vtden(k)
+! vtmax = amax1(vt(k), vtmax)
+ rdzw(k) = 1./dz8w(i,k,j)
+ crmax = amax1(vt(k)*dt*rdzw(k),crmax)
+ ENDDO
+ DO k = 1, kte-1
+ rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j))
+ ENDDO
+ rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j))
+
+ nfall = max(1,nint(0.5+crmax/max_cr_sedimentation)) ! courant number for big timestep.
+ dtfall = dt / float(nfall) ! splitting so courant number for sedimentation
+ time_sediment = dt ! is stable
+
+!------------------------------------------------------------------------------
+! Terminal velocity calculation and advection
+! Do a time split loop on this for stability.
+!------------------------------------------------------------------------------
+
+ column_sedimentation: DO WHILE ( nfall > 0 )
+
+ time_sediment = time_sediment - dtfall
+ DO k = 1, kte-1
+ factor(k) = dtfall*rdzk(k)/rhok(k)
+ ENDDO
+ factor(kte) = dtfall*rdzk(kte)
+
+ ppt=0.
+
+ k = 1
+ ppt=rhok(k)*prodk(k)*vt(k)*dtfall/rhowater
+ RAINNCV(i,j)=ppt*1000.
+ RAINNC(i,j)=RAINNC(i,j)+ppt*1000. ! unit = mm
+
+!------------------------------------------------------------------------------
+! Time split loop, Fallout done with flux upstream
+!------------------------------------------------------------------------------
+
+ DO k = kts, kte-1
+ prodk(k) = prodk(k) - factor(k) &
+ * (rhok(k)*prodk(k)*vt(k) &
+ -rhok(k+1)*prodk(k+1)*vt(k+1))
+ ENDDO
+
+ k = kte
+ prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
+
+!------------------------------------------------------------------------------
+! compute new sedimentation velocity, and check/recompute new
+! sedimentation timestep if this isn't the last split step.
+!------------------------------------------------------------------------------
+
+ IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep
+
+ nfall = nfall - 1
+ crmax = 0.
+ DO k = kts, kte
+ qrr = amax1(0.,prodk(k)*0.001*rhok(k))
+ vt(k) = 36.34*(qrr**0.1364) * vtden(k)
+! vtmax = amax1(vt(k), vtmax)
+ crmax = amax1(vt(k)*time_sediment*rdzw(k),crmax)
+ ENDDO
+
+ nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation))
+ if (nfall_new /= nfall ) then
+ nfall = nfall_new
+ dtfall = time_sediment/nfall
+ end if
+
+ ELSE ! this was the last timestep
+
+ DO k=kts,kte
+ prod(i,k,j) = prodk(k)
+ ENDDO
+ nfall = 0 ! exit condition for sedimentation loop
+
+ END IF
+
+ ENDDO column_sedimentation
+
+ ENDDO sedimentation_outer_i_loop
+
+!------------------------------------------------------------------------------
+! Production of rain and deletion of qc
+! Production of qc from supersaturation
+! Evaporation of QR
+!------------------------------------------------------------------------------
+
+ DO k = kts, kte
+ DO i = its, ite
+ factorn = 1.0 / (1.+c3*dt*amax1(0.,qr(i,k,j))**c4)
+ qrprod = qc(i,k,j) * (1.0 - factorn) &
+ + factorn*c1*dt*amax1(qc(i,k,j)-c2,0.)
+ rcgs = 0.001*rho(i,k,j)
+
+ qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.)
+ qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j))
+ qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.)
+
+ temp = pii(i,k,j)*t(i,k,j)
+ pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.))
+ gam = 2.5e+06/(1004.*pii(i,k,j))
+! qvs = 380.*exp(17.27*(temp-273.)/(temp- 36.))/pressure
+ es = 1000.*svp1*exp(svp2*(temp-svpt0)/(temp-svp3))
+ qvs = ep2*es/(pressure-es)
+! prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvs*f5/(temp-36.)**2)
+ prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)*qvs*f5/(temp-svp3)**2)
+ ern = amin1(dt*(((1.6+124.9*(rcgs*qr(i,k,j))**.2046) &
+ *(rcgs*qr(i,k,j))**.525)/(2.55e8/(pressure*qvs) &
+ +5.4e5))*(dim(qvs,qv(i,k,j))/(rcgs*qvs)), &
+ amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j))
+
+! Update all variables
+
+ product = amax1(prod(i,k,j),-qc(i,k,j))
+ t (i,k,j) = t(i,k,j) + gam*(product - ern)
+ qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.)
+ qc(i,k,j) = qc(i,k,j) + product
+ qr(i,k,j) = qr(i,k,j) - ern
+
+ ENDDO
+ ENDDO
+
+ ENDDO microphysics_outer_j_loop
+
+ RETURN
+
+ END SUBROUTINE kessler
+
+END MODULE module_mp_kessler
</font>
</pre>