← Back to Documentation Home
CCPP SciDoc v8.0.0  v8.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
progomega_calc.f90
1
2
7
8
9module progomega
10
11 implicit none
12
13 public progomega_calc
14
15contains
16
22 subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, &
23 zi,cnvflg,omegaout,grav,buo,drag,wush,lbb1,lbb2,lbb3,dt_decay)
24
25 use machine, only : kind_phys
26
27 implicit none
28
29 integer, intent(in) :: im,km
30 integer, intent(in) :: kbcon1(im),ktcon(im)
31
32 real(kind=kind_phys), intent(in) :: delt,grav
33 real(kind=kind_phys), intent(in) :: lbb1,lbb2,lbb3,dt_decay
34 real(kind=kind_phys), intent(in) :: omegain(im,km)
35 real(kind=kind_phys), intent(in) :: del(im,km),zi(im,km)
36 real(kind=kind_phys), intent(in) :: drag(im,km)
37 real(kind=kind_phys), intent(in) :: buo(im,km)
38 real(kind=kind_phys), intent(in) :: wush(im,km)
39
40 real(kind=kind_phys), intent(inout) :: omegaout(im,km)
41
42 logical, intent(in) :: cnvflg(im)
43 logical, intent(in) :: first_time_step
44 logical, intent(in) :: flag_restart
45
46 !--------------------------------------------------------------------
47 ! Local arrays
48 !
49 ! omega = state at beginning of current internal substep
50 ! omega_new = state at end of current internal substep
51 !--------------------------------------------------------------------
52
53 real(kind=kind_phys) :: omega(im,km)
54 real(kind=kind_phys) :: omega_new(im,km)
55 real(kind=kind_phys) :: omega_start(im,km)
56
57 real(kind=kind_phys) :: terma(im,km)
58 real(kind=kind_phys) :: termb(im,km)
59 real(kind=kind_phys) :: termc(im,km)
60 real(kind=kind_phys) :: memory_term,buoy_term,adv_term
61 real(kind=kind_phys) :: decay_fac
62 !--------------------------------------------------------------------
63 ! Scalars
64 !--------------------------------------------------------------------
65
66 real(kind=kind_phys) :: dp,dz,pi_conv,discr
67 real(kind=kind_phys) :: rhs_exp
68
69 real(kind=kind_phys) :: omega_eps
70 real(kind=kind_phys) :: a_eps,b_eps,disc_eps
71
72 real(kind=kind_phys) :: cfl_target
73 real(kind=kind_phys) :: dt_sub
74 real(kind=kind_phys) :: dt_remaining
75 real(kind=kind_phys) :: dt_cfl
76 real(kind=kind_phys) :: time_done
77
78 integer :: i,k
79
80 logical :: active_point_found
81
82 !--------------------------------------------------------------------
83 ! Numerical parameters
84 !--------------------------------------------------------------------
85
86 ! Target CFL for explicit pressure-coordinate vertical advection.
87 cfl_target = 0.8_kind_phys
88
89 omega_eps = 1.0e-5_kind_phys
90 a_eps = 1.0e-12_kind_phys
91 b_eps = 1.0e-12_kind_phys
92 disc_eps = 1.0e-12_kind_phys
93
94 ! Decay applied over one host physics timestep.
95 decay_fac = exp(-delt/dt_decay)
96
97 !--------------------------------------------------------------------
98 ! Initialize from incoming prognostic tracer
99 !--------------------------------------------------------------------
100
101 do k = 1,km
102 do i = 1,im
103
104 terma(i,k) = 0.0_kind_phys
105 termb(i,k) = 0.0_kind_phys
106 termc(i,k) = 0.0_kind_phys
107
108 omega(i,k) = omegain(i,k)
109 omega_new(i,k) = omegain(i,k)
110 omegaout(i,k) = omegain(i,k)
111
112 enddo
113 enddo
114
115 !--------------------------------------------------------------------
116 ! Retain memory when convection is inactive.
117 !--------------------------------------------------------------------
118
119 do k = 1,km
120 do i = 1,im
121
122 if (.not. cnvflg(i)) then
123
124 omega(i,k) = omegain(i,k)
125 omega_new(i,k) = omega(i,k)
126 omegaout(i,k) = omega(i,k)
127
128 endif
129
130 ! Remove numerically negligible values.
131 if (abs(omega(i,k)) < omega_eps) then
132
133 omega(i,k) = 0.0_kind_phys
134 omega_new(i,k) = 0.0_kind_phys
135 omegaout(i,k) = 0.0_kind_phys
136
137 endif
138
139 enddo
140 enddo
141
142 !--------------------------------------------------------------------
143 ! Decay prognostic updraft memory when convection is inactive.
144 !
145 ! The decay is applied once over the full host physics timestep.
146 !--------------------------------------------------------------------
147
148 do k = 1,km
149 do i = 1,im
150
151 if (.not. cnvflg(i)) then
152
153 omega(i,k) = omegain(i,k) * decay_fac
154 omega_new(i,k) = omega(i,k)
155 omegaout(i,k) = omega(i,k)
156
157 endif
158
159 ! Remove numerically negligible values.
160 if (abs(omega(i,k)) < omega_eps) then
161
162 omega(i,k) = 0.0_kind_phys
163 omega_new(i,k) = 0.0_kind_phys
164 omegaout(i,k) = 0.0_kind_phys
165
166 endif
167
168 enddo
169 enddo
170
171 !--------------------------------------------------------------------
172 ! Cold-start initialization
173 !--------------------------------------------------------------------
174
175 if (first_time_step .and. .not. flag_restart) then
176
177 do k = 1,km
178 do i = 1,im
179
180 if (cnvflg(i)) then
181
182 if (k >= kbcon1(i) .and. k < ktcon(i)) then
183
184 omega(i,k) = -1.2_kind_phys
185 omega_new(i,k) = -1.2_kind_phys
186 omegaout(i,k) = -1.2_kind_phys
187
188 endif
189
190 endif
191
192 enddo
193 enddo
194
195 endif
196
197 !--------------------------------------------------------------------
198 ! Adaptive CFL-controlled subcycling
199 !--------------------------------------------------------------------
200
201 time_done = 0.0_kind_phys
202
203 ! Save the state entering the prognostic integration so that the
204 ! host-timestep local tendency can be returned after all substeps.
205 omega_start(:,:) = omega(:,:)
206
207 do while (time_done < delt)
208
209 dt_remaining = delt - time_done
210
211 !---------------------------------------------------------------
212 ! Determine CFL-limited timestep from current omega profile.
213 !
214 ! The explicit vertical-advection term is zero at cloud base,
215 ! so only levels above kbcon1 are included in the CFL estimate.
216 !---------------------------------------------------------------
217
218 dt_cfl = dt_remaining
219 active_point_found = .false.
220
221 do k = 2,km
222 do i = 1,im
223
224 if (cnvflg(i)) then
225 if (k > kbcon1(i) .and. k < ktcon(i)) then
226 dp = 1000.0_kind_phys * del(i,k)
227 if (dp > 0.0_kind_phys) then
228 active_point_found = .true.
229 if (abs(omega(i,k)) > omega_eps) then
230 dt_cfl = min(dt_cfl, &
231 cfl_target * dp / abs(omega(i,k)))
232 endif
233 endif
234 endif
235 endif
236 enddo
237 enddo
238
239 ! If there are no active plume points above cloud base,
240 ! no prognostic subcycling is needed.
241 if (.not. active_point_found) exit
242
243 dt_sub = min(dt_remaining,dt_cfl)
244
245 ! Numerical protection against pathological tiny timesteps.
246 if (dt_sub <= 1.0e-8_kind_phys) then
247 dt_sub = dt_remaining
248 endif
249
250 !---------------------------------------------------------------
251 ! Start new substep from previous-substep profile.
252 !
253 ! This ensures every vertical level uses the same
254 ! profile for the explicit vertical-advection term.
255 !---------------------------------------------------------------
256
257 omega_new(:,:) = omega(:,:)
258
259 !---------------------------------------------------------------
260 ! Solve prognostic momentum equation following
261 ! Bengtsson et al. 2026
262 !---------------------------------------------------------------
263
264 do k = 2,km
265 do i = 1,im
266
267 if (cnvflg(i)) then
268 if (k >= kbcon1(i) .and. k < ktcon(i)) then
269 ! Cloud-base boundary condition.
270 omega_new(i,kbcon1(i)) = 0.0_kind_phys
271
272 dp = 1000.0_kind_phys * del(i,k)
273 dz = zi(i,k+1) - zi(i,k)
274
275 if (dp <= 0.0_kind_phys) cycle
276 if (abs(dz) <= 1.0e-12_kind_phys) cycle
277
278 pi_conv = dp/dz
279
280 !----------------------------------------------------
281 ! Explicit contributions
282 !----------------------------------------------------
283
284 memory_term = omega(i,k)
285
286 buoy_term = -0.5_kind_phys * dt_sub * lbb2 &
287 * buo(i,k) * pi_conv
288
289 if (k == kbcon1(i)) then
290 adv_term = 0.0_kind_phys
291 else
292 adv_term = -dt_sub * omega(i,k) &
293 * (omega(i,k-1)-omega(i,k)) / dp
294
295 endif
296 rhs_exp = memory_term + buoy_term + adv_term
297
298 !----------------------------------------------------
299 ! Quadratic coefficients
300 !
301 ! A * omega_new**2 + B * omega_new + C = 0
302 !----------------------------------------------------
303
304 terma(i,k) = -0.5_kind_phys * dt_sub * lbb1 &
305 * drag(i,k) / pi_conv
306
307 termb(i,k) = 1.0_kind_phys &
308 + 0.5_kind_phys * dt_sub * wush(i,k)
309
310 termc(i,k) = -rhs_exp
311
312 !----------------------------------------------------
313 ! Robust quadratic / linear solution
314 !----------------------------------------------------
315
316 if (abs(terma(i,k)) < a_eps) then
317
318 if (abs(termb(i,k)) > b_eps) then
319 omega_new(i,k) = -termc(i,k) / termb(i,k)
320 else
321 ! Degenerate case: retain previous state.
322 omega_new(i,k) = omega(i,k)
323 endif
324
325 else
326
327 discr = termb(i,k)**2 &
328 - 4.0_kind_phys * terma(i,k) * termc(i,k)
329
330 if (discr >= -disc_eps) then
331
332 discr = max(discr,0.0_kind_phys)
333
334 omega_new(i,k) = &
335 (-termb(i,k) + sqrt(discr)) &
336 / (2.0_kind_phys * terma(i,k))
337
338 else
339
340 ! No real solution: retain previous state.
341 omega_new(i,k) = omega(i,k)
342
343 endif
344
345 endif
346
347 !----------------------------------------------------
348 ! Physical bounds for active updrafts
349 !----------------------------------------------------
350
351 omega_new(i,k) = max( &
352 min(omega_new(i,k),-1.2_kind_phys), &
353 -80.0_kind_phys)
354
355 endif
356
357 endif
358
359 enddo
360 enddo
361
362 !---------------------------------------------------------------
363 ! Advance entire vertical profile simultaneously
364 !---------------------------------------------------------------
365
366 omega(:,:) = omega_new(:,:)
367
368 time_done = time_done + dt_sub
369
370 ! Prevent tiny floating-point remainder from generating an
371 ! unnecessary additional substep.
372 if (delt-time_done < 1.0e-8_kind_phys*delt) then
373 time_done = delt
374 endif
375
376 enddo
377
378 !--------------------------------------------------------------------
379 ! Return final prognostic state and host-timestep local tendency
380 !--------------------------------------------------------------------
381
382 omegaout(:,:) = omega(:,:)
383
384 end subroutine progomega_calc
385
386end module progomega
progomega
This module contains the subroutine that calculates the prognostic updraft velocity that is used for ...
Definition progomega_calc.f90:9