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 1495368 : subroutine zm_conv_momtran_run(ncol, pver, pverp, &
19 1495368 : domomtran,windu, windv,num_winds, mu, md, &
20 : momcu, momcd, &
21 1495368 : du, eu, ed, dp, dsubcld , &
22 1495368 : jt, mx, ideep , il1g, il2g, &
23 1495368 : nstep, windu_tend, windv_tend, pguallu, pguallv, pgdallu, pgdallv, &
24 2990736 : 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
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 2990736 : real(kind_phys) chat(ncol,pver) ! Mix ratio in env at interfaces
96 2990736 : real(kind_phys) cond(ncol,pver) ! Mix ratio in downdraft at interfaces
97 2990736 : real(kind_phys) const(ncol,pver) ! Gathered wind array
98 2990736 : real(kind_phys) conu(ncol,pver) ! Mix ratio in updraft at interfaces
99 2990736 : 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 2990736 : real(kind_phys) mududp(ncol,pver) ! working variable
113 2990736 : real(kind_phys) mddudp(ncol,pver) ! working variable
114 :
115 2990736 : real(kind_phys) pgu(ncol,pver) ! Pressure gradient term for updraft
116 2990736 : 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 2990736 : real(kind_phys) gseten(ncol,pver) ! Gathered dry static energy tendency
130 :
131 2990736 : real(kind_phys) :: winds(ncol,pver,num_winds) ! combined winds array
132 2990736 : real(kind_phys) :: wind_tends(ncol,pver,num_winds) ! combined tendency array
133 2990736 : real(kind_phys) :: pguall(ncol,pver,num_winds) ! Combined apparent force from updraft PG on U winds
134 2990736 : real(kind_phys) :: pgdall(ncol,pver,num_winds) ! Combined apparent force from downdraft PG on U winds
135 2990736 : real(kind_phys) :: icwu(ncol,pver,num_winds) ! Combined In-cloud winds in updraft
136 2990736 : real(kind_phys) :: icwd(ncol,pver,num_winds) ! Combined In-cloud winds in downdraft
137 :
138 2990736 : real(kind_phys) mflux(ncol,pverp,num_winds) ! Gathered momentum flux
139 :
140 2990736 : real(kind_phys) wind0(ncol,pver,num_winds) ! gathered wind before time step
141 2990736 : 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 652189104 : winds(:,:,1) = windu(:,:)
149 650693736 : winds(:,:,2) = windv(:,:)
150 :
151 : ! Initialize outgoing fields
152 1302882840 : pguall(:,:,:) = 0.0_kind_phys
153 1302882840 : pgdall(:,:,:) = 0.0_kind_phys
154 : ! Initialize in-cloud winds to environmental wind
155 1302882840 : icwu(:ncol,:,:) = winds(:ncol,:,:)
156 1302882840 : icwd(:ncol,:,:) = winds(:ncol,:,:)
157 :
158 : ! Initialize momentum flux and final winds
159 1352821176 : mflux(:,:,:) = 0.0_kind_phys
160 1302882840 : wind0(:,:,:) = 0.0_kind_phys
161 1302882840 : windf(:,:,:) = 0.0_kind_phys
162 :
163 : ! Initialize dry static energy
164 :
165 650693736 : seten(:,:) = 0.0_kind_phys
166 650693736 : gseten(:,:) = 0.0_kind_phys
167 :
168 : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
169 1495368 : mbsth = 1.e-15_kind_phys
170 :
171 : ! Find the highest level top and bottom levels of convection
172 1495368 : ktm = pver
173 1495368 : kbm = pver
174 13077368 : do i = il1g, il2g
175 11582000 : ktm = min(ktm,jt(i))
176 13077368 : kbm = min(kbm,mx(i))
177 : end do
178 :
179 : ! Loop ever each wind component
180 4486104 : do m = 1, num_winds !start at m = 1 to transport momentum
181 4486104 : if (domomtran(m)) then
182 :
183 : ! Gather up the winds and set tend to zero
184 80749872 : do k = 1,pver
185 683013872 : do i =il1g,il2g
186 602264000 : const(i,k) = winds(ideep(i),k,m)
187 680023136 : 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 80749872 : do k = 1,pver
197 77759136 : km1 = max(1,k-1)
198 683013872 : do i = il1g, il2g
199 :
200 : ! use arithmetic mean
201 602264000 : chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
202 :
203 : ! Provisional up and down draft values
204 602264000 : conu(i,k) = chat(i,k)
205 602264000 : cond(i,k) = chat(i,k)
206 :
207 : ! provisional tends
208 680023136 : 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 2990736 : k=1
221 29145472 : pgu(:il2g,k) = 0.0_kind_phys
222 26154736 : pgd(:il2g,k) = 0.0_kind_phys
223 :
224 74768400 : do k=2,pver-1
225 71777664 : km1 = max(1,k-1)
226 71777664 : kp1 = min(pver,k+1)
227 630704400 : do i = il1g,il2g
228 :
229 : !interior points
230 :
231 1667808000 : mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
232 2223744000 : + mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
233 :
234 555936000 : pgu(i,k) = - momcu * 0.5_kind_phys * mududp(i,k)
235 :
236 :
237 555936000 : mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
238 555936000 : + md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
239 :
240 627713664 : pgd(i,k) = - momcd * 0.5_kind_phys * mddudp(i,k)
241 :
242 :
243 : end do
244 : end do
245 :
246 : ! bottom boundary
247 2990736 : k = pver
248 2990736 : km1 = max(1,k-1)
249 26154736 : do i=il1g,il2g
250 :
251 23164000 : mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
252 23164000 : pgu(i,k) = - momcu * mududp(i,k)
253 :
254 23164000 : mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
255 :
256 26154736 : 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 26154736 : k = 2
267 26154736 : km1 = 1
268 26154736 : kk = pver
269 26154736 : kkm1 = max(1,kk-1)
270 26154736 : do i = il1g,il2g
271 23164000 : mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
272 23164000 : if (mupdudp > mbsth) then
273 :
274 22223084 : conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
275 : endif
276 26154736 : 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 77759136 : do kk = pver-1,1,-1
287 74768400 : kkm1 = max(1,kk-1)
288 74768400 : kkp1 = min(pver,kk+1)
289 656859136 : do i = il1g,il2g
290 579100000 : mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
291 653868400 : if (mupdudp > mbsth) then
292 :
293 236049804 : conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* &
294 236049804 : 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 74768400 : do k = 3,pver
303 71777664 : km1 = max(1,k-1)
304 630704400 : do i = il1g,il2g
305 627713664 : if (md(i,k) < -mbsth) then
306 :
307 196768404 : cond(i,k) = ( md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) &
308 196768404 : *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 18397672 : do k = ktm,pver
320 15406936 : km1 = max(1,k-1)
321 15406936 : kp1 = min(pver,k+1)
322 186633220 : do i = il1g,il2g
323 168235548 : ii = ideep(i)
324 :
325 : ! version 1 hard to check for roundoff errors
326 168235548 : dcondt(i,k) = &
327 168235548 : +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
328 168235548 : -mu(i,k)* (conu(i,k)-chat(i,k)) &
329 168235548 : +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
330 168235548 : -md(i,k)* (cond(i,k)-chat(i,k)) &
331 520113580 : )/dp(i,k)
332 :
333 : end do
334 : end do
335 :
336 : ! dcont for bottom layer
337 : !
338 6246608 : do k = kbm,pver
339 28967488 : km1 = max(1,k-1)
340 31958224 : do i = il1g,il2g
341 28967488 : if (k == mx(i)) then
342 :
343 : ! version 1
344 23164000 : dcondt(i,k) = (1._kind_phys/dp(i,k))* &
345 23164000 : (-mu(i,k)*(conu(i,k)-chat(i,k)) &
346 23164000 : -md(i,k)*(cond(i,k)-chat(i,k)) &
347 46328000 : )
348 : end if
349 : end do
350 : end do
351 :
352 : ! Initialize to zero everywhere, then scatter tendency back to full array
353 1301387472 : wind_tends(:,:,m) = 0._kind_phys
354 :
355 80749872 : do k = 1,pver
356 683013872 : do i = il1g,il2g
357 602264000 : ii = ideep(i)
358 602264000 : wind_tends(ii,k,m) = dcondt(i,k)
359 : ! Output apparent force on the mean flow from pressure gradient
360 602264000 : pguall(ii,k,m) = -pgu(i,k)
361 602264000 : pgdall(ii,k,m) = -pgd(i,k)
362 602264000 : icwu(ii,k,m) = conu(i,k)
363 680023136 : 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 18397672 : do k = ktm,pver
370 186633220 : do i = il1g,il2g
371 168235548 : ii = ideep(i)
372 168235548 : mflux(i,k,m) = &
373 168235548 : -mu(i,k)* (conu(i,k)-chat(i,k)) &
374 351878032 : -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 18397672 : do k = ktm,pver
382 186633220 : do i = il1g,il2g
383 168235548 : ii = ideep(i)
384 168235548 : km1 = max(1,k-1)
385 168235548 : kp1 = k+1
386 183642484 : 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 9198836 : do k = ktm,pver
399 7703468 : km1 = max(1,k-1)
400 7703468 : kp1 = min(pver,k+1)
401 93316610 : do i = il1g,il2g
402 :
403 84117774 : 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 84117774 : utop = (wind0(i,k,1)+wind0(i,km1,1))/2._kind_phys
408 84117774 : vtop = (wind0(i,k,2)+wind0(i,km1,2))/2._kind_phys
409 84117774 : ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2._kind_phys
410 84117774 : vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2._kind_phys
411 84117774 : fket = utop*mflux(i,k,1) + vtop*mflux(i,k,2) ! top of layer
412 84117774 : 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 84117774 : ketend_cons = (fket-fkeb)/dp(i,k)
416 :
417 : ! tendency in kinetic energy resulting from the momentum transport
418 84117774 : ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))/dt
419 :
420 : ! the difference should be the dissipation
421 84117774 : gset2 = ketend_cons - ketend
422 91821242 : gseten(i,k) = gset2
423 :
424 : end do
425 :
426 : end do
427 :
428 : ! Scatter dry static energy to full array
429 40374936 : do k = 1,pver
430 341506936 : do i = il1g,il2g
431 301132000 : ii = ideep(i)
432 340011568 : seten(ii,k) = gseten(i,k)
433 :
434 : end do
435 : end do
436 :
437 : ! Split out the wind tendencies
438 652189104 : windu_tend(:,:) = wind_tends(:,:,1)
439 652189104 : windv_tend(:,:) = wind_tends(:,:,2)
440 :
441 652189104 : pguallu(:,:) = pguall(:,:,1)
442 652189104 : pguallv(:,:) = pguall(:,:,2)
443 652189104 : pgdallu(:,:) = pgdall(:,:,1)
444 652189104 : pgdallv(:,:) = pgdall(:,:,2)
445 652189104 : icwuu(:ncol,:) = icwu(:,:,1)
446 652189104 : icwuv(:ncol,:) = icwu(:,:,2)
447 652189104 : icwdu(:ncol,:) = icwd(:,:,1)
448 652189104 : icwdv(:ncol,:) = icwd(:,:,2)
449 :
450 1495368 : return
451 : end subroutine zm_conv_momtran_run
452 :
453 :
454 : end module zm_conv_momtran
|
