Line data Source code
1 : module zm_conv_convtran
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 : !
10 : ! PUBLIC: interfaces
11 : !
12 : public zm_conv_convtran_run ! convective transport
13 :
14 :
15 : contains
16 :
17 : !===============================================================================
18 : !> \section arg_table_zm_conv_convtran_run Argument Table
19 : !! \htmlinclude zm_conv_convtran_run.html
20 : !!
21 123840 : subroutine zm_conv_convtran_run(ncol, pver, &
22 123840 : doconvtran,q ,ncnst ,mu ,md , &
23 123840 : du ,eu ,ed ,dp ,dsubcld , &
24 247680 : jt ,mx ,ideep ,il1g ,il2g , &
25 371520 : nstep ,fracis ,dqdt ,dpdry ,dt )
26 : !-----------------------------------------------------------------------
27 : !
28 : ! Purpose:
29 : ! Convective transport of trace species
30 : !
31 : ! Mixing ratios may be with respect to either dry or moist air
32 : !
33 : ! Method:
34 : ! <Describe the algorithm(s) used in the routine.>
35 : ! <Also include any applicable external references.>
36 : !
37 : ! Author: P. Rasch
38 : !
39 : !-----------------------------------------------------------------------
40 : ! CACNOTE - replace with CCPP constituents
41 : use constituents, only: cnst_get_type_byind
42 :
43 : implicit none
44 : !-----------------------------------------------------------------------
45 : !
46 : ! Input arguments
47 : !
48 : integer, intent(in) :: ncol
49 : integer, intent(in) :: pver
50 : integer, intent(in) :: ncnst ! number of tracers to transport
51 : logical, intent(in) :: doconvtran(:) ! flag for doing convective transport (ncnst)
52 : real(kind_phys), intent(in) :: q(:,:,:) ! Tracer array including moisture (ncol,pver,ncnst)
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) :: du(:,:) ! Mass detraining from updraft (ncol,pver)
56 : real(kind_phys), intent(in) :: eu(:,:) ! Mass entraining from updraft (ncol,pver)
57 : real(kind_phys), intent(in) :: ed(:,:) ! Mass entraining from downdraft (ncol,pver)
58 : real(kind_phys), intent(in) :: dp(:,:) ! Delta pressure between interfaces (ncol,pver)
59 : real(kind_phys), intent(in) :: dsubcld(:) ! Delta pressure from cloud base to sfc (ncol)
60 : real(kind_phys), intent(in) :: fracis(:,:,:) ! fraction of tracer that is insoluble (ncol,pver,ncnst)
61 :
62 : integer, intent(in) :: jt(:) ! Index of cloud top for each column (ncol)
63 : integer, intent(in) :: mx(:) ! Index of cloud top for each column (ncol)
64 : integer, intent(in) :: ideep(:) ! Gathering array (ncol)
65 : integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
66 : integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
67 : integer, intent(in) :: nstep ! Time step index
68 :
69 : real(kind_phys), intent(in) :: dpdry(:,:) ! Delta pressure between interfaces (ncol,pver)
70 :
71 : real(kind_phys), intent(in) :: dt ! 2 delta t (model time increment)
72 :
73 :
74 : ! input/output
75 :
76 : real(kind_phys), intent(out) :: dqdt(:,:,:) ! Tracer tendency array (ncol,pver,ncnst)
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 km1 ! Work index
86 : integer kp1 ! Work index
87 : integer ktm ! Highest altitude index of cloud top
88 : integer m ! Work index
89 :
90 : real(kind_phys) cabv ! Mix ratio of constituent above
91 : real(kind_phys) cbel ! Mix ratio of constituent below
92 : real(kind_phys) cdifr ! Normalized diff between cabv and cbel
93 247680 : real(kind_phys) chat(ncol,pver) ! Mix ratio in env at interfaces
94 247680 : real(kind_phys) cond(ncol,pver) ! Mix ratio in downdraft at interfaces
95 247680 : real(kind_phys) const(ncol,pver) ! Gathered tracer array
96 247680 : real(kind_phys) fisg(ncol,pver) ! gathered insoluble fraction of tracer
97 247680 : real(kind_phys) conu(ncol,pver) ! Mix ratio in updraft at interfaces
98 247680 : real(kind_phys) dcondt(ncol,pver) ! Gathered tend array
99 : real(kind_phys) small ! A small number
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 247680 : real(kind_phys) dutmp(ncol,pver) ! Mass detraining from updraft
109 247680 : real(kind_phys) eutmp(ncol,pver) ! Mass entraining from updraft
110 247680 : real(kind_phys) edtmp(ncol,pver) ! Mass entraining from downdraft
111 247680 : real(kind_phys) dptmp(ncol,pver) ! Delta pressure between interfaces
112 : real(kind_phys) total(ncol)
113 : real(kind_phys) negadt,qtmp
114 :
115 : !-----------------------------------------------------------------------
116 : !
117 123840 : small = 1.e-36_kind_phys
118 : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
119 123840 : mbsth = 1.e-15_kind_phys
120 :
121 : ! Find the highest level top and bottom levels of convection
122 123840 : ktm = pver
123 123840 : kbm = pver
124 470817 : do i = il1g, il2g
125 346977 : ktm = min(ktm,jt(i))
126 470817 : kbm = min(kbm,mx(i))
127 : end do
128 :
129 : ! Loop ever each constituent
130 5077440 : do m = 2, ncnst
131 5077440 : if (doconvtran(m)) then
132 :
133 1297224 : if (cnst_get_type_byind(m).eq.'dry') then
134 60824016 : do k = 1,pver
135 229490118 : do i =il1g,il2g
136 168666102 : dptmp(i,k) = dpdry(i,k)
137 168666102 : dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
138 168666102 : eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
139 228843054 : edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
140 : end do
141 : end do
142 : else
143 62412264 : do k = 1,pver
144 230470830 : do i =il1g,il2g
145 169355790 : dptmp(i,k) = dp(i,k)
146 169355790 : dutmp(i,k) = du(i,k)
147 169355790 : eutmp(i,k) = eu(i,k)
148 229820670 : edtmp(i,k) = ed(i,k)
149 : end do
150 : end do
151 : endif
152 : ! dptmp = dp
153 :
154 : ! Gather up the constituent and set tend to zero
155 121939056 : do k = 1,pver
156 459960948 : do i =il1g,il2g
157 338021892 : const(i,k) = q(ideep(i),k,m)
158 458663724 : fisg(i,k) = fracis(ideep(i),k,m)
159 : end do
160 : end do
161 :
162 : ! From now on work only with gathered data
163 :
164 : ! Interpolate environment tracer values to interfaces
165 121939056 : do k = 1,pver
166 120641832 : km1 = max(1,k-1)
167 459960948 : do i = il1g, il2g
168 338021892 : minc = min(const(i,km1),const(i,k))
169 338021892 : maxc = max(const(i,km1),const(i,k))
170 338021892 : if (minc < 0) then
171 : cdifr = 0._kind_phys
172 : else
173 338021892 : cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small)
174 : endif
175 :
176 : ! If the two layers differ significantly use a geometric averaging
177 : ! procedure
178 338021892 : if (cdifr > 1.E-6_kind_phys) then
179 244939541 : cabv = max(const(i,km1),maxc*1.e-12_kind_phys)
180 244939541 : cbel = max(const(i,k),maxc*1.e-12_kind_phys)
181 244939541 : chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
182 :
183 : else ! Small diff, so just arithmetic mean
184 93082351 : chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
185 : end if
186 :
187 : ! Provisional up and down draft values
188 338021892 : conu(i,k) = chat(i,k)
189 338021892 : cond(i,k) = chat(i,k)
190 :
191 : ! provisional tends
192 458663724 : dcondt(i,k) = 0._kind_phys
193 :
194 : end do
195 : end do
196 :
197 : ! Do levels adjacent to top and bottom
198 4931868 : k = 2
199 4931868 : km1 = 1
200 4931868 : kk = pver
201 4931868 : do i = il1g,il2g
202 3634644 : mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
203 3634644 : if (mupdudp > mbsth) then
204 0 : conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp
205 : endif
206 4931868 : if (md(i,k) < -mbsth) then
207 0 : cond(i,k) = (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k)
208 : endif
209 : end do
210 :
211 : ! Updraft from bottom to top
212 120641832 : do kk = pver-1,1,-1
213 119344608 : kkp1 = min(pver,kk+1)
214 455029080 : do i = il1g,il2g
215 334387248 : mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
216 453731856 : if (mupdudp > mbsth) then
217 55763550 : conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
218 111527100 : const(i,kk)*dptmp(i,kk) )/mupdudp
219 : endif
220 : end do
221 : end do
222 :
223 : ! Downdraft from top to bottom
224 119344608 : do k = 3,pver
225 118047384 : km1 = max(1,k-1)
226 450097212 : do i = il1g,il2g
227 448799988 : if (md(i,k) < -mbsth) then
228 85525742 : cond(i,k) = ( md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
229 85525742 : *dptmp(i,km1) )/md(i,k)
230 : endif
231 : end do
232 : end do
233 :
234 :
235 12905354 : do k = ktm,pver
236 11608130 : km1 = max(1,k-1)
237 11608130 : kp1 = min(pver,k+1)
238 97812107 : do i = il1g,il2g
239 :
240 : ! limit fluxes outside convection to mass in appropriate layer
241 : ! these limiters are probably only safe for positive definite quantitities
242 : ! it assumes that mu and md already satify a courant number limit of 1
243 339627012 : fluxin = mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
244 424533765 : -(md(i,k) *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1)))
245 : fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) &
246 84906753 : -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
247 :
248 84906753 : netflux = fluxin - fluxout
249 84906753 : if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
250 4735199 : netflux = 0._kind_phys
251 : endif
252 96514883 : dcondt(i,k) = netflux/dptmp(i,k)
253 : end do
254 : end do
255 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256 : !
257 5087271 : do k = kbm,pver
258 3790047 : km1 = max(1,k-1)
259 27743084 : do i = il1g,il2g
260 26445860 : if (k == mx(i)) then
261 :
262 3634644 : fluxin = mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
263 3634644 : fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
264 :
265 3634644 : netflux = fluxin - fluxout
266 3634644 : if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
267 0 : netflux = 0._kind_phys
268 : endif
269 3634644 : dcondt(i,k) = netflux/dptmp(i,k)
270 19021169 : else if (k > mx(i)) then
271 16242147 : dcondt(i,k) = 0._kind_phys
272 : end if
273 : end do
274 : end do
275 :
276 : ! Initialize to zero everywhere, then scatter tendency back to full array
277 2015735256 : dqdt(:,:,m) = 0._kind_phys
278 121939056 : do k = 1,pver
279 120641832 : kp1 = min(pver,k+1)
280 459960948 : do i = il1g,il2g
281 458663724 : dqdt(ideep(i),k,m) = dcondt(i,k)
282 : end do
283 : end do
284 :
285 : end if ! for doconvtran
286 :
287 : end do
288 :
289 123840 : return
290 : end subroutine zm_conv_convtran_run
291 :
292 :
293 : end module zm_conv_convtran
|
