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 2325123360 : winds(:,:,1) = windu(:,:)
149 2323627992 : winds(:,:,2) = windv(:,:)
150 :
151 : ! Initialize outgoing fields
152 4648751352 : pguall(:,:,:) = 0.0_kind_phys
153 4648751352 : pgdall(:,:,:) = 0.0_kind_phys
154 : ! Initialize in-cloud winds to environmental wind
155 4648751352 : icwu(:ncol,:,:) = winds(:ncol,:,:)
156 4648751352 : icwd(:ncol,:,:) = winds(:ncol,:,:)
157 :
158 : ! Initialize momentum flux and final winds
159 4698689688 : mflux(:,:,:) = 0.0_kind_phys
160 4648751352 : wind0(:,:,:) = 0.0_kind_phys
161 4648751352 : windf(:,:,:) = 0.0_kind_phys
162 :
163 : ! Initialize dry static energy
164 :
165 2323627992 : seten(:,:) = 0.0_kind_phys
166 2323627992 : 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 5720188 : do i = il1g, il2g
175 4224820 : ktm = min(ktm,jt(i))
176 5720188 : 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 281129184 : do k = 1,pver
185 1066945704 : do i =il1g,il2g
186 785816520 : const(i,k) = winds(ideep(i),k,m)
187 1063954968 : 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 281129184 : do k = 1,pver
197 278138448 : km1 = max(1,k-1)
198 1066945704 : do i = il1g, il2g
199 :
200 : ! use arithmetic mean
201 785816520 : chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
202 :
203 : ! Provisional up and down draft values
204 785816520 : conu(i,k) = chat(i,k)
205 785816520 : cond(i,k) = chat(i,k)
206 :
207 : ! provisional tends
208 1063954968 : 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 14431112 : pgu(:il2g,k) = 0.0_kind_phys
222 11440376 : pgd(:il2g,k) = 0.0_kind_phys
223 :
224 275147712 : do k=2,pver-1
225 272156976 : km1 = max(1,k-1)
226 272156976 : kp1 = min(pver,k+1)
227 1044064952 : do i = il1g,il2g
228 :
229 : !interior points
230 :
231 2306751720 : mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
232 3075668960 : + mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
233 :
234 768917240 : pgu(i,k) = - momcu * 0.5_kind_phys * mududp(i,k)
235 :
236 :
237 768917240 : mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
238 768917240 : + md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
239 :
240 1041074216 : 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 11440376 : do i=il1g,il2g
250 :
251 8449640 : mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
252 8449640 : pgu(i,k) = - momcu * mududp(i,k)
253 :
254 8449640 : mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
255 :
256 11440376 : 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 11440376 : k = 2
267 11440376 : km1 = 1
268 11440376 : kk = pver
269 11440376 : kkm1 = max(1,kk-1)
270 11440376 : do i = il1g,il2g
271 8449640 : mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
272 8449640 : 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 11440376 : 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 278138448 : do kk = pver-1,1,-1
287 275147712 : kkm1 = max(1,kk-1)
288 275147712 : kkp1 = min(pver,kk+1)
289 1055505328 : do i = il1g,il2g
290 777366880 : mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
291 1052514592 : if (mupdudp > mbsth) then
292 :
293 270170456 : conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* &
294 270170456 : 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 275147712 : do k = 3,pver
303 272156976 : km1 = max(1,k-1)
304 1044064952 : do i = il1g,il2g
305 1041074216 : if (md(i,k) < -mbsth) then
306 :
307 209102332 : cond(i,k) = ( md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) &
308 209102332 : *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 30634632 : do k = ktm,pver
320 27643896 : km1 = max(1,k-1)
321 27643896 : kp1 = min(pver,k+1)
322 233243994 : do i = il1g,il2g
323 202609362 : ii = ideep(i)
324 :
325 : ! version 1 hard to check for roundoff errors
326 202609362 : dcondt(i,k) = &
327 202609362 : +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
328 202609362 : -mu(i,k)* (conu(i,k)-chat(i,k)) &
329 202609362 : +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
330 202609362 : -md(i,k)* (cond(i,k)-chat(i,k)) &
331 635471982 : )/dp(i,k)
332 :
333 : end do
334 : end do
335 :
336 : ! dcont for bottom layer
337 : !
338 11914548 : do k = kbm,pver
339 62676426 : km1 = max(1,k-1)
340 65667162 : do i = il1g,il2g
341 62676426 : if (k == mx(i)) then
342 :
343 : ! version 1
344 8449640 : dcondt(i,k) = (1._kind_phys/dp(i,k))* &
345 8449640 : (-mu(i,k)*(conu(i,k)-chat(i,k)) &
346 8449640 : -md(i,k)*(cond(i,k)-chat(i,k)) &
347 16899280 : )
348 : end if
349 : end do
350 : end do
351 :
352 : ! Initialize to zero everywhere, then scatter tendency back to full array
353 4647255984 : wind_tends(:,:,m) = 0._kind_phys
354 :
355 281129184 : do k = 1,pver
356 1066945704 : do i = il1g,il2g
357 785816520 : ii = ideep(i)
358 785816520 : wind_tends(ii,k,m) = dcondt(i,k)
359 : ! Output apparent force on the mean flow from pressure gradient
360 785816520 : pguall(ii,k,m) = -pgu(i,k)
361 785816520 : pgdall(ii,k,m) = -pgd(i,k)
362 785816520 : icwu(ii,k,m) = conu(i,k)
363 1063954968 : 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 30634632 : do k = ktm,pver
370 233243994 : do i = il1g,il2g
371 202609362 : ii = ideep(i)
372 202609362 : mflux(i,k,m) = &
373 202609362 : -mu(i,k)* (conu(i,k)-chat(i,k)) &
374 432862620 : -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 30634632 : do k = ktm,pver
382 233243994 : do i = il1g,il2g
383 202609362 : ii = ideep(i)
384 202609362 : km1 = max(1,k-1)
385 202609362 : kp1 = k+1
386 230253258 : 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 15317316 : do k = ktm,pver
399 13821948 : km1 = max(1,k-1)
400 13821948 : kp1 = min(pver,k+1)
401 116621997 : do i = il1g,il2g
402 :
403 101304681 : 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 101304681 : utop = (wind0(i,k,1)+wind0(i,km1,1))/2._kind_phys
408 101304681 : vtop = (wind0(i,k,2)+wind0(i,km1,2))/2._kind_phys
409 101304681 : ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2._kind_phys
410 101304681 : vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2._kind_phys
411 101304681 : fket = utop*mflux(i,k,1) + vtop*mflux(i,k,2) ! top of layer
412 101304681 : 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 101304681 : ketend_cons = (fket-fkeb)/dp(i,k)
416 :
417 : ! tendency in kinetic energy resulting from the momentum transport
418 101304681 : 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 101304681 : gset2 = ketend_cons - ketend
422 115126629 : gseten(i,k) = gset2
423 :
424 : end do
425 :
426 : end do
427 :
428 : ! Scatter dry static energy to full array
429 140564592 : do k = 1,pver
430 533472852 : do i = il1g,il2g
431 392908260 : ii = ideep(i)
432 531977484 : seten(ii,k) = gseten(i,k)
433 :
434 : end do
435 : end do
436 :
437 : ! Split out the wind tendencies
438 2325123360 : windu_tend(:,:) = wind_tends(:,:,1)
439 2325123360 : windv_tend(:,:) = wind_tends(:,:,2)
440 :
441 2325123360 : pguallu(:,:) = pguall(:,:,1)
442 2325123360 : pguallv(:,:) = pguall(:,:,2)
443 2325123360 : pgdallu(:,:) = pgdall(:,:,1)
444 2325123360 : pgdallv(:,:) = pgdall(:,:,2)
445 2325123360 : icwuu(:ncol,:) = icwu(:,:,1)
446 2325123360 : icwuv(:ncol,:) = icwu(:,:,2)
447 2325123360 : icwdu(:ncol,:) = icwd(:,:,1)
448 2325123360 : icwdv(:ncol,:) = icwd(:,:,2)
449 :
450 1495368 : return
451 : end subroutine zm_conv_momtran_run
452 :
453 :
454 : end module zm_conv_momtran
|
