1 | MODULE dcmip2016_kessler_physic_mod |
---|
2 | |
---|
3 | |
---|
4 | |
---|
5 | CONTAINS |
---|
6 | !----------------------------------------------------------------------- |
---|
7 | ! |
---|
8 | ! Version: 2.0 |
---|
9 | ! |
---|
10 | ! Date: January 22nd, 2015 |
---|
11 | ! |
---|
12 | ! Change log: |
---|
13 | ! v2 - Added sub-cycling of rain sedimentation so as not to violate |
---|
14 | ! CFL condition. |
---|
15 | ! |
---|
16 | ! The KESSLER subroutine implements the Kessler (1969) microphysics |
---|
17 | ! parameterization as described by Soong and Ogura (1973) and Klemp |
---|
18 | ! and Wilhelmson (1978, KW). KESSLER is called at the end of each |
---|
19 | ! time step and makes the final adjustments to the potential |
---|
20 | ! temperature and moisture variables due to microphysical processes |
---|
21 | ! occurring during that time step. KESSLER is called once for each |
---|
22 | ! vertical column of grid cells. Increments are computed and added |
---|
23 | ! into the respective variables. The Kessler scheme contains three |
---|
24 | ! moisture categories: water vapor, cloud water (liquid water that |
---|
25 | ! moves with the flow), and rain water (liquid water that falls |
---|
26 | ! relative to the surrounding air). There are no ice categories. |
---|
27 | ! Variables in the column are ordered from the surface to the top. |
---|
28 | ! |
---|
29 | ! SUBROUTINE KESSLER(theta, qv, qc, qr, rho, pk, dt, z, nz, rainnc) |
---|
30 | ! |
---|
31 | ! Input variables: |
---|
32 | ! theta - potential temperature (K) |
---|
33 | ! qv - water vapor mixing ratio (gm/gm) |
---|
34 | ! qc - cloud water mixing ratio (gm/gm) |
---|
35 | ! qr - rain water mixing ratio (gm/gm) |
---|
36 | ! rho - dry air density (not mean state as in KW) (kg/m^3) |
---|
37 | ! pk - Exner function (not mean state as in KW) (p/p0)**(R/cp) |
---|
38 | ! dt - time step (s) |
---|
39 | ! z - heights of thermodynamic levels in the grid column (m) |
---|
40 | ! nz - number of thermodynamic levels in the column |
---|
41 | ! precl - Precipitation rate (m_water/s) |
---|
42 | ! |
---|
43 | ! Output variables: |
---|
44 | ! Increments are added into t, qv, qc, qr, and rainnc which are |
---|
45 | ! returned to the routine from which KESSLER was called. To obtain |
---|
46 | ! the total precip qt, after calling the KESSLER routine, compute: |
---|
47 | ! |
---|
48 | ! qt = sum over surface grid cells of (rainnc * cell area) (kg) |
---|
49 | ! [here, the conversion to kg uses (10^3 kg/m^3)*(10^-3 m/mm) = 1] |
---|
50 | ! |
---|
51 | ! |
---|
52 | ! Authors: Paul Ullrich |
---|
53 | ! University of California, Davis |
---|
54 | ! Email: paullrich@ucdavis.edu |
---|
55 | ! |
---|
56 | ! Based on a code by Joseph Klemp |
---|
57 | ! (National Center for Atmospheric Research) |
---|
58 | ! |
---|
59 | ! Reference: |
---|
60 | ! |
---|
61 | ! Klemp, J. B., W. C. Skamarock, W. C., and S.-H. Park, 2015: |
---|
62 | ! Idealized Global Nonhydrostatic Atmospheric Test Cases on a Reduced |
---|
63 | ! Radius Sphere. Journal of Advances in Modeling Earth Systems. |
---|
64 | ! doi:10.1002/2015MS000435 |
---|
65 | ! |
---|
66 | !======================================================================= |
---|
67 | |
---|
68 | SUBROUTINE KESSLER(theta, qv, qc, qr, rho, pk, dt, z, nz, precl) |
---|
69 | |
---|
70 | IMPLICIT NONE |
---|
71 | |
---|
72 | !------------------------------------------------ |
---|
73 | ! Input / output parameters |
---|
74 | !------------------------------------------------ |
---|
75 | INTEGER, INTENT(IN) :: nz ! Number of thermodynamic levels in the column |
---|
76 | REAL(8), DIMENSION(nz), INTENT(INOUT) :: & |
---|
77 | theta , & ! Potential temperature (K) |
---|
78 | qv , & ! Water vapor mixing ratio (gm/gm) |
---|
79 | qc , & ! Cloud water mixing ratio (gm/gm) |
---|
80 | qr ! Rain water mixing ratio (gm/gm) |
---|
81 | |
---|
82 | REAL(8), DIMENSION(nz), INTENT(IN) :: & |
---|
83 | rho ! Dry air density (not mean state as in KW) (kg/m^3) |
---|
84 | |
---|
85 | REAL(8), INTENT(OUT) :: & |
---|
86 | precl ! Precipitation rate (m_water / s) |
---|
87 | |
---|
88 | REAL(8), DIMENSION(nz), INTENT(IN) :: & |
---|
89 | z , & ! Heights of thermo. levels in the grid column (m) |
---|
90 | pk ! Exner function (p/p0)**(R/cp) |
---|
91 | |
---|
92 | REAL(8), INTENT(IN) :: & |
---|
93 | dt ! Time step (s) |
---|
94 | |
---|
95 | !------------------------------------------------ |
---|
96 | ! Local variables |
---|
97 | !------------------------------------------------ |
---|
98 | REAL, DIMENSION(nz) :: r, rhalf, velqr, sed, pc |
---|
99 | |
---|
100 | REAL(8) :: f5, f2x, xk, ern, qrprod, prod, qvs, psl, rhoqr, dt_max, dt0 |
---|
101 | |
---|
102 | INTEGER :: k, rainsplit, nt |
---|
103 | |
---|
104 | !------------------------------------------------ |
---|
105 | ! Begin calculation |
---|
106 | !------------------------------------------------ |
---|
107 | f2x = 17.27d0 |
---|
108 | f5 = 237.3d0 * f2x * 2500000.d0 / 1003.d0 |
---|
109 | xk = .2875d0 ! kappa (r/cp) |
---|
110 | psl = 1000.d0 ! pressure at sea level (mb) |
---|
111 | rhoqr = 1000.d0 ! density of liquid water (kg/m^3) |
---|
112 | |
---|
113 | do k=1,nz |
---|
114 | r(k) = 0.001d0*rho(k) |
---|
115 | rhalf(k) = sqrt(rho(1)/rho(k)) |
---|
116 | pc(k) = 3.8d0/(pk(k)**(1./xk)*psl) |
---|
117 | |
---|
118 | ! Liquid water terminal velocity (m/s) following KW eq. 2.15 |
---|
119 | velqr(k) = 36.34d0*(qr(k)*r(k))**0.1364*rhalf(k) |
---|
120 | |
---|
121 | end do |
---|
122 | |
---|
123 | ! Maximum time step size in accordance with CFL condition |
---|
124 | if (dt .le. 0.d0) then |
---|
125 | write(*,*) 'kessler.f90 called with nonpositive dt' |
---|
126 | stop |
---|
127 | end if |
---|
128 | |
---|
129 | dt_max = dt |
---|
130 | do k=1,nz-1 |
---|
131 | if (velqr(k) .ne. 0.d0) then |
---|
132 | dt_max = min(dt_max, 0.8d0*(z(k+1)-z(k))/velqr(k)) |
---|
133 | end if |
---|
134 | end do |
---|
135 | |
---|
136 | ! Number of subcycles |
---|
137 | rainsplit = ceiling(dt / dt_max) |
---|
138 | dt0 = dt / real(rainsplit,8) |
---|
139 | |
---|
140 | ! Subcycle through rain process |
---|
141 | precl = 0.d0 |
---|
142 | |
---|
143 | do nt=1,rainsplit |
---|
144 | |
---|
145 | ! Precipitation rate (m/s) |
---|
146 | precl = precl + rho(1) * qr(1) * velqr(1) / rhoqr |
---|
147 | |
---|
148 | ! Sedimentation term using upstream differencing |
---|
149 | do k=1,nz-1 |
---|
150 | sed(k) = dt0*(r(k+1)*qr(k+1)*velqr(k+1)-r(k)*qr(k)*velqr(k))/(r(k)*(z(k+1)-z(k))) |
---|
151 | end do |
---|
152 | sed(nz) = -dt0*qr(nz)*velqr(nz)/(.5*(z(nz)-z(nz-1))) |
---|
153 | |
---|
154 | ! Adjustment terms |
---|
155 | do k=1,nz |
---|
156 | |
---|
157 | ! Autoconversion and accretion rates following KW eq. 2.13a,b |
---|
158 | qrprod = qc(k) - (qc(k)-dt0*amax1(.001*(qc(k)-.001d0),0.))/(1.d0+dt0*2.2d0*qr(k)**.875) |
---|
159 | qc(k) = amax1(qc(k)-qrprod,0.) |
---|
160 | qr(k) = amax1(qr(k)+qrprod+sed(k),0.) |
---|
161 | |
---|
162 | ! Saturation vapor mixing ratio (gm/gm) following KW eq. 2.11 |
---|
163 | qvs = pc(k)*exp(f2x*(pk(k)*theta(k)-273.d0) & |
---|
164 | /(pk(k)*theta(k)- 36.d0)) |
---|
165 | prod = (qv(k)-qvs)/(1.d0+qvs*f5/(pk(k)*theta(k)-36.d0)**2) |
---|
166 | |
---|
167 | ! Evaporation rate following KW eq. 2.14a,b |
---|
168 | ern = amin1(dt0*(((1.6d0+124.9d0*(r(k)*qr(k))**.2046) & |
---|
169 | *(r(k)*qr(k))**.525)/(2550000d0*pc(k) & |
---|
170 | /(3.8d0 *qvs)+540000d0))*(dim(qvs,qv(k)) & |
---|
171 | /(r(k)*qvs)),amax1(-prod-qc(k),0.),qr(k)) |
---|
172 | |
---|
173 | ! Saturation adjustment following KW eq. 3.10 |
---|
174 | theta(k)= theta(k) + 2500000d0/(1003.d0*pk(k))*(amax1( prod,-qc(k))-ern) |
---|
175 | qv(k) = amax1(qv(k)-max(prod,-qc(k))+ern,0.) |
---|
176 | qc(k) = qc(k)+max(prod,-qc(k)) |
---|
177 | qr(k) = qr(k)-ern |
---|
178 | end do |
---|
179 | |
---|
180 | ! Recalculate liquid water terminal velocity |
---|
181 | if (nt .ne. rainsplit) then |
---|
182 | do k=1,nz |
---|
183 | velqr(k) = 36.34d0*(qr(k)*r(k))**0.1364*rhalf(k) |
---|
184 | end do |
---|
185 | end if |
---|
186 | end do |
---|
187 | |
---|
188 | precl = precl / dble(rainsplit) |
---|
189 | |
---|
190 | END SUBROUTINE KESSLER |
---|
191 | |
---|
192 | !======================================================================= |
---|
193 | |
---|
194 | |
---|
195 | |
---|
196 | END MODULE dcmip2016_kessler_physic_mod |
---|