Line data Source code
1 : module pumas_stochastic_collect_tau
2 : ! From Morrison (Lebo, originally TAU bin code)
3 : ! Gettelman and Chen 2018
4 : !the subroutines take in air density, air temperature, and the bin mass boundaries, and
5 : !output the mass and number mixing ratio tendencies in each bin directly.
6 : !this is then wrapped for CAM.
7 :
8 : ! note, this is now coded locally. Want the CAM interface to be over i,k I think.
9 :
10 : #ifndef HAVE_GAMMA_INTRINSICS
11 : use shr_spfn_mod, only: gamma => shr_spfn_gamma
12 : #endif
13 :
14 : use shr_kind_mod, only: r8=>shr_kind_r8
15 : use cam_history, only: addfld
16 : use micro_pumas_utils, only: pi, rhow, qsmall, VLENS
17 : use cam_logfile, only: iulog
18 :
19 : implicit none
20 : private
21 : save
22 :
23 : ! Subroutines
24 : public :: pumas_stochastic_kernel_init, pumas_stochastic_collect_tau_tend
25 :
26 : !In the module top, declare the following so that these can be used throughout the module:
27 :
28 : integer, parameter, public :: ncd = 35
29 : integer, parameter, public :: ncdp = ncd + 1
30 : integer, parameter, public :: ncdl = ncd
31 : integer, parameter, public :: ncdpl = ncdl+1
32 :
33 : ! for Zach's collision-coalescence code
34 :
35 : real(r8), private :: knn(ncd,ncd)
36 :
37 : real(r8), private :: mmean(ncd), diammean(ncd) ! kg & m at bin mid-points
38 : real(r8), private :: medge(ncdp), diamedge(ncdp) ! kg & m at bin edges
39 : integer, private :: cutoff_id ! cutoff between cloud water and rain drop, D = 40 microns
40 :
41 : ! Assume 6 microns for each...
42 : real(r8), parameter :: m1 = 4._r8/3._r8*pi*rhow*(6.e-6_r8)**3
43 :
44 : !$acc declare create(knn,cutoff_id,mmean,diammean,medge,diamedge)
45 :
46 : !===============================================================================
47 : contains
48 : !===============================================================================
49 :
50 0 : subroutine calc_bins
51 :
52 : real(r8) :: DIAM(ncdp)
53 : real(r8) :: X(ncdp)
54 : real(r8) :: radsl(ncdp)
55 : real(r8) :: radl(ncd)
56 : integer :: L, lcl
57 : real(r8) :: kkfac
58 :
59 : !Then before doing any calculations you'll need to calculate the bin mass grid
60 : ! (note this code could be cleaned up, I'm just taking it as it's used in our bin scheme).
61 : ! This only needs to be done once, since we'll use the same bin mass grid for all calculations.
62 :
63 : ! use mass doubling bins from Graham Feingold (note cgs units)
64 :
65 0 : DIAM(1)=1.5625*2.E-04_r8 ! cm
66 0 : X(1)=PI/6._r8*DIAM(1)**3*rhow/1000._r8 ! rhow kg/m3 --> g/cm3
67 0 : radsl(1) = X(1) ! grams
68 :
69 0 : DO l=2,ncdp
70 0 : X(l)=2._r8*X(l-1)
71 0 : DIAM(l)=(6._r8/pi*X(l)*1000._r8/rhow)**(1._r8/3._r8) ! cm
72 0 : radsl(l)=X(l)
73 : ENDDO
74 :
75 : ! now get bin mid-points
76 :
77 0 : do l=1,ncd
78 0 : radl(l)=(radsl(l)+radsl(l+1))/2._r8 ! grams
79 0 : diammean(l) = (6._r8/pi*radl(l)*1000._r8/rhow)**(1._r8/3._r8) ! cm
80 : end do
81 :
82 : ! set bin grid for method of moments
83 :
84 : ! for method of moments
85 :
86 0 : do lcl = 1,ncd+1
87 0 : medge(lcl) = radsl(lcl) ! grams
88 0 : diamedge(lcl) = DIAM(lcl) ! cm
89 : enddo
90 :
91 0 : do lcl = 1,ncd
92 0 : mmean(lcl) = radl(lcl)
93 0 : diammean(lcl) = diammean(lcl)
94 : enddo
95 :
96 0 : do lcl = ncdp,1,-1
97 0 : if( diamedge(lcl).ge.40.e-4_r8 ) cutoff_id = lcl
98 : end do
99 :
100 0 : end subroutine calc_bins
101 :
102 0 : subroutine pumas_stochastic_kernel_init(kernel_filename)
103 :
104 : use cam_history_support, only: add_hist_coord
105 :
106 : character(len=*), intent(in) :: kernel_filename ! Full pathname to kernel file
107 :
108 : integer :: iunit ! unit number of opened file for collection kernel code from a lookup table.
109 :
110 : integer :: idd, jdd
111 : real(r8) :: kkfac
112 :
113 0 : call calc_bins
114 :
115 : ! Read in the collection kernel code from a lookup table. Again, this only needs to be done once.
116 : ! use kernel from Zach (who got it from Jerry)
117 :
118 0 : KNN(:,:)=0._r8 ! initialize values
119 0 : kkfac=1.5_r8 ! from Zach
120 :
121 0 : open(newunit=iunit,file=kernel_filename,status='old')
122 :
123 0 : do idd=1,ncd
124 0 : do jdd=1,idd
125 0 : READ(iunit,941) KNN(IDD,JDD)
126 : 941 FORMAT(2X,E12.5)
127 :
128 0 : KNN(IDD,JDD)=(mmean(IDD)*kkfac+mmean(JDD)*kkfac)*KNN(IDD,JDD)
129 0 : if (knn(idd,jdd) < 0._r8) knn(idd,jdd)=0._r8
130 : end do
131 : end do
132 :
133 : !$acc update device(knn,cutoff_id,mmean,diammean,medge,diamedge)
134 :
135 0 : end subroutine pumas_stochastic_kernel_init
136 :
137 : !main driver routine
138 : !needs to pull in i,k fields (so might need dimensions here too)
139 :
140 0 : subroutine pumas_stochastic_collect_tau_tend(deltatin,t,rho,qc,qr,qcin, &
141 0 : ncin,qrin,nrin,lcldm,precip_frac,mu_c,lambda_c, &
142 0 : n0r,lambda_r,qcin_new,ncin_new,qrin_new,nrin_new, &
143 0 : qctend,nctend,qrtend,nrtend,qctend_TAU,nctend_TAU, &
144 0 : qrtend_TAU,nrtend_TAU,scale_qc,scale_nc,scale_qr, &
145 0 : scale_nr,amk_c,ank_c,amk_r,ank_r,amk,ank,amk_out, &
146 0 : ank_out,gmnnn_lmnnn_TAU,mgncol,nlev)
147 :
148 : !inputs: mgncol,nlev,t,rho,qcin,ncin,qrin,nrin
149 : !outputs: qctend,nctend,qrtend,nrtend
150 : !not sure if we want to output bins (extra dimension). Good for testing?
151 :
152 : integer, intent(in) :: mgncol,nlev
153 : real(r8), intent(in) :: deltatin
154 : real(r8), intent(in) :: t(mgncol,nlev)
155 : real(r8), intent(in) :: rho(mgncol,nlev)
156 : real(r8), intent(in) :: qc(mgncol,nlev)
157 : real(r8), intent(in) :: qr(mgncol,nlev)
158 : real(r8), intent(in) :: qcin(mgncol,nlev)
159 : real(r8), intent(in) :: ncin(mgncol,nlev)
160 : real(r8), intent(in) :: qrin(mgncol,nlev)
161 : real(r8), intent(in) :: nrin(mgncol,nlev)
162 : real(r8), intent(in) :: lcldm(mgncol,nlev)
163 : real(r8), intent(in) :: precip_frac(mgncol,nlev)
164 : real(r8), intent(inout) :: qctend(mgncol,nlev)
165 : real(r8), intent(inout) :: nctend(mgncol,nlev)
166 : real(r8), intent(inout) :: qrtend(mgncol,nlev)
167 : real(r8), intent(inout) :: nrtend(mgncol,nlev)
168 : real(r8), intent(out) :: qctend_TAU(mgncol,nlev)
169 : real(r8), intent(out) :: nctend_TAU(mgncol,nlev)
170 : real(r8), intent(out) :: qrtend_TAU(mgncol,nlev)
171 : real(r8), intent(out) :: nrtend_TAU(mgncol,nlev)
172 :
173 : real(r8), intent(out) :: scale_qc(mgncol,nlev)
174 : real(r8), intent(out) :: scale_nc(mgncol,nlev)
175 : real(r8), intent(out) :: scale_qr(mgncol,nlev)
176 : real(r8), intent(out) :: scale_nr(mgncol,nlev)
177 :
178 : real(r8), intent(out) :: amk_c(mgncol,nlev,ncd)
179 : real(r8), intent(out) :: ank_c(mgncol,nlev,ncd)
180 : real(r8), intent(out) :: amk_r(mgncol,nlev,ncd)
181 : real(r8), intent(out) :: ank_r(mgncol,nlev,ncd)
182 : real(r8), intent(out) :: amk(mgncol,nlev,ncd)
183 : real(r8), intent(out) :: ank(mgncol,nlev,ncd)
184 : real(r8), intent(out) :: amk_out(mgncol,nlev,ncd)
185 : real(r8), intent(out) :: ank_out(mgncol,nlev,ncd)
186 :
187 : real(r8), intent(out) :: mu_c(mgncol,nlev)
188 : real(r8), intent(out) :: lambda_c(mgncol,nlev)
189 : real(r8), intent(out) :: lambda_r(mgncol,nlev)
190 : real(r8), intent(out) :: n0r(mgncol,nlev)
191 :
192 : real(r8), intent(out) :: qcin_new(mgncol,nlev)
193 : real(r8), intent(out) :: ncin_new(mgncol,nlev)
194 : real(r8), intent(out) :: qrin_new(mgncol,nlev)
195 : real(r8), intent(out) :: nrin_new(mgncol,nlev)
196 : real(r8), intent(out) :: gmnnn_lmnnn_TAU(mgncol,nlev)
197 :
198 : ! Local variables
199 :
200 : integer :: i,k,n,lcl
201 0 : integer :: cutoff_amk(mgncol,nlev),cutoff(mgncol,nlev)
202 :
203 : real(r8) :: all_gmnnn,all_lmnnn
204 : real(r8) :: qscl
205 :
206 0 : real(r8) :: qcin_old(mgncol,nlev)
207 0 : real(r8) :: ncin_old(mgncol,nlev)
208 0 : real(r8) :: qrin_old(mgncol,nlev)
209 0 : real(r8) :: nrin_old(mgncol,nlev)
210 :
211 0 : real(r8) :: amk0(mgncol,nlev,ncd)
212 0 : real(r8) :: ank0(mgncol,nlev,ncd)
213 0 : real(r8) :: gnnnn(mgncol,nlev,ncd)
214 0 : real(r8) :: gmnnn(mgncol,nlev,ncd)
215 0 : real(r8) :: lnnnn(mgncol,nlev,ncd)
216 0 : real(r8) :: lmnnn(mgncol,nlev,ncd)
217 0 : real(r8) :: gnnnn0(mgncol,nlev,ncd)
218 0 : real(r8) :: gmnnn0(mgncol,nlev,ncd)
219 0 : real(r8) :: lnnnn0(mgncol,nlev,ncd)
220 0 : real(r8) :: lmnnn0(mgncol,nlev,ncd)
221 :
222 : integer, parameter :: sub_step = 60
223 :
224 : !$acc data create (cutoff_amk,cutoff,qcin_old,ncin_old,qrin_old, &
225 : !$acc nrin_old,amk0,ank0,gnnnn,gmnnn,lnnnn,lmnnn, &
226 : !$acc gnnnn0,gmnnn0,lnnnn0,lmnnn0)
227 :
228 : !$acc parallel vector_length(VLENS) default(present)
229 : !$acc loop gang vector collapse(2)
230 0 : do k=1,nlev
231 0 : do i=1,mgncol
232 0 : cutoff(i,k) = cutoff_id - 1
233 : end do
234 : end do
235 : !$acc end parallel
236 :
237 : ! First make bins from cam size distribution (bins are diagnostic)
238 :
239 : call cam_bin_distribute(qc,qr,qcin,ncin,qrin,nrin,mu_c,lambda_c, &
240 : lambda_r,n0r,lcldm,precip_frac,scale_qc, &
241 : scale_nc,scale_qr,scale_nr,amk_c,ank_c, &
242 0 : amk_r,ank_r,amk,ank,cutoff_amk,mgncol,nlev)
243 :
244 : !$acc parallel vector_length(VLENS) default(present)
245 : !$acc loop gang vector collapse(2)
246 0 : do k=1,nlev
247 0 : do i=1,mgncol
248 0 : if ( cutoff_amk(i,k) > 0 ) then
249 0 : cutoff(i,k) = cutoff_amk(i,k)
250 : end if
251 : end do
252 : end do
253 : !$acc end parallel
254 :
255 : !Then call the subroutines that actually do the calculations. The inputs/ouputs are described in comments below.
256 :
257 : !This part of the code needs to be called each time for each process rate calculation
258 : ! (i.e., for each sampled cloud/rain gamma distribution):
259 :
260 : ! note: variables passed to compute_column_params are all inputs,
261 : ! outputs from this subroutine are stored as global variables
262 :
263 : ! inputs: t --> input air temperature (K)
264 : ! rho --> input air density (kg/m^3)
265 : ! medge --> bin mass boundary (g)
266 : ! amk --> array of bin mass mixing ratio, i.e., the input drop mass distribution (kg/kg)
267 : ! ank --> array of bin number mixing ratio, i.e., the input drop number distribution (kg^-1)
268 :
269 : ! inputs: medge --> bin mass boundary (g), same as above
270 :
271 : ! outputs: gnnnn --> bin number mass mixing tendency gain, array in bins (#/cm^3/s)
272 : ! lnnnn --> bin number mass mixing tendency loss, array in bins (#/cm^3/s)
273 : ! gmnnn --> bin mass mixing ratio tendency gain, array in bins (g/cm^3/s)
274 : ! lmnnn --> bin mass mixing ratio tendency loss, array in bins (g/cm^3/s)
275 :
276 :
277 : ! Call Kernel
278 :
279 : !$acc parallel vector_length(VLENS) default(present)
280 : !$acc loop gang vector collapse(2)
281 0 : do k=1,nlev
282 0 : do i=1,mgncol
283 0 : qcin_new(i,k) = 0._r8
284 0 : ncin_new(i,k) = 0._r8
285 0 : qrin_new(i,k) = 0._r8
286 0 : nrin_new(i,k) = 0._r8
287 :
288 0 : qcin_old(i,k) = 0._r8
289 0 : ncin_old(i,k) = 0._r8
290 0 : qrin_old(i,k) = 0._r8
291 0 : nrin_old(i,k) = 0._r8
292 :
293 0 : qctend_TAU(i,k) = 0._r8
294 0 : nctend_TAU(i,k) = 0._r8
295 0 : qrtend_TAU(i,k) = 0._r8
296 0 : nrtend_TAU(i,k) = 0._r8
297 : end do
298 : end do
299 : !$acc end parallel
300 :
301 : !$acc parallel vector_length(VLENS) default(present)
302 : !$acc loop gang vector collapse(3)
303 0 : do lcl=1,ncd
304 0 : do k=1,nlev
305 0 : do i=1,mgncol
306 0 : gnnnn(i,k,lcl) = 0._r8
307 0 : gmnnn(i,k,lcl) = 0._r8
308 0 : lnnnn(i,k,lcl) = 0._r8
309 0 : lmnnn(i,k,lcl) = 0._r8
310 : end do
311 : end do
312 : end do
313 : !$acc end parallel
314 :
315 : ! update qc, nc, qr, nr
316 :
317 : !$acc parallel vector_length(VLENS) default(present)
318 : !$acc loop gang vector collapse(2)
319 0 : do k=1,nlev
320 0 : do i=1,mgncol
321 : !$acc loop seq
322 0 : do lcl=1,ncd
323 0 : amk0(i,k,lcl) = amk(i,k,lcl)
324 0 : ank0(i,k,lcl) = ank(i,k,lcl)
325 : end do
326 : ! substep bin code
327 : !$acc loop seq
328 0 : do n=1,sub_step
329 0 : call compute_coll_params(rho(i,k),medge,amk0(i,k,1:ncd),ank0(i,k,1:ncd),gnnnn0(i,k,1:ncd),gmnnn0(i,k,1:ncd),lnnnn0(i,k,1:ncd),lmnnn0(i,k,1:ncd))
330 :
331 0 : all_gmnnn=0._r8
332 0 : all_lmnnn=0._r8
333 : !scaling gmnnn, lmnnn
334 : !$acc loop seq
335 0 : do lcl=1,ncd
336 0 : all_gmnnn = all_gmnnn+gmnnn0(i,k,lcl)
337 0 : all_lmnnn = all_lmnnn+lmnnn0(i,k,lcl)
338 : end do
339 :
340 0 : if ( (all_gmnnn == 0._r8) .or. (all_lmnnn == 0._r8) ) then
341 : !$acc loop seq
342 0 : do lcl=1,ncd
343 0 : gmnnn0(i,k,lcl) = 0._r8
344 0 : lmnnn0(i,k,lcl) = 0._r8
345 : end do
346 : else
347 : !$acc loop seq
348 0 : do lcl=1,ncd
349 0 : lmnnn0(i,k,lcl) = lmnnn0(i,k,lcl)*(all_gmnnn/all_lmnnn)
350 : end do
351 : end if
352 :
353 : !$acc loop seq
354 0 : do lcl=1,ncd
355 0 : amk0(i,k,lcl) = amk0(i,k,lcl)+(gmnnn0(i,k,lcl)-lmnnn0(i,k,lcl))*1.e3_r8/ &
356 0 : rho(i,k)*deltatin/dble(sub_step)
357 : ank0(i,k,lcl) = ank0(i,k,lcl)+(gnnnn0(i,k,lcl)-lnnnn0(i,k,lcl))*1.e6_r8/ &
358 0 : rho(i,k)*deltatin/dble(sub_step)
359 0 : gmnnn(i,k,lcl) = gmnnn(i,k,lcl)+gmnnn0(i,k,lcl)/sub_step
360 0 : gnnnn(i,k,lcl) = gnnnn(i,k,lcl)+gnnnn0(i,k,lcl)/sub_step
361 0 : lmnnn(i,k,lcl) = lmnnn(i,k,lcl)+lmnnn0(i,k,lcl)/sub_step
362 0 : lnnnn(i,k,lcl) = lnnnn(i,k,lcl)+lnnnn0(i,k,lcl)/sub_step
363 : end do
364 : end do ! end of loop "sub_step"
365 :
366 : ! cloud water
367 : !$acc loop seq
368 0 : do lcl=1,cutoff(i,k)
369 0 : qcin_old(i,k) = qcin_old(i,k)+amk(i,k,lcl)
370 0 : ncin_old(i,k) = ncin_old(i,k)+ank(i,k,lcl)
371 0 : qcin_new(i,k) = qcin_new(i,k)+(gmnnn(i,k,lcl)-lmnnn(i,k,lcl))*1.e3_r8/rho(i,k)*deltatin
372 0 : ncin_new(i,k) = ncin_new(i,k)+(gnnnn(i,k,lcl)-lnnnn(i,k,lcl))*1.e6_r8/rho(i,k)*deltatin
373 0 : qctend_TAU(i,k) = qctend_TAU(i,k)+(amk0(i,k,lcl)-amk(i,k,lcl))/deltatin
374 0 : nctend_TAU(i,k) = nctend_TAU(i,k)+(ank0(i,k,lcl)-ank(i,k,lcl))/deltatin
375 0 : gmnnn_lmnnn_TAU(i,k) = gmnnn_lmnnn_TAU(i,k)+gmnnn(i,k,lcl)-lmnnn(i,k,lcl)
376 : end do
377 :
378 : ! rain
379 : !$acc loop seq
380 0 : do lcl=cutoff(i,k)+1,ncd
381 0 : qrin_old(i,k) = qrin_old(i,k)+amk(i,k,lcl)
382 0 : nrin_old(i,k) = nrin_old(i,k)+ank(i,k,lcl)
383 0 : qrin_new(i,k) = qrin_new(i,k)+(gmnnn(i,k,lcl)-lmnnn(i,k,lcl))*1.e3_r8/rho(i,k)*deltatin
384 0 : nrin_new(i,k) = nrin_new(i,k)+(gnnnn(i,k,lcl)-lnnnn(i,k,lcl))*1.e6_r8/rho(i,k)*deltatin
385 0 : qrtend_TAU(i,k) = qrtend_TAU(i,k)+(amk0(i,k,lcl)-amk(i,k,lcl))/deltatin
386 0 : nrtend_TAU(i,k) = nrtend_TAU(i,k)+(ank0(i,k,lcl)-ank(i,k,lcl))/deltatin
387 0 : gmnnn_lmnnn_TAU(i,k) = gmnnn_lmnnn_TAU(i,k)+gmnnn(i,k,lcl)-lmnnn(i,k,lcl)
388 : end do
389 :
390 : !$acc loop seq
391 0 : do lcl=1,ncd
392 0 : amk_out(i,k,lcl) = amk(i,k,lcl) + (gmnnn(i,k,lcl)-lmnnn(i,k,lcl))*1.e3_r8/rho(i,k)*deltatin
393 0 : ank_out(i,k,lcl) = ank(i,k,lcl) + (gnnnn(i,k,lcl)-lnnnn(i,k,lcl))*1.e6_r8/rho(i,k)*deltatin
394 : end do
395 :
396 0 : qcin_new(i,k) = qcin_new(i,k)+qcin_old(i,k)
397 0 : ncin_new(i,k) = ncin_new(i,k)+ncin_old(i,k)
398 0 : qrin_new(i,k) = qrin_new(i,k)+qrin_old(i,k)
399 0 : nrin_new(i,k) = nrin_new(i,k)+nrin_old(i,k)
400 : end do
401 : end do
402 : !$acc end parallel
403 :
404 : ! Conservation checks
405 : ! AG: Added May 2023
406 :
407 : !$acc parallel vector_length(VLENS) default(present)
408 : !$acc loop gang vector collapse(2)
409 0 : do k=1,nlev
410 0 : do i=1,mgncol
411 :
412 : ! First make sure all not negative
413 0 : qcin_new(i,k)=max(qcin_new(i,k),0._r8)
414 0 : ncin_new(i,k)=max(ncin_new(i,k),0._r8)
415 0 : qrin_new(i,k)=max(qrin_new(i,k),0._r8)
416 0 : nrin_new(i,k)=max(nrin_new(i,k),0._r8)
417 :
418 : ! Now adjust so that sign is correct. qc_new,nc_new <= input, qr_new >= input
419 : ! NOte that due to self collection nr can be larger or smaller than input....
420 : ! Makes above check redundant I think.
421 :
422 0 : qcin_new(i,k)=min(qcin_new(i,k),qcin(i,k))
423 0 : ncin_new(i,k)=min(ncin_new(i,k),ncin(i,k))
424 0 : qrin_new(i,k)=max(qrin_new(i,k),qrin(i,k))
425 :
426 : ! Next scale mass...so output qc+qr is the same as input
427 :
428 0 : if ( (qcin_new(i,k)+qrin_new(i,k)) > 0._r8 ) then
429 0 : qscl = (qcin(i,k)+qrin(i,k))/(qcin_new(i,k)+qrin_new(i,k))
430 : else
431 : qscl = 0._r8
432 : end if
433 0 : qcin_new(i,k) = qcin_new(i,k) * qscl
434 0 : qrin_new(i,k) = qrin_new(i,k) * qscl
435 :
436 : ! Now zero nr,nc if either small or no mass?
437 :
438 0 : if ( qcin_new(i,k) < qsmall ) then
439 0 : ncin_new(i,k) = 0._r8
440 : end if
441 :
442 0 : if ( qrin_new(i,k) < qsmall ) then
443 0 : nrin_new(i,k) = 0._r8
444 : end if
445 :
446 : !Finally add number if mass but no (or super small) number
447 :
448 0 : if ( qcin_new(i,k) > qsmall .and. ncin_new(i,k) < qsmall ) then
449 0 : ncin_new(i,k) = qcin_new(i,k)/m1
450 : end if
451 :
452 0 : if ( qrin_new(i,k) > qsmall .and. nrin_new(i,k) < qsmall) then
453 0 : nrin_new(i,k) = qrin_new(i,k)/m1
454 : end if
455 :
456 : ! Then recalculate tendencies based on difference
457 : ! Clip tendencies for cloud (qc,nc) to be <= 0.
458 : ! Qrtend is not used in pumas (-qctend is used) but clip that too).
459 : ! Nr tend can be muliply signed.
460 :
461 0 : qctend_TAU(i,k)= min((qcin_new(i,k) - qcin(i,k)) / deltatin,0._r8)
462 0 : nctend_TAU(i,k)= min((ncin_new(i,k) - ncin(i,k)) / deltatin,0._r8)
463 0 : qrtend_TAU(i,k)= max((qrin_new(i,k) - qrin(i,k)) / deltatin,0._r8)
464 0 : nrtend_TAU(i,k)= (nrin_new(i,k) - nrin(i,k)) / deltatin
465 :
466 : end do
467 : end do
468 : !$acc end parallel
469 :
470 : !$acc end data
471 :
472 0 : end subroutine pumas_stochastic_collect_tau_tend
473 :
474 0 : subroutine cam_bin_distribute(qc_all,qr_all,qc,nc,qr,nr,mu_c,lambda_c, &
475 0 : lambda_r,n0r,lcldm,precip_frac,scale_qc, &
476 0 : scale_nc,scale_qr,scale_nr,amk_c,ank_c, &
477 0 : amk_r,ank_r,amk,ank,cutoff_amk,mgncol,nlev)
478 :
479 : implicit none
480 :
481 : integer, intent(in) :: mgncol,nlev
482 : real(r8), dimension(mgncol,nlev), intent(in) :: qc_all,qr_all,qc,nc,qr,nr,mu_c, &
483 : lambda_c,lambda_r,n0r,lcldm, &
484 : precip_frac
485 : real(r8), dimension(mgncol,nlev,ncd), intent(out) :: amk_c,ank_c,amk_r,ank_r,amk,ank
486 : real(r8), dimension(mgncol,nlev), intent(out) :: scale_nc,scale_qc,scale_nr,scale_qr
487 : integer, dimension(mgncol,nlev), intent(out) :: cutoff_amk
488 :
489 : ! Local variables
490 :
491 : integer :: i,j,k
492 : real(r8) :: phi
493 : integer :: id_max_qc, id_max_qr
494 : real(r8) :: max_qc, max_qr, min_amk
495 :
496 : !$acc parallel vector_length(VLENS) default(present)
497 : !$acc loop gang vector collapse(3)
498 0 : do j=1,ncd
499 0 : do k=1,nlev
500 0 : do i=1,mgncol
501 0 : ank_c(i,k,j) = 0._r8
502 0 : amk_c(i,k,j) = 0._r8
503 0 : ank_r(i,k,j) = 0._r8
504 0 : amk_r(i,k,j) = 0._r8
505 0 : ank(i,k,j) = 0._r8
506 0 : amk(i,k,j) = 0._r8
507 : end do
508 : end do
509 : end do
510 : !$acc end parallel
511 :
512 : !$acc parallel vector_length(VLENS) default(present)
513 : !$acc loop gang vector collapse(2)
514 0 : do k=1,nlev
515 0 : do i=1,mgncol
516 0 : scale_nc(i,k) = 0._r8
517 0 : scale_qc(i,k) = 0._r8
518 0 : scale_nr(i,k) = 0._r8
519 0 : scale_qr(i,k) = 0._r8
520 0 : cutoff_amk(i,k) = 0
521 :
522 0 : id_max_qc = 0
523 0 : id_max_qr = 0
524 0 : max_qc = 0._r8
525 0 : max_qr = 0._r8
526 :
527 : ! cloud water, nc in #/m3 --> #/cm3
528 0 : if ( (qc_all(i,k) > qsmall) .and. (qc(i,k) > qsmall) ) then
529 : !$acc loop seq
530 0 : do j=1,ncd
531 : phi = nc(i,k)*lambda_c(i,k)**(mu_c(i,k)+1._r8)/ &
532 0 : gamma(mu_c(i,k)+1._r8)*(diammean(j)*1.e-2_r8)**mu_c(i,k)* &
533 0 : exp(-lambda_c(i,k)*diammean(j)*1.e-2_r8) ! D cm --> m
534 0 : ank_c(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8 ! D cm --> m
535 0 : amk_c(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8*mmean(j)*1.e-3_r8 ! mass in bin g --> kg
536 0 : scale_nc(i,k) = scale_nc(i,k)+ank_c(i,k,j)
537 0 : scale_qc(i,k) = scale_qc(i,k)+amk_c(i,k,j)
538 : end do
539 0 : scale_nc(i,k) = scale_nc(i,k)/nc(i,k)
540 0 : scale_qc(i,k) = scale_qc(i,k)/qc(i,k)
541 :
542 : !$acc loop seq
543 0 : do j=1,ncd
544 0 : ank_c(i,k,j) = ank_c(i,k,j)/scale_nc(i,k)*lcldm(i,k)
545 0 : amk_c(i,k,j) = amk_c(i,k,j)/scale_qc(i,k)*lcldm(i,k)
546 0 : if ( amk_c(i,k,j) > max_qc ) then
547 0 : id_max_qc = j
548 0 : max_qc = amk_c(i,k,j)
549 : end if
550 : end do
551 : end if
552 :
553 : ! rain drop
554 0 : if ( (qr_all(i,k) > qsmall) .and. (qr(i,k) > qsmall) ) then
555 : !$acc loop seq
556 0 : do j=1,ncd
557 0 : phi = n0r(i,k)*exp(-lambda_r(i,k)*diammean(j)*1.e-2_r8) ! D cm --> m
558 0 : ank_r(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8 ! D cm --> m
559 0 : amk_r(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8*mmean(j)*1.e-3_r8
560 0 : scale_nr(i,k) = scale_nr(i,k) + ank_r(i,k,j)
561 0 : scale_qr(i,k) = scale_qr(i,k) + amk_r(i,k,j)
562 : end do
563 0 : scale_nr(i,k) = scale_nr(i,k)/nr(i,k)
564 0 : scale_qr(i,k) = scale_qr(i,k)/qr(i,k)
565 :
566 : !$acc loop seq
567 0 : do j=1,ncd
568 0 : ank_r(i,k,j) = ank_r(i,k,j)/scale_nr(i,k)*precip_frac(i,k)
569 0 : amk_r(i,k,j) = amk_r(i,k,j)/scale_qr(i,k)*precip_frac(i,k)
570 0 : if ( amk_r(i,k,j) > max_qr ) then
571 0 : id_max_qr = j
572 0 : max_qr = amk_r(i,k,j)
573 : end if
574 : end do
575 : end if
576 :
577 : !$acc loop seq
578 0 : do j=1,ncd
579 0 : amk(i,k,j) = amk_c(i,k,j) + amk_r(i,k,j)
580 0 : ank(i,k,j) = ank_c(i,k,j) + ank_r(i,k,j)
581 : end do
582 :
583 0 : if ( (id_max_qc > 0) .and. (id_max_qr > 0) ) then
584 0 : if ( (max_qc/max_qr < 10._r8) .or. (max_qc/max_qr > 0.1_r8) ) then
585 0 : min_amk = amk(i,k,id_max_qc)
586 : !$acc loop seq
587 0 : do j=id_max_qc,id_max_qr
588 0 : if ( amk(i,k,j) <= min_amk ) then
589 0 : cutoff_amk(i,k) = j
590 0 : min_amk = amk(i,k,j)
591 : end if
592 : end do
593 : end if
594 : end if
595 : end do ! end of loop "mgncol"
596 : end do ! end of loop "nlev"
597 : !$acc end parallel
598 :
599 : !input: qc,nc,qr,nr, medge (bin edges). May also need # bins?
600 : !output: amk, ank (mixing ratio and number in each bin)
601 :
602 : !this part will take a bit of thinking about.
603 : !use size distribution parameters (mu, lambda) to generate the values at discrete size points
604 : !need to also ensure mass conservation
605 :
606 0 : end subroutine cam_bin_distribute
607 :
608 : ! here are the subroutines called above that actually do the collision-coalescence calculations:
609 :
610 : ! The Kernel is from Jerry from many moons ago (included)
611 :
612 : ! I read in the file data and multiply by the summed mass of the individual bins
613 : ! (with a factor of 1.5 so that the values represent the middle of the bin
614 :
615 : ! 941 FORMAT(2X,E12.5)
616 : ! READ(iunit,941) KNN(IDD,JDD)
617 : ! KNN(IDD,JDD)=(XK_GR(IDD)*kkfac+XK_GR(JDD)*kkfac)*KNN(IDD,JDD)
618 :
619 : !where idd and jdd are the indexes for the bins and xk_gr is the mass of drops in a bin in grams
620 : !
621 :
622 : !************************************************************************************
623 : ! Setup variables needed for collection
624 : ! Either pass in or define globally the following variables
625 : ! tbase(height) - temperature in K as a function of height
626 : ! rhon(height) - air density as a function of height in kg/m^3
627 : ! xk_gr(bins) - mass of single drop in each bin in grams
628 : ! lsmall - small number
629 : ! QC - mass mixing ratio in kg/kg
630 : ! QN - number mixing ratio in #/kg
631 : ! All parameters are defined to be global in my version so that they are readily available throughout the code:
632 : ! SMN0,SNN0,SMCN,APN,AMN2,AMN3,PSIN,FN,FPSIN,XPSIN,HPSIN,FN2,XXPSIN (all arrays of drop bins)
633 : !************************************************************************************
634 :
635 : !AG: Global arrays need to be passed around I think? Right now at the module level. Is that okay?
636 :
637 0 : SUBROUTINE COMPUTE_COLL_PARAMS(rhon,xk_gr,qc,qn,gnnnn,gmnnn,lnnnn,lmnnn)
638 :
639 : !$acc routine seq
640 :
641 : IMPLICIT NONE
642 :
643 : ! variable declarations (added by hm, 020118)
644 : ! note: vertical array looping is stripped out, this subroutine operates
645 : ! only on LOCAL values
646 :
647 : real(r8), dimension(ncd) :: qc,qn
648 : real(r8), dimension(ncdp) :: xk_gr
649 : real(r8) :: tbase,rhon
650 : integer :: lk
651 : integer :: l
652 : real(r8), parameter :: lsmall = 1.e-12_r8
653 : real(r8), dimension(ncd) :: smn0,snn0,smcn,amn2,amn3,psin,fn,fpsin, &
654 : xpsin,hpsin,fn2,xxpsin
655 : real(r8) :: apn
656 :
657 : real(r8), dimension(ncd) :: gnnnn,gmnnn,lnnnn,lmnnn
658 : integer :: lm1,ll
659 :
660 0 : lk=ncd
661 :
662 0 : DO L=1,LK
663 0 : SMN0(L)=QC(L)*RHON/1.E3_r8
664 0 : SNN0(L)=QN(L)*RHON/1.E6_r8
665 :
666 0 : IF(SMN0(L).LT.lsmall.OR.SNN0(L).LT.lsmall)THEN
667 0 : SMN0(L)=0.0_r8
668 0 : SNN0(L)=0.0_r8
669 : ENDIF
670 : ENDDO
671 :
672 0 : DO L=1,LK
673 0 : IF(SMN0(L) .gt. 0._r8.AND.SNN0(L) .gt. 0._r8)THEN
674 0 : SMCN(L)=SMN0(L)/SNN0(L)
675 0 : IF((SMCN(L) .GT. 2._r8*XK_GR(L)))THEN
676 0 : SMCN(L) = (2._r8*XK_GR(L))
677 : ENDIF
678 0 : IF((SMCN(L) .LT. XK_GR(L)))THEN
679 0 : SMCN(L) = XK_GR(L)
680 : ENDIF
681 : ELSE
682 0 : SMCN(L)=0._r8
683 : ENDIF
684 0 : IF (SMCN(L).LT.XK_GR(L).OR.SMCN(L).GT.(2._r8*XK_GR(L)).OR.SMCN(L).EQ.0.0_r8)THEN
685 : APN=1.0_r8
686 : ELSE
687 0 : APN=0.5_r8*(1._r8+3._r8*(XK_GR(L)/SMCN(L))-2._r8*((XK_GR(L)/SMCN(L))**2._r8))
688 : ENDIF
689 :
690 0 : IF(SNN0(L) .GT. LSMALL)THEN
691 0 : AMN2(L)=APN*SMN0(L)*SMN0(L)/SNN0(L)
692 0 : AMN3(L)=APN*APN*APN*SMN0(L)*SMN0(L)*SMN0(L)/(SNN0(L)*SNN0(L))
693 : ELSE
694 0 : AMN2(L)=0._r8
695 0 : AMN3(L)=0._r8
696 : ENDIF
697 :
698 0 : IF(SMCN(L).LT.XK_GR(L))THEN
699 0 : PSIN(L)=0.0_r8
700 0 : FN(L)=2._r8*SNN0(L)/XK_GR(L)
701 : ELSE
702 0 : IF(SMCN(L).GT.(2._r8*XK_GR(L)))THEN
703 0 : FN(L)=0.0_r8
704 0 : PSIN(L)=2._r8*SNN0(L)/XK_GR(L)
705 : ELSE
706 0 : PSIN(L)=2._r8/XK_GR(L)*(SMN0(L)/XK_GR(L)-SNN0(L))
707 0 : FN(L)=2._r8/XK_GR(L)*(2._r8*SNN0(L)-SMN0(L)/XK_GR(L))
708 : ENDIF
709 : ENDIF
710 :
711 0 : IF(SNN0(L).LT.LSMALL.OR.SMN0(L).LT.LSMALL)THEN
712 0 : PSIN(L)=0.0_r8
713 0 : FN(L)=0.0_r8
714 : ENDIF
715 :
716 0 : FPSIN(L)=0.5_r8/XK_GR(L)*(PSIN(L)-FN(L))
717 0 : XPSIN(L)=2._r8*XK_GR(L)*PSIN(L)
718 0 : HPSIN(L)=PSIN(L)-0.5_r8*FN(L)
719 0 : FN2(L)=FN(L)/2._r8
720 :
721 0 : IF(L.GT.1)THEN
722 0 : XXPSIN(L)=XK_GR(L)*PSIN(L-1)
723 : ENDIF
724 : ENDDO
725 :
726 : !************************************************************************************
727 : ! Compute collision coalescence
728 : ! Either pass in or define globally the following variables
729 : ! Gain terms begin with G, loss terms begin with L
730 : ! Second letter defines mass (M) or number (N)
731 : ! Third and fourth letters define the types of particles colling, i.e., NN means drops colliding with drops
732 : ! Last letter defines the category the new particles go into, in this case just N for liquid drops
733 : ! The resulting rates are in units of #/cm^3/s and g/cm^3/s
734 : ! Relies on predefined kernel array KNN(bins,bins) - see top of this file
735 : !************************************************************************************
736 :
737 0 : GMNNN = 0._r8
738 0 : GNNNN = 0._r8
739 0 : LMNNN = 0._r8
740 0 : LNNNN = 0._r8
741 : ! remove verical array index, calculate gain/loss terms locally
742 :
743 0 : DO L=3,LK-1
744 0 : LM1=L-1
745 0 : DO LL=1,L-2
746 0 : GNNNN(L)=GNNNN(L)+(PSIN(LM1)*SMN0(LL)-FPSIN(LM1)*AMN2(LL))*KNN(LM1,LL)
747 0 : GMNNN(L)=GMNNN(L)+(XK_GR(L)*PSIN(LM1)*SMN0(LL)+FN2(LM1)*AMN2(LL)-FPSIN(LM1)*AMN3(LL))*KNN(LM1,LL)
748 : ENDDO
749 : ENDDO
750 :
751 0 : DO L=2,LK-1
752 0 : LM1=L-1
753 0 : GNNNN(L)=GNNNN(L)+0.5_r8*SNN0(LM1)*SNN0(LM1)*KNN(LM1,LM1)
754 0 : GMNNN(L)=GMNNN(L)+0.5_r8*(SNN0(LM1)*SMN0(LM1)+SMN0(LM1)*SNN0(LM1))*KNN(LM1,LM1)
755 0 : DO LL=1,L-1
756 0 : LNNNN(L)=LNNNN(L)+(PSIN(L)*SMN0(LL)-FPSIN(L)*AMN2(LL))*KNN(L,LL)
757 0 : GMNNN(L)=GMNNN(L)+(SMN0(LL)*SNN0(L)-PSIN(L)*AMN2(LL)+FPSIN(L)*AMN3(LL))*KNN(L,LL)
758 0 : LMNNN(L)=LMNNN(L)+(XPSIN(L)*SMN0(LL)-HPSIN(L)*AMN2(LL))*KNN(L,LL)
759 : ENDDO
760 : ENDDO
761 :
762 0 : DO L=1,LK-1
763 0 : DO LL=L,LK-1
764 0 : LNNNN(L)=LNNNN(L)+(SNN0(LL)*SNN0(L))*KNN(LL,L)
765 0 : LMNNN(L)=LMNNN(L)+(SNN0(LL)*SMN0(L))*KNN(LL,L)
766 : ENDDO
767 : ENDDO
768 :
769 0 : END SUBROUTINE COMPUTE_COLL_PARAMS
770 :
771 :
772 : end module pumas_stochastic_collect_tau
773 :
774 :
|
