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 1495368 : subroutine zm_conv_convtran_run(ncol, pver, &
22 1495368 : doconvtran,q ,ncnst ,mu ,md , &
23 1495368 : du ,eu ,ed ,dp ,dsubcld , &
24 2990736 : jt ,mx ,ideep ,il1g ,il2g , &
25 4486104 : nstep ,fracis ,dqdt ,dpdry)
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 : ! input/output
72 :
73 : real(kind_phys), intent(out) :: dqdt(:,:,:) ! Tracer tendency array (ncol,pver,ncnst)
74 :
75 : !--------------------------Local Variables------------------------------
76 :
77 : integer i ! Work index
78 : integer k ! Work index
79 : integer kbm ! Highest altitude index of cloud base
80 : integer kk ! Work index
81 : integer kkp1 ! Work index
82 : integer km1 ! Work index
83 : integer kp1 ! Work index
84 : integer ktm ! Highest altitude index of cloud top
85 : integer m ! Work index
86 :
87 : real(kind_phys) cabv ! Mix ratio of constituent above
88 : real(kind_phys) cbel ! Mix ratio of constituent below
89 : real(kind_phys) cdifr ! Normalized diff between cabv and cbel
90 2990736 : real(kind_phys) chat(ncol,pver) ! Mix ratio in env at interfaces
91 2990736 : real(kind_phys) cond(ncol,pver) ! Mix ratio in downdraft at interfaces
92 2990736 : real(kind_phys) const(ncol,pver) ! Gathered tracer array
93 2990736 : real(kind_phys) fisg(ncol,pver) ! gathered insoluble fraction of tracer
94 2990736 : real(kind_phys) conu(ncol,pver) ! Mix ratio in updraft at interfaces
95 2990736 : real(kind_phys) dcondt(ncol,pver) ! Gathered tend array
96 : real(kind_phys) small ! A small number
97 : real(kind_phys) mbsth ! Threshold for mass fluxes
98 : real(kind_phys) mupdudp ! A work variable
99 : real(kind_phys) minc ! A work variable
100 : real(kind_phys) maxc ! A work variable
101 : real(kind_phys) fluxin ! A work variable
102 : real(kind_phys) fluxout ! A work variable
103 : real(kind_phys) netflux ! A work variable
104 :
105 2990736 : real(kind_phys) dutmp(ncol,pver) ! Mass detraining from updraft
106 2990736 : real(kind_phys) eutmp(ncol,pver) ! Mass entraining from updraft
107 2990736 : real(kind_phys) edtmp(ncol,pver) ! Mass entraining from downdraft
108 2990736 : real(kind_phys) dptmp(ncol,pver) ! Delta pressure between interfaces
109 : real(kind_phys) total(ncol)
110 : real(kind_phys) negadt,qtmp
111 :
112 : !-----------------------------------------------------------------------
113 : !
114 1495368 : small = 1.e-36_kind_phys
115 : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
116 1495368 : mbsth = 1.e-15_kind_phys
117 :
118 : ! Find the highest level top and bottom levels of convection
119 1495368 : ktm = pver
120 1495368 : kbm = pver
121 5720188 : do i = il1g, il2g
122 4224820 : ktm = min(ktm,jt(i))
123 5720188 : kbm = min(kbm,mx(i))
124 : end do
125 :
126 : ! Loop ever each constituent
127 16449048 : do m = 2, ncnst
128 16449048 : if (doconvtran(m)) then
129 :
130 14953680 : if (cnst_get_type_byind(m).eq.'dry') then
131 0 : do k = 1,pver
132 0 : do i =il1g,il2g
133 0 : dptmp(i,k) = dpdry(i,k)
134 0 : dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
135 0 : eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
136 0 : edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
137 : end do
138 : end do
139 : else
140 1420599600 : do k = 1,pver
141 5334728520 : do i =il1g,il2g
142 3929082600 : dptmp(i,k) = dp(i,k)
143 3929082600 : dutmp(i,k) = du(i,k)
144 3929082600 : eutmp(i,k) = eu(i,k)
145 5319774840 : edtmp(i,k) = ed(i,k)
146 : end do
147 : end do
148 : endif
149 : ! dptmp = dp
150 :
151 : ! Gather up the constituent and set tend to zero
152 1405645920 : do k = 1,pver
153 5334728520 : do i =il1g,il2g
154 3929082600 : const(i,k) = q(ideep(i),k,m)
155 5319774840 : fisg(i,k) = fracis(ideep(i),k,m)
156 : end do
157 : end do
158 :
159 : ! From now on work only with gathered data
160 :
161 : ! Interpolate environment tracer values to interfaces
162 1405645920 : do k = 1,pver
163 1390692240 : km1 = max(1,k-1)
164 5334728520 : do i = il1g, il2g
165 3929082600 : minc = min(const(i,km1),const(i,k))
166 3929082600 : maxc = max(const(i,km1),const(i,k))
167 3929082600 : if (minc < 0) then
168 : cdifr = 0._kind_phys
169 : else
170 3929082600 : cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small)
171 : endif
172 :
173 : ! If the two layers differ significantly use a geometric averaging
174 : ! procedure
175 3929082600 : if (cdifr > 1.E-6_kind_phys) then
176 2011955509 : cabv = max(const(i,km1),maxc*1.e-12_kind_phys)
177 2011955509 : cbel = max(const(i,k),maxc*1.e-12_kind_phys)
178 2011955509 : chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
179 :
180 : else ! Small diff, so just arithmetic mean
181 1917127091 : chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
182 : end if
183 :
184 : ! Provisional up and down draft values
185 3929082600 : conu(i,k) = chat(i,k)
186 3929082600 : cond(i,k) = chat(i,k)
187 :
188 : ! provisional tends
189 5319774840 : dcondt(i,k) = 0._kind_phys
190 :
191 : end do
192 : end do
193 :
194 : ! Do levels adjacent to top and bottom
195 57201880 : k = 2
196 57201880 : km1 = 1
197 57201880 : kk = pver
198 57201880 : do i = il1g,il2g
199 42248200 : mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
200 42248200 : if (mupdudp > mbsth) then
201 0 : conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp
202 : endif
203 57201880 : if (md(i,k) < -mbsth) then
204 0 : cond(i,k) = (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k)
205 : endif
206 : end do
207 :
208 : ! Updraft from bottom to top
209 1390692240 : do kk = pver-1,1,-1
210 1375738560 : kkp1 = min(pver,kk+1)
211 5277526640 : do i = il1g,il2g
212 3886834400 : mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
213 5262572960 : if (mupdudp > mbsth) then
214 675426140 : conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
215 1350852280 : const(i,kk)*dptmp(i,kk) )/mupdudp
216 : endif
217 : end do
218 : end do
219 :
220 : ! Downdraft from top to bottom
221 1375738560 : do k = 3,pver
222 1360784880 : km1 = max(1,k-1)
223 5220324760 : do i = il1g,il2g
224 5205371080 : if (md(i,k) < -mbsth) then
225 1045511660 : cond(i,k) = ( md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
226 1045511660 : *dptmp(i,km1) )/md(i,k)
227 : endif
228 : end do
229 : end do
230 :
231 :
232 153173160 : do k = ktm,pver
233 138219480 : km1 = max(1,k-1)
234 138219480 : kp1 = min(pver,k+1)
235 1166219970 : do i = il1g,il2g
236 :
237 : ! limit fluxes outside convection to mass in appropriate layer
238 : ! these limiters are probably only safe for positive definite quantitities
239 : ! it assumes that mu and md already satify a courant number limit of 1
240 4052187240 : fluxin = mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
241 5065234050 : -(md(i,k) *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1)))
242 : fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) &
243 1013046810 : -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
244 :
245 1013046810 : netflux = fluxin - fluxout
246 1013046810 : if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
247 146297763 : netflux = 0._kind_phys
248 : endif
249 1151266290 : dcondt(i,k) = netflux/dptmp(i,k)
250 : end do
251 : end do
252 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
253 : !
254 59572740 : do k = kbm,pver
255 44619060 : km1 = max(1,k-1)
256 328335810 : do i = il1g,il2g
257 313382130 : if (k == mx(i)) then
258 :
259 42248200 : fluxin = mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
260 42248200 : fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
261 :
262 42248200 : netflux = fluxin - fluxout
263 42248200 : if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
264 0 : netflux = 0._kind_phys
265 : endif
266 42248200 : dcondt(i,k) = netflux/dptmp(i,k)
267 226514870 : else if (k > mx(i)) then
268 189933590 : dcondt(i,k) = 0._kind_phys
269 : end if
270 : end do
271 : end do
272 :
273 : ! Initialize to zero everywhere, then scatter tendency back to full array
274 23236279920 : dqdt(:,:,m) = 0._kind_phys
275 1405645920 : do k = 1,pver
276 1390692240 : kp1 = min(pver,k+1)
277 5334728520 : do i = il1g,il2g
278 5319774840 : dqdt(ideep(i),k,m) = dcondt(i,k)
279 : end do
280 : end do
281 :
282 : end if ! for doconvtran
283 :
284 : end do
285 :
286 1495368 : return
287 : end subroutine zm_conv_convtran_run
288 :
289 :
290 : end module zm_conv_convtran
|
