Line data Source code
1 : module zm_conv_momtran
2 :
3 : use ccpp_kinds, only: kind_phys
4 :
5 : implicit none
6 :
7 : save
8 : private ! Make default type private to the module
9 : public zm_conv_momtran_run ! convective momentum transport
10 :
11 :
12 : contains
13 :
14 : !===============================================================================
15 : !> \section arg_table_zm_conv_momtran_run Argument Table
16 : !! \htmlinclude zm_conv_momtran_run.html
17 : !!
18 65016 : subroutine zm_conv_momtran_run(ncol, pver, pverp, &
19 65016 : domomtran,windu, windv,num_winds, mu, md, &
20 : momcu, momcd, &
21 65016 : du, eu, ed, dp, dsubcld , &
22 65016 : jt, mx, ideep , il1g, il2g, &
23 65016 : nstep, windu_tend, windv_tend, pguallu, pguallv, pgdallu, pgdallv, &
24 130032 : icwuu, icwuv, icwdu, icwdv, dt, seten)
25 : !-----------------------------------------------------------------------
26 : !
27 : ! Purpose:
28 : ! Convective transport of momentum
29 : !
30 : ! Mixing ratios may be with respect to either dry or moist air
31 : !
32 : ! Method:
33 : ! Based on the convtran subroutine by P. Rasch
34 : ! <Also include any applicable external references.>
35 : !
36 : ! Author: J. Richter and P. Rasch
37 : !
38 : !-----------------------------------------------------------------------
39 : ! CACNOTE - use CCPP constituents object
40 : use constituents, only: cnst_get_type_byind
41 :
42 : implicit none
43 : !-----------------------------------------------------------------------
44 : !
45 : ! Input arguments
46 : !
47 : integer, intent(in) :: ncol ! number of atmospheric columns
48 : integer, intent(in) :: num_winds ! number of wind directions
49 : integer, intent(in) :: pver, pverp
50 : logical, intent(in) :: domomtran(:) ! flag for doing convective transport (num_winds)
51 : real(kind_phys), intent(in) :: windu(:,:) ! U Wind array (ncol,pver)
52 : real(kind_phys), intent(in) :: windv(:,:) ! V Wind array (ncol,pver)
53 : real(kind_phys), intent(in) :: mu(:,:) ! Mass flux up (ncol,pver)
54 : real(kind_phys), intent(in) :: md(:,:) ! Mass flux down (ncol,pver)
55 : real(kind_phys), intent(in) :: momcu
56 : real(kind_phys), intent(in) :: momcd
57 : real(kind_phys), intent(in) :: du(:,:) ! Mass detraining from updraft (ncol,pver)
58 : real(kind_phys), intent(in) :: eu(:,:) ! Mass entraining from updraft (ncol,pver)
59 : real(kind_phys), intent(in) :: ed(:,:) ! Mass entraining from downdraft (ncol,pver)
60 : real(kind_phys), intent(in) :: dp(:,:) ! Delta pressure between interfaces (ncol,pver)
61 : real(kind_phys), intent(in) :: dsubcld(:) ! Delta pressure from cloud base to sfc (ncol)
62 : real(kind_phys), intent(in) :: dt ! time step in seconds : 2*delta_t
63 :
64 : integer, intent(in) :: jt(:) ! Index of cloud top for each column (ncol)
65 : integer, intent(in) :: mx(:) ! Index of cloud top for each column (ncol)
66 : integer, intent(in) :: ideep(:) ! Gathering array (ncol)
67 : integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
68 : integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
69 : integer, intent(in) :: nstep ! Time step index
70 :
71 :
72 :
73 : ! input/output
74 :
75 : real(kind_phys), intent(out) :: windu_tend(:,:) ! U wind tendency
76 : real(kind_phys), intent(out) :: windv_tend(:,:) ! V wind tendency
77 :
78 : !--------------------------Local Variables------------------------------
79 :
80 : integer i ! Work index
81 : integer k ! Work index
82 : integer kbm ! Highest altitude index of cloud base
83 : integer kk ! Work index
84 : integer kkp1 ! Work index
85 : integer kkm1 ! Work index
86 : integer km1 ! Work index
87 : integer kp1 ! Work index
88 : integer ktm ! Highest altitude index of cloud top
89 : integer m ! Work index
90 : integer ii ! Work index
91 :
92 : real(kind_phys) cabv ! Mix ratio of constituent above
93 : real(kind_phys) cbel ! Mix ratio of constituent below
94 : real(kind_phys) cdifr ! Normalized diff between cabv and cbel
95 130032 : real(kind_phys) chat(ncol,pver) ! Mix ratio in env at interfaces
96 130032 : real(kind_phys) cond(ncol,pver) ! Mix ratio in downdraft at interfaces
97 130032 : real(kind_phys) const(ncol,pver) ! Gathered wind array
98 130032 : real(kind_phys) conu(ncol,pver) ! Mix ratio in updraft at interfaces
99 130032 : real(kind_phys) dcondt(ncol,pver) ! Gathered tend array
100 : real(kind_phys) mbsth ! Threshold for mass fluxes
101 : real(kind_phys) mupdudp ! A work variable
102 : real(kind_phys) minc ! A work variable
103 : real(kind_phys) maxc ! A work variable
104 : real(kind_phys) fluxin ! A work variable
105 : real(kind_phys) fluxout ! A work variable
106 : real(kind_phys) netflux ! A work variable
107 :
108 :
109 : real(kind_phys) sum ! sum
110 : real(kind_phys) sum2 ! sum2
111 :
112 130032 : real(kind_phys) mududp(ncol,pver) ! working variable
113 130032 : real(kind_phys) mddudp(ncol,pver) ! working variable
114 :
115 130032 : real(kind_phys) pgu(ncol,pver) ! Pressure gradient term for updraft
116 130032 : real(kind_phys) pgd(ncol,pver) ! Pressure gradient term for downdraft
117 :
118 : real(kind_phys),intent(out) :: pguallu(:,:) ! Apparent force from updraft PG on U winds ! (ncol,pver)
119 : real(kind_phys),intent(out) :: pguallv(:,:) ! Apparent force from updraft PG on V winds ! (ncol,pver)
120 : real(kind_phys),intent(out) :: pgdallu(:,:) ! Apparent force from downdraft PG on U winds! (ncol,pver)
121 : real(kind_phys),intent(out) :: pgdallv(:,:) ! Apparent force from downdraft PG on V winds! (ncol,pver)
122 :
123 : real(kind_phys),intent(out) :: icwuu(:,:) ! In-cloud U winds in updraft ! (ncol,pver)
124 : real(kind_phys),intent(out) :: icwuv(:,:) ! In-cloud V winds in updraft ! (ncol,pver)
125 : real(kind_phys),intent(out) :: icwdu(:,:) ! In-cloud U winds in downdraft ! (ncol,pver)
126 : real(kind_phys),intent(out) :: icwdv(:,:) ! In-cloud V winds in downdraft ! (ncol,pver)
127 :
128 : real(kind_phys),intent(out) :: seten(:,:) ! Dry static energy tendency ! (ncol,pver)
129 130032 : real(kind_phys) gseten(ncol,pver) ! Gathered dry static energy tendency
130 :
131 130032 : real(kind_phys) :: winds(ncol,pver,num_winds) ! combined winds array
132 130032 : real(kind_phys) :: wind_tends(ncol,pver,num_winds) ! combined tendency array
133 130032 : real(kind_phys) :: pguall(ncol,pver,num_winds) ! Combined apparent force from updraft PG on U winds
134 130032 : real(kind_phys) :: pgdall(ncol,pver,num_winds) ! Combined apparent force from downdraft PG on U winds
135 130032 : real(kind_phys) :: icwu(ncol,pver,num_winds) ! Combined In-cloud winds in updraft
136 130032 : real(kind_phys) :: icwd(ncol,pver,num_winds) ! Combined In-cloud winds in downdraft
137 :
138 130032 : real(kind_phys) mflux(ncol,pverp,num_winds) ! Gathered momentum flux
139 :
140 130032 : real(kind_phys) wind0(ncol,pver,num_winds) ! gathered wind before time step
141 130032 : real(kind_phys) windf(ncol,pver,num_winds) ! gathered wind after time step
142 : real(kind_phys) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2
143 :
144 :
145 : !-----------------------------------------------------------------------
146 : !
147 : ! Combine winds in single array
148 101092320 : winds(:,:,1) = windu(:,:)
149 101027304 : winds(:,:,2) = windv(:,:)
150 :
151 : ! Initialize outgoing fields
152 202119624 : pguall(:,:,:) = 0.0_kind_phys
153 202119624 : pgdall(:,:,:) = 0.0_kind_phys
154 : ! Initialize in-cloud winds to environmental wind
155 202119624 : icwu(:ncol,:,:) = winds(:ncol,:,:)
156 202119624 : icwd(:ncol,:,:) = winds(:ncol,:,:)
157 :
158 : ! Initialize momentum flux and final winds
159 204290856 : mflux(:,:,:) = 0.0_kind_phys
160 202119624 : wind0(:,:,:) = 0.0_kind_phys
161 202119624 : windf(:,:,:) = 0.0_kind_phys
162 :
163 : ! Initialize dry static energy
164 :
165 101027304 : seten(:,:) = 0.0_kind_phys
166 101027304 : gseten(:,:) = 0.0_kind_phys
167 :
168 : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
169 65016 : mbsth = 1.e-15_kind_phys
170 :
171 : ! Find the highest level top and bottom levels of convection
172 65016 : ktm = pver
173 65016 : kbm = pver
174 247119 : do i = il1g, il2g
175 182103 : ktm = min(ktm,jt(i))
176 247119 : kbm = min(kbm,mx(i))
177 : end do
178 :
179 : ! Loop ever each wind component
180 195048 : do m = 1, num_winds !start at m = 1 to transport momentum
181 195048 : if (domomtran(m)) then
182 :
183 : ! Gather up the winds and set tend to zero
184 12223008 : do k = 1,pver
185 46094166 : do i =il1g,il2g
186 33871158 : const(i,k) = winds(ideep(i),k,m)
187 45964134 : wind0(i,k,m) = const(i,k)
188 : end do
189 : end do
190 :
191 :
192 : ! From now on work only with gathered data
193 :
194 : ! Interpolate winds to interfaces
195 :
196 12223008 : do k = 1,pver
197 12092976 : km1 = max(1,k-1)
198 46094166 : do i = il1g, il2g
199 :
200 : ! use arithmetic mean
201 33871158 : chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
202 :
203 : ! Provisional up and down draft values
204 33871158 : conu(i,k) = chat(i,k)
205 33871158 : cond(i,k) = chat(i,k)
206 :
207 : ! provisional tends
208 45964134 : dcondt(i,k) = 0._kind_phys
209 :
210 : end do
211 : end do
212 :
213 :
214 : !
215 : ! Pressure Perturbation Term
216 : !
217 :
218 : !Top boundary: assume mu is zero
219 :
220 130032 : k=1
221 624270 : pgu(:il2g,k) = 0.0_kind_phys
222 494238 : pgd(:il2g,k) = 0.0_kind_phys
223 :
224 11962944 : do k=2,pver-1
225 11832912 : km1 = max(1,k-1)
226 11832912 : kp1 = min(pver,k+1)
227 45105690 : do i = il1g,il2g
228 :
229 : !interior points
230 :
231 99428238 : mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
232 132570984 : + mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
233 :
234 33142746 : pgu(i,k) = - momcu * 0.5_kind_phys * mududp(i,k)
235 :
236 :
237 33142746 : mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
238 33142746 : + md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
239 :
240 44975658 : pgd(i,k) = - momcd * 0.5_kind_phys * mddudp(i,k)
241 :
242 :
243 : end do
244 : end do
245 :
246 : ! bottom boundary
247 130032 : k = pver
248 130032 : km1 = max(1,k-1)
249 494238 : do i=il1g,il2g
250 :
251 364206 : mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
252 364206 : pgu(i,k) = - momcu * mududp(i,k)
253 :
254 364206 : mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
255 :
256 494238 : pgd(i,k) = - momcd * mddudp(i,k)
257 :
258 : end do
259 :
260 :
261 : !
262 : ! In-cloud velocity calculations
263 : !
264 :
265 : ! Do levels adjacent to top and bottom
266 494238 : k = 2
267 494238 : km1 = 1
268 494238 : kk = pver
269 494238 : kkm1 = max(1,kk-1)
270 494238 : do i = il1g,il2g
271 364206 : mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
272 364206 : if (mupdudp > mbsth) then
273 :
274 0 : conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
275 : endif
276 494238 : if (md(i,k) < -mbsth) then
277 0 : cond(i,k) = (-ed(i,km1)*const(i,km1)*dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k)
278 : endif
279 :
280 :
281 : end do
282 :
283 :
284 :
285 : ! Updraft from bottom to top
286 12092976 : do kk = pver-1,1,-1
287 11962944 : kkm1 = max(1,kk-1)
288 11962944 : kkp1 = min(pver,kk+1)
289 45599928 : do i = il1g,il2g
290 33506952 : mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
291 45469896 : if (mupdudp > mbsth) then
292 :
293 11159912 : conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* &
294 11159912 : const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
295 : endif
296 : end do
297 :
298 : end do
299 :
300 :
301 : ! Downdraft from top to bottom
302 11962944 : do k = 3,pver
303 11832912 : km1 = max(1,k-1)
304 45105690 : do i = il1g,il2g
305 44975658 : if (md(i,k) < -mbsth) then
306 :
307 8556560 : cond(i,k) = ( md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) &
308 8556560 : *dp(i,km1)-pgd(i,km1)*dp(i,km1) )/md(i,k)
309 :
310 : endif
311 : end do
312 : end do
313 :
314 :
315 : sum = 0._kind_phys
316 : sum2 = 0._kind_phys
317 :
318 :
319 1293600 : do k = ktm,pver
320 1163568 : km1 = max(1,k-1)
321 1163568 : kp1 = min(pver,k+1)
322 9799444 : do i = il1g,il2g
323 8505844 : ii = ideep(i)
324 :
325 : ! version 1 hard to check for roundoff errors
326 8505844 : dcondt(i,k) = &
327 8505844 : +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
328 8505844 : -mu(i,k)* (conu(i,k)-chat(i,k)) &
329 8505844 : +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
330 8505844 : -md(i,k)* (cond(i,k)-chat(i,k)) &
331 26681100 : )/dp(i,k)
332 :
333 : end do
334 : end do
335 :
336 : ! dcont for bottom layer
337 : !
338 510682 : do k = kbm,pver
339 2657928 : km1 = max(1,k-1)
340 2787960 : do i = il1g,il2g
341 2657928 : if (k == mx(i)) then
342 :
343 : ! version 1
344 364206 : dcondt(i,k) = (1._kind_phys/dp(i,k))* &
345 364206 : (-mu(i,k)*(conu(i,k)-chat(i,k)) &
346 364206 : -md(i,k)*(cond(i,k)-chat(i,k)) &
347 728412 : )
348 : end if
349 : end do
350 : end do
351 :
352 : ! Initialize to zero everywhere, then scatter tendency back to full array
353 202054608 : wind_tends(:,:,m) = 0._kind_phys
354 :
355 12223008 : do k = 1,pver
356 46094166 : do i = il1g,il2g
357 33871158 : ii = ideep(i)
358 33871158 : wind_tends(ii,k,m) = dcondt(i,k)
359 : ! Output apparent force on the mean flow from pressure gradient
360 33871158 : pguall(ii,k,m) = -pgu(i,k)
361 33871158 : pgdall(ii,k,m) = -pgd(i,k)
362 33871158 : icwu(ii,k,m) = conu(i,k)
363 45964134 : icwd(ii,k,m) = cond(i,k)
364 : end do
365 : end do
366 :
367 : ! Calculate momentum flux in units of mb*m/s2
368 :
369 1293600 : do k = ktm,pver
370 9799444 : do i = il1g,il2g
371 8505844 : ii = ideep(i)
372 8505844 : mflux(i,k,m) = &
373 8505844 : -mu(i,k)* (conu(i,k)-chat(i,k)) &
374 18175256 : -md(i,k)* (cond(i,k)-chat(i,k))
375 : end do
376 : end do
377 :
378 :
379 : ! Calculate winds at the end of the time step
380 :
381 1293600 : do k = ktm,pver
382 9799444 : do i = il1g,il2g
383 8505844 : ii = ideep(i)
384 8505844 : km1 = max(1,k-1)
385 8505844 : kp1 = k+1
386 9669412 : windf(i,k,m) = const(i,k) - (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k)
387 :
388 : end do
389 : end do
390 :
391 : end if ! for domomtran
392 : end do
393 :
394 : ! Need to add an energy fix to account for the dissipation of kinetic energy
395 : ! Formulation follows from Boville and Bretherton (2003)
396 : ! formulation by PJR
397 :
398 646800 : do k = ktm,pver
399 581784 : km1 = max(1,k-1)
400 581784 : kp1 = min(pver,k+1)
401 4899722 : do i = il1g,il2g
402 :
403 4252922 : ii = ideep(i)
404 :
405 : ! calculate the KE fluxes at top and bot of layer
406 : ! based on a discrete approximation to b&b eq(35) F_KE = u*F_u + v*F_v at interface
407 4252922 : utop = (wind0(i,k,1)+wind0(i,km1,1))/2._kind_phys
408 4252922 : vtop = (wind0(i,k,2)+wind0(i,km1,2))/2._kind_phys
409 4252922 : ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2._kind_phys
410 4252922 : vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2._kind_phys
411 4252922 : fket = utop*mflux(i,k,1) + vtop*mflux(i,k,2) ! top of layer
412 4252922 : fkeb = ubot*mflux(i,k+1,1) + vbot*mflux(i,k+1,2) ! bot of layer
413 :
414 : ! divergence of these fluxes should give a conservative redistribution of KE
415 4252922 : ketend_cons = (fket-fkeb)/dp(i,k)
416 :
417 : ! tendency in kinetic energy resulting from the momentum transport
418 4252922 : ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))*0.5_kind_phys/dt
419 :
420 : ! the difference should be the dissipation
421 4252922 : gset2 = ketend_cons - ketend
422 4834706 : gseten(i,k) = gset2
423 :
424 : end do
425 :
426 : end do
427 :
428 : ! Scatter dry static energy to full array
429 6111504 : do k = 1,pver
430 23047083 : do i = il1g,il2g
431 16935579 : ii = ideep(i)
432 22982067 : seten(ii,k) = gseten(i,k)
433 :
434 : end do
435 : end do
436 :
437 : ! Split out the wind tendencies
438 101092320 : windu_tend(:,:) = wind_tends(:,:,1)
439 101092320 : windv_tend(:,:) = wind_tends(:,:,2)
440 :
441 101092320 : pguallu(:,:) = pguall(:,:,1)
442 101092320 : pguallv(:,:) = pguall(:,:,2)
443 101092320 : pgdallu(:,:) = pgdall(:,:,1)
444 101092320 : pgdallv(:,:) = pgdall(:,:,2)
445 101092320 : icwuu(:ncol,:) = icwu(:,:,1)
446 101092320 : icwuv(:ncol,:) = icwu(:,:,2)
447 101092320 : icwdu(:ncol,:) = icwd(:,:,1)
448 101092320 : icwdv(:ncol,:) = icwd(:,:,2)
449 :
450 65016 : return
451 : end subroutine zm_conv_momtran_run
452 :
453 :
454 : end module zm_conv_momtran
|
