| blob | 4413b02961bc4734914dc9ea248f17cd23466a3b |
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for references and terms of use
8 !**********************************************************************************
10 module module_mosaic_cloudchem
14 integer, parameter :: l_so4_aqyy = 1
15 integer, parameter :: l_no3_aqyy = 2
16 integer, parameter :: l_cl_aqyy = 3
17 integer, parameter :: l_nh4_aqyy = 4
18 integer, parameter :: l_na_aqyy = 5
19 integer, parameter :: l_oin_aqyy = 6
20 integer, parameter :: l_bc_aqyy = 7
21 integer, parameter :: l_oc_aqyy = 8
23 integer, parameter :: nyyy = 8
27 contains
31 !-----------------------------------------------------------------------
32 subroutine mosaic_cloudchem_driver( &
33 id, ktau, ktauc, dtstepc, config_flags, &
34 p_phy, t_phy, rho_phy, alt, &
35 cldfra, ph_no2, &
36 moist, chem, ph_cw, &
37 gas_aqfrac, numgas_aqfrac, &
38 ids,ide, jds,jde, kds,kde, &
39 ims,ime, jms,jme, kms,kme, &
40 its,ite, jts,jte, kts,kte )
42 use module_state_description, only: &
43 num_moist, num_chem, p_qc
45 use module_configure, only: grid_config_rec_type
47 use module_data_mosaic_asect, only: cw_phase, nphase_aer
49 use module_data_mosaic_other, only: k_pegbegin, name
51 use module_mosaic_driver, only: mapaer_tofrom_host
54 implicit none
56 ! subr arguments
57 integer, intent(in) :: &
58 id, ktau, ktauc, &
59 numgas_aqfrac, &
60 ids, ide, jds, jde, kds, kde, &
61 ims, ime, jms, jme, kms, kme, &
62 its, ite, jts, jte, kts, kte
63 ! id - domain index
64 ! ktau - time step number
65 ! ktauc - gas and aerosol chemistry time step number
66 ! numgas_aqfrac - last dimension of gas_aqfrac
68 ! [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for 'domain'
69 ! [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for 'memory'
70 ! Most arrays that are arguments to chem_driver
71 ! are dimensioned with these spatial indices.
72 ! [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for 'tile'
73 ! chem_driver and routines under it do calculations
74 ! over these spatial indices.
76 type(grid_config_rec_type), intent(in) :: config_flags
77 ! config_flags - configuration and control parameters
79 real, intent(in) :: &
80 dtstepc
81 ! dtstepc - time step for gas and aerosol chemistry(s)
83 real, intent(in), &
84 dimension( ims:ime, kms:kme, jms:jme ) :: &
85 p_phy, t_phy, rho_phy, alt, cldfra, ph_no2
86 ! p_phy - air pressure (Pa)
87 ! t_phy - temperature (K)
88 ! rho_phy - moist air density (kg/m^3)
89 ! alt - dry air specific volume (m^3/kg)
90 ! cldfra - cloud fractional area (0-1)
91 ! ph_no2 - no2 photolysis rate (1/min)
93 real, intent(in), &
94 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
95 moist
96 ! moist - mixing ratios of moisture species (water vapor,
97 ! cloud water, ...) (kg/kg for mass species, #/kg for number species)
99 real, intent(inout), &
100 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
101 chem
102 ! chem - mixing ratios of trace gas and aerosol species (ppm for gases,
103 ! ug/kg for aerosol mass species, #/kg for aerosol number species)
105 real, intent(inout), &
106 dimension( ims:ime, kms:kme, jms:jme, numgas_aqfrac ) :: &
107 gas_aqfrac
108 ! gas_aqfrac - fraction (0-1) of gas that is dissolved in cloud water
110 real, intent(out), &
111 OPTIONAL, &
112 dimension( ims:ime, kms:kme, jms:jme ) :: &
113 ph_cw
115 ! local variables
116 integer :: it, jt, kt, kpeg, k_pegshift, l, mpeg
117 integer :: icase
118 integer :: igaschem_onoff, iphotol_onoff, iradical_onoff
120 real :: gas_aqfrac_box(numgas_aqfrac)
121 real :: ph_aq_box
122 real, parameter :: qcldwtr_cutoff = 1.0e-6
123 real :: qcldwtr
126 ! check that cw_phase is active
127 if ((cw_phase .le. 0) .or. (cw_phase .gt. nphase_aer)) then
128 print *, '*** mosaic_cloudchem_driver - cw_phase not active'
129 return
130 end if
132 print 93010, 'entering mosaic_cloudchem_driver - ktau =', ktau
134 icase = 0
136 ! iphotol_onoff = 1 if photolysis rate calcs are on; 0 if off
137 iphotol_onoff = 0
138 if (config_flags%phot_opt .gt. 0) iphotol_onoff = 1
139 ! igaschem_onoff = 1 if gas-phase chemistry is on; 0 if off
140 igaschem_onoff = 0
141 if (config_flags%gaschem_onoff .gt. 0) igaschem_onoff = 1
143 ! iradical_onoff turns aqueous radical chemistry on/off
144 ! set iradical_onoff=0 if either photolysis or gas-phase chem are off
145 if ((igaschem_onoff .le. 0) .or. (iphotol_onoff .le. 0)) then
146 iradical_onoff = 0
147 else
148 iradical_onoff = 1
149 end if
150 ! following line turns aqueous radical chem off unconditionally
151 iradical_onoff = 0
152 if ( config_flags%mozart_ph_diag .eq. 1 ) then
153 ! Initialize pH of CW ph_cw to a FillValue value
154 do jt = jts, jte
155 do kt = kts, kte
156 do it = its, ite
157 ph_cw(it,kt,jt) = -9999.
158 end do
159 end do
160 end do
161 end if
162 do 3920 jt = jts, jte
163 do 3910 it = its, ite
165 do 3800 kt = kts, kte
167 qcldwtr = moist(it,kt,jt,p_qc)
168 if (qcldwtr .le. qcldwtr_cutoff) goto 3800
171 k_pegshift = k_pegbegin - kts
172 kpeg = kt + k_pegshift
173 mpeg = 1
174 icase = icase + 1
176 ! detailed dump for debugging
177 if (ktau .eq. -13579) then
178 ! if ((ktau .eq. 30) .and. (it .eq. 23) .and. &
179 ! (jt .eq. 1) .and. (kt .eq. 11)) then
180 call mosaic_cloudchem_dumpaa( &
181 id, ktau, ktauc, dtstepc, config_flags, &
182 p_phy, t_phy, rho_phy, alt, &
183 cldfra, ph_no2, &
184 moist, chem, &
185 gas_aqfrac, numgas_aqfrac, &
186 ids,ide, jds,jde, kds,kde, &
187 ims,ime, jms,jme, kms,kme, &
188 its,ite, jts,jte, kts,kte, &
189 qcldwtr_cutoff, &
190 it, jt, kt )
191 end if
193 ! map from wrf-chem 3d arrays to pegasus clm & sub arrays
194 call mapaer_tofrom_host( 0, &
195 ims,ime, jms,jme, kms,kme, &
196 its,ite, jts,jte, kts,kte, &
197 it, jt, kt, kt, &
198 num_moist, num_chem, moist, chem, &
199 t_phy, p_phy, rho_phy )
201 ! make call '1box' cloudchem routine
202 ! print 93010, 'calling mosaic_cloudchem_1 at ijk =', it, jt, kt
203 call mosaic_cloudchem_1box( &
204 id, ktau, ktauc, dtstepc, &
205 iphotol_onoff, iradical_onoff, &
206 ph_no2(it,kt,jt), &
207 ph_aq_box, gas_aqfrac_box, &
208 numgas_aqfrac, it, jt, kt, kpeg, mpeg, icase )
210 ! map back to wrf-chem 3d arrays
211 call mapaer_tofrom_host( 1, &
212 ims,ime, jms,jme, kms,kme, &
213 its,ite, jts,jte, kts,kte, &
214 it, jt, kt, kt, &
215 num_moist, num_chem, moist, chem, &
216 t_phy, p_phy, rho_phy )
218 gas_aqfrac(it,kt,jt,:) = gas_aqfrac_box(:)
219 if ( config_flags%mozart_ph_diag .eq. 1 ) then
220 ph_cw(it,kt,jt) = ph_aq_box
221 end if
222 3800 continue
224 3910 continue
225 3920 continue
227 print 93010, 'leaving mosaic_cloudchem_driver - ktau =', ktau, icase
228 93010 format( a, 8(1x,i6) )
230 return
231 end subroutine mosaic_cloudchem_driver
235 !-----------------------------------------------------------------------
236 subroutine mosaic_cloudchem_1box( &
237 id, ktau, ktauc, dtstepc, &
238 iphotol_onoff, iradical_onoff, &
239 photol_no2_box, &
240 ph_aq_box, gas_aqfrac_box, &
241 numgas_aqfrac, it, jt, kt, kpeg, mpeg, icase )
243 use module_state_description, only: &
244 num_moist, num_chem
246 use module_data_mosaic_asect, only: &
247 msectional, &
248 maxd_asize, maxd_atype, &
249 cw_phase, nsize_aer, ntype_aer, &
250 lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
251 lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
252 lptr_oin_aer, lptr_bc_aer, lptr_oc_aer
254 use module_data_mosaic_other, only: &
255 l2maxd, ltot2, rsub
257 use module_data_cmu_bulkaqchem, only: &
258 meqn1max
261 implicit none
263 ! subr arguments
264 integer, intent(in) :: &
265 id, ktau, ktauc, &
266 numgas_aqfrac, it, jt, kt, kpeg, mpeg, &
267 icase, iphotol_onoff, iradical_onoff
269 real, intent(in) :: &
270 dtstepc, photol_no2_box
272 real, intent(inout) :: ph_aq_box
274 real, intent(inout), dimension( numgas_aqfrac ) :: gas_aqfrac_box
276 ! local variables
277 integer :: iphase
278 integer :: icase_in, idecomp_hmsa_hso5, &
279 iradical_in, istat_aqop
281 integer :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
283 real :: co2_mixrat_in
284 real :: ph_cmuaq_cur
285 real :: photol_no2_in
286 real :: xprescribe_ph
288 real :: yaq_beg(meqn1max), yaq_end(meqn1max)
289 real :: rbox(l2maxd), rbox_sv1(l2maxd)
290 real :: rbulk_cwaer(nyyy,2)
292 real, dimension( maxd_asize, maxd_atype ) :: fr_partit_cw
295 !
296 ! set the lptr_yyy_cwaer
297 !
298 iphase = cw_phase
299 lptr_yyy_cwaer(:,:,l_so4_aqyy) = lptr_so4_aer(:,:,iphase)
300 lptr_yyy_cwaer(:,:,l_no3_aqyy) = lptr_no3_aer(:,:,iphase)
301 lptr_yyy_cwaer(:,:,l_cl_aqyy ) = lptr_cl_aer( :,:,iphase)
302 lptr_yyy_cwaer(:,:,l_nh4_aqyy) = lptr_nh4_aer(:,:,iphase)
303 lptr_yyy_cwaer(:,:,l_na_aqyy ) = lptr_na_aer( :,:,iphase)
304 lptr_yyy_cwaer(:,:,l_oin_aqyy) = lptr_oin_aer(:,:,iphase)
305 lptr_yyy_cwaer(:,:,l_bc_aqyy ) = lptr_bc_aer( :,:,iphase)
306 lptr_yyy_cwaer(:,:,l_oc_aqyy ) = lptr_oc_aer( :,:,iphase)
308 !
309 ! xfer from rsub to rbox
310 !
311 rbox(1:ltot2) = max( 0.0, rsub(1:ltot2,kpeg,mpeg) )
312 rbox_sv1(1:ltot2) = rbox(1:ltot2)
314 !
315 !
316 ! do bulk cloud-water chemistry
317 !
318 !
319 icase_in = icase
320 iradical_in = 1
321 idecomp_hmsa_hso5 = 1
323 co2_mixrat_in = 350.0
325 photol_no2_in = photol_no2_box
327 ! when xprescribe_ph >= 0, ph is forced to its value
328 xprescribe_ph = -1.0e31
330 ! turn off aqueous phase photolytic and radical chemistry
331 ! if either of the iphotol_onoff and iradical_onoff flags are 0
332 if ((iphotol_onoff .le. 0) .or. (iradical_onoff .le. 0)) then
333 photol_no2_in = 0.0
334 iradical_in = 0
335 end if
337 #if defined ( ccboxtest_box_testing_active)
338 ! following is for off-line box testing only
339 call ccboxtest_extra_args_aa( 'get', &
340 co2_mixrat_in, iradical_in, &
341 idecomp_hmsa_hso5, icase_in, &
342 xprescribe_ph )
343 #endif
345 gas_aqfrac_box(:) = 0.0
348 ! make call to interface_to_aqoperator1
349 call interface_to_aqoperator1( &
350 istat_aqop, &
351 dtstepc, &
352 rbox, gas_aqfrac_box, &
353 rbulk_cwaer, lptr_yyy_cwaer, &
354 co2_mixrat_in, photol_no2_in, xprescribe_ph, &
355 iradical_in, idecomp_hmsa_hso5, &
356 yaq_beg, yaq_end, ph_cmuaq_cur, &
357 numgas_aqfrac, id, it, jt, kt, kpeg, mpeg, ktau, icase_in )
359 ph_aq_box = ph_cmuaq_cur
362 #if defined ( ccboxtest_box_testing_active)
363 ! following is for off-line box testing only
364 call ccboxtest_extra_args_bb( 'put', &
365 yaq_beg, yaq_end, ph_cmuaq_cur )
366 #endif
369 !
370 !
371 ! calculate fraction of cloud-water associated with each activated aerosol bin
372 !
373 !
374 call partition_cldwtr( &
375 rbox, fr_partit_cw, &
376 it, jt, kt, kpeg, mpeg, icase_in )
378 !
379 !
380 ! distribute changes in bulk cloud-water composition among size bins
381 !
382 !
383 call distribute_bulk_changes( &
384 rbox, rbox_sv1, fr_partit_cw, &
385 rbulk_cwaer, lptr_yyy_cwaer, &
386 it, jt, kt, kpeg, mpeg, icase_in )
389 !
390 ! xfer back to rsub
391 !
392 rsub(1:ltot2,kpeg,mpeg) = max( 0.0, rbox(1:ltot2) )
395 !
396 ! do move-sections
397 !
398 if (msectional .lt. 1000000000) then
399 call cloudchem_apply_move_sections( &
400 rbox, rbox_sv1, &
401 it, jt, kt, kpeg, mpeg, icase_in )
402 end if
406 return
407 end subroutine mosaic_cloudchem_1box
411 !-----------------------------------------------------------------------
412 subroutine interface_to_aqoperator1( &
413 istat_aqop, &
414 dtstepc, &
415 rbox, gas_aqfrac_box, &
416 rbulk_cwaer, lptr_yyy_cwaer, &
417 co2_mixrat_in, photol_no2_in, xprescribe_ph, &
418 iradical_in, idecomp_hmsa_hso5, &
419 yaq_beg, yaq_end, ph_cmuaq_cur, &
420 numgas_aqfrac, id, it, jt, kt, kpeg, mpeg, ktau, icase )
422 use module_state_description, only: &
423 num_chem, param_first_scalar, p_qc, &
424 p_nh3, p_hno3, p_hcl, p_sulf, p_h2so4, p_hcho, &
425 p_ora1, p_so2, p_h2o2, p_o3, p_ho, &
426 p_ho2, p_no3, p_no, p_no2, p_hono, &
427 ! p_pan, p_ch3o2, p_ch3oh, p_op1
428 p_pan, p_ch3o2, p_ch3oh, p_op1, &
429 p_hcooh, p_ch3ooh
431 use module_data_cmu_bulkaqchem, only: &
432 meqn1max, naers, ngas, &
433 na4, naa, nac, nae, nah, nahmsa, nahso5, &
434 nan, nao, nar, nas, naw, &
435 ng4, nga, ngc, ngch3co3h, ngch3o2, ngch3o2h, ngch3oh, &
436 ngh2o2, nghcho, nghcooh, nghno2, ngho2, &
437 ngn, ngno, ngno2, ngno3, ngo3, ngoh, ngpan, ngso2
439 use module_cmu_bulkaqchem, only: aqoperator1
441 use module_data_mosaic_asect, only: &
442 maxd_asize, maxd_atype, &
443 cw_phase, nsize_aer, ntype_aer, &
444 lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
445 lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
446 lptr_oin_aer, lptr_bc_aer, lptr_oc_aer, &
447 mw_cl_aer, mw_na_aer, mw_nh4_aer, mw_no3_aer, mw_so4_aer
449 use module_data_mosaic_other, only: &
450 aboxtest_units_convert, cairclm, &
451 ktemp, l2maxd, ptotclm, rcldwtr_sub
454 implicit none
456 ! subr arguments
457 integer, intent(in) :: &
458 iradical_in, idecomp_hmsa_hso5, &
459 numgas_aqfrac, id, it, jt, kt, kpeg, mpeg, ktau, icase
460 integer, intent(inout) :: &
461 istat_aqop
463 integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
465 real, intent(in) :: &
466 dtstepc, co2_mixrat_in, &
467 photol_no2_in, xprescribe_ph
469 real, intent(inout) :: ph_cmuaq_cur
471 real, intent(inout), dimension( 1:l2maxd ) :: rbox
473 real, intent(inout), dimension( nyyy, 2 ) :: rbulk_cwaer
475 real, intent(inout), dimension( 1:numgas_aqfrac ) :: gas_aqfrac_box
477 real, intent(inout), dimension( meqn1max ) :: yaq_beg, yaq_end
480 ! local variables
481 integer :: i, iphase, isize, itype
482 integer :: iaq, istat_fatal, istat_warn
483 integer :: l, lunxx, lyyy
484 integer :: p1st
486 real, parameter :: eps=0.622 ! (mw h2o)/(mw air)
489 real :: cair_moleperm3
490 real :: dum, dumb
491 real :: factgas, factlwc, factpatm, factphoto
492 real :: factaerbc, factaercl, factaerna, factaernh4, &
493 factaerno3, factaeroc, factaeroin, factaerso4
494 real :: lwc
495 real :: p_atm, photo_in
496 real :: rh
497 real :: temp, tstep_beg_sec, tstep_end_sec
498 real :: totsulf_beg, totsulf_end
499 real :: gas(ngas), aerosol(naers)
500 real :: gas_aqfrac_cmu(ngas)
502 double precision tstep_beg_sec_dp, tstep_end_sec_dp, &
503 temp_dp, p_atm_dp, lwc_dp, rh_dp, &
504 co2_mixrat_in_dp, photo_in_dp, ph_cmuaq_cur_dp, &
505 xprescribe_ph_dp
506 double precision gas_dp(ngas), gas_aqfrac_cmu_dp(ngas), &
507 aerosol_dp(naers), yaq_beg_dp(meqn1max), yaq_end_dp(meqn1max)
511 p1st = param_first_scalar
513 !
514 ! units conversion factors
515 ! 'cmuaq-bulk' value = pegasus value X factor
516 !
517 ! [pres in atmospheres] = [pres in dynes/cm2] * factpatm
518 factpatm = 1.0/1.01325e6
519 ! [cldwtr in g-h2o/m3-air] = [cldwtr in mole-h2o/mole-air] * factlwc
520 factlwc = 28.966*eps*1.0e6*cairclm(kpeg)
521 ! [aq photolysis rate scaling factor in --] = [jno2 in 1/min] * factphoto
522 factphoto = 1.6
524 ! [gas in ppm] = [gas in mole/mole-air] * factgas
525 factgas = 1.0e6
527 ! [aerosol in ug/m3-air] = [aerosol in mole/mole-air] * factaer
528 dum = cairclm(kpeg)*1.0e12
529 factaerso4 = dum*mw_so4_aer
530 factaerno3 = dum*mw_no3_aer
531 factaercl = dum*mw_cl_aer
532 factaernh4 = dum*mw_nh4_aer
533 factaerna = dum*mw_na_aer
534 factaeroin = dum
535 factaeroc = dum
536 factaerbc = dum
538 ! rce 2005-jul-11 - use same molecular weights here as in cmu code
539 ! factaerso4 = dum*96.0
540 ! factaerno3 = dum*62.0
541 ! factaercl = dum*35.5
542 ! factaernh4 = dum*18.0
543 ! factaerna = dum*23.0
545 ! If aboxtest_units_convert=10, turn off units conversions both here
546 ! and in module_mosaic. This is for testing, to allow exact agreements.
547 if (aboxtest_units_convert .eq. 10) then
548 factpatm = 1.0
549 factlwc = 1.0
550 factphoto = 1.0
551 factgas = 1.0
552 factaerso4 = 1.0
553 factaerno3 = 1.0
554 factaercl = 1.0
555 factaernh4 = 1.0
556 factaerna = 1.0
557 factaeroin = 1.0
558 factaeroc = 1.0
559 factaerbc = 1.0
560 end if
562 !
563 ! map from rbox to gas,aerosol
564 !
565 temp = rbox(ktemp)
567 lwc = rcldwtr_sub(kpeg,mpeg) * factlwc
568 p_atm = ptotclm(kpeg) * factpatm
570 ! rce 2005-jul-11 - set p_atm so that cmu code's cair will match cairclm
571 p_atm = cairclm(kpeg)*1.0e3*0.082058e0*temp
573 photo_in = photol_no2_in * factphoto
575 rh = 1.0
576 iaq = 1
578 tstep_beg_sec = 0.0
579 tstep_end_sec = dtstepc
581 ! map gases and convert to ppm
582 gas(:) = 0.0
584 gas(nga ) = rbox(p_nh3 ) * factgas
585 gas(ngn ) = rbox(p_hno3 ) * factgas
586 if(p_hcl > param_first_scalar ) gas(ngc ) = rbox(p_hcl ) * factgas
587 if(p_sulf > param_first_scalar ) gas(ng4 ) = rbox(p_sulf ) * factgas
588 if(p_h2so4 > param_first_scalar ) gas(ng4 ) = rbox(p_h2so4 ) * factgas
591 gas(nghcho ) = rbox(p_hcho ) * factgas
592 if(p_ora1 > param_first_scalar ) gas(nghcooh ) = rbox(p_ora1 ) * factgas
593 if(p_hcooh > param_first_scalar ) gas(nghcooh ) = rbox(p_hcooh ) * factgas
594 gas(ngso2 ) = rbox(p_so2 ) * factgas
595 gas(ngh2o2 ) = rbox(p_h2o2 ) * factgas
596 gas(ngo3 ) = rbox(p_o3 ) * factgas
597 gas(ngoh ) = rbox(p_ho ) * factgas
598 gas(ngho2 ) = rbox(p_ho2 ) * factgas
599 gas(ngno3 ) = rbox(p_no3 ) * factgas
601 gas(ngno ) = rbox(p_no ) * factgas
602 gas(ngno2 ) = rbox(p_no2 ) * factgas
603 gas(nghno2 ) = rbox(p_hono ) * factgas
604 gas(ngpan ) = rbox(p_pan ) * factgas
605 gas(ngch3o2 ) = rbox(p_ch3o2 ) * factgas
606 gas(ngch3oh ) = rbox(p_ch3oh ) * factgas
607 if(p_op1 > param_first_scalar ) gas(ngch3o2h) = rbox(p_op1 ) * factgas
608 if(p_ch3ooh > param_first_scalar ) gas(ngch3o2h) = rbox(p_ch3ooh) * factgas
610 ! compute bulk activated-aerosol mixing ratios
611 aerosol(:) = 0.0
612 rbulk_cwaer(:,:) = 0.0
614 iphase = cw_phase
615 do itype = 1, ntype_aer
616 do isize = 1, nsize_aer(itype)
618 do lyyy = 1, nyyy
620 l = lptr_yyy_cwaer(isize,itype,lyyy)
621 if (l .ge. p1st) rbulk_cwaer(lyyy,1) = rbulk_cwaer(lyyy,1) + rbox(l)
623 end do
625 end do
626 end do
628 ! map them to 'aerosol' array and convert to ug/m3
629 aerosol(na4) = rbulk_cwaer(l_so4_aqyy,1) * factaerso4
630 aerosol(nan) = rbulk_cwaer(l_no3_aqyy,1) * factaerno3
631 aerosol(nac) = rbulk_cwaer(l_cl_aqyy, 1) * factaercl
632 aerosol(naa) = rbulk_cwaer(l_nh4_aqyy,1) * factaernh4
633 aerosol(nas) = rbulk_cwaer(l_na_aqyy, 1) * factaerna
634 aerosol(nar) = rbulk_cwaer(l_oin_aqyy,1) * factaeroin
635 aerosol(nae) = rbulk_cwaer(l_bc_aqyy, 1) * factaerbc
636 aerosol(nao) = rbulk_cwaer(l_oc_aqyy, 1) * factaeroc
639 !
640 ! make call to aqoperator1
641 !
642 !#if defined ( ccboxtest_box_testing_active)
643 #if defined ( temp_box_testing_active)
644 lunxx = 87
645 lunxx = -1
646 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
647 if (lunxx .gt. 0) then
648 write(lunxx,*)
649 write(lunxx,*)
650 write(lunxx,*) 'interface_to_aqoperator1 - icase, irad, idecomp'
651 write(lunxx,9870) icase, iradical_in, idecomp_hmsa_hso5
652 write(lunxx,*) 'it, jt, kt, kpeg, mpeg, ktau'
653 write(lunxx,9870) it, jt, kt, kpeg, mpeg, ktau
654 write(lunxx,*) 'temp, p_atm, lwc, photo, co2, xprescribe_ph'
655 write(lunxx,9875) temp, p_atm, lwc, photo_in, co2_mixrat_in, xprescribe_ph
656 write(lunxx,*) 'ptot, cair, rcldwtr_clm, dt_sec'
657 write(lunxx,9875) ptotclm(kpeg), cairclm(kpeg), &
658 rcldwtr_sub(kpeg,mpeg), (tstep_end_sec-tstep_beg_sec)
659 write(lunxx,*) 'gas (1=nh3, 2=hno3, 3=hcl, 4=h2so4, 11=so2, 12=h2o2, 18=o3)'
660 write(lunxx,9875) gas
661 write(lunxx,*) 'rbox(nh3, hno3, hcl, h2so4, so2, h2o2, o3)'
662 write(lunxx,9875) rbox(p_nh3), rbox(p_hno3), rbox(p_hcl), &
663 rbox(p_sulf),rbox(p_h2so4), rbox(p_so2), rbox(p_h2o2), rbox(p_o3)
664 write(lunxx,*) 'aerosol (1=na, 3=nh4, 4=no3, 5=cl, 6=so4, 8=ec, 9=oc, 10=crus)'
665 write(lunxx,9875) aerosol
666 write(lunxx,*) 'rbulk_cwaer (1=so4, 2=no3, 3-cl, 4=nh4, 5=na, 6=oin, 7=bc, 8=oc)'
667 write(lunxx,9875) rbulk_cwaer(:,1)
668 ! if (icase .ge. 3) then
669 ! write(*,*) &
670 ! '*** stopping in interface_to_aqop1 at icase =', icase
671 ! stop
672 ! end if
673 end if
674 9870 format( 8i5 )
675 9875 format( 5(1pe14.6) )
676 #endif
678 #if 0
679 ! Print outs for debugging of aqoperator1... wig, 26-Oct-2005
680 !!$ if( (id == 1 .and. ktau >= 207 ) .or. &
681 !!$ (id == 2 .and. ktau >= 610 ) .or. &
682 !!$ (id == 3 .and. ktau >= 1830 ) ) then
683 write(6,'(a)') '---Begin input for aqoperator1---'
684 write(6,'(a,4i)') 'id, it, jt, kt =', id, it, jt, kt
685 write(6,'(a,1p,2e20.12)') 'tstep_beg_sec, tstep_end_sec = ',tstep_beg_sec, tstep_end_sec
686 do l=1,ngas
687 write(6,'("gas(",i2,") = ",1p,1e20.12)') l, gas(l)
688 end do
689 do l=1,naers
690 write(6,'("aerosol(",i2,") = ",1p,1e20.12)') l, aerosol(l)
691 end do
692 write(6,'(a,1p,4e20.12)') "temp, p_atm, lwc, rh = ", temp, p_atm, lwc, rh
693 write(6,'(a,1p,3e20.12)') "co2_mixrat_in, photo_in, xprescribe_ph = ", &
694 co2_mixrat_in, photo_in, xprescribe_ph
695 write(6,'(a,3i)') " iradical_in, idecomp_hmsa_hso5, iaq = ", iradical_in, idecomp_hmsa_hso5, iaq
696 write(6,'(a)') "---End input for aqoperator1---"
697 !!$ end if
698 #endif
701 ! convert arguments to double prec
702 tstep_beg_sec_dp = 0.0d0
703 if (tstep_beg_sec .ne. 0.0) tstep_beg_sec_dp = tstep_beg_sec
704 tstep_end_sec_dp = 0.0d0
705 if (tstep_end_sec .ne. 0.0) tstep_end_sec_dp = tstep_end_sec
706 temp_dp = 0.0d0
707 if (temp .ne. 0.0) temp_dp = temp
708 p_atm_dp = 0.0d0
709 if (p_atm .ne. 0.0) p_atm_dp = p_atm
710 lwc_dp = 0.0d0
711 if (lwc .ne. 0.0) lwc_dp = lwc
712 rh_dp = 0.0d0
713 if (rh .ne. 0.0) rh_dp = rh
714 co2_mixrat_in_dp = 0.0d0
715 if (co2_mixrat_in .ne. 0.0) co2_mixrat_in_dp = co2_mixrat_in
716 photo_in_dp = 0.0d0
717 if (photo_in .ne. 0.0) photo_in_dp = photo_in
718 xprescribe_ph_dp = 0.0d0
719 if (xprescribe_ph .ne. 0.0) xprescribe_ph_dp = xprescribe_ph
720 ph_cmuaq_cur_dp = 0.0d0
721 if (ph_cmuaq_cur .ne. 0.0) ph_cmuaq_cur_dp = ph_cmuaq_cur
723 do i = 1, ngas
724 gas_dp(i) = 0.0d0
725 if (gas(i) .ne. 0.0) gas_dp(i) = gas(i)
726 end do
727 do i = 1, naers
728 aerosol_dp(i) = 0.0d0
729 if (aerosol(i) .ne. 0.0) aerosol_dp(i) = aerosol(i)
730 end do
731 do i = 1, ngas
732 gas_aqfrac_cmu_dp(i) = 0.0d0
733 if (gas_aqfrac_cmu(i) .ne. 0.0) gas_aqfrac_cmu_dp(i) = gas_aqfrac_cmu(i)
734 end do
735 do i = 1, meqn1max
736 yaq_beg_dp(i) = 0.0d0
737 if (yaq_beg(i) .ne. 0.0) yaq_beg_dp(i) = yaq_beg(i)
738 end do
739 do i = 1, meqn1max
740 yaq_end_dp(i) = 0.0d0
741 if (yaq_end(i) .ne. 0.0) yaq_end_dp(i) = yaq_end(i)
742 end do
745 ! total sulfur species conc as sulfate (ug/m3)
746 cair_moleperm3 = 1.0e3*p_atm_dp/(0.082058e0*temp_dp)
747 totsulf_beg = ( aerosol_dp(na4)/96. &
748 + aerosol_dp(nahso5)/113. + aerosol_dp(nahmsa)/111. &
749 + (gas_dp(ngso2) + gas_dp(ng4))*cair_moleperm3 )*96.0
751 ! call aqoperator1( &
752 ! istat_fatal, istat_warn, &
753 ! tstep_beg_sec, tstep_end_sec, &
754 ! gas, aerosol, gas_aqfrac_cmu, &
755 ! temp, p_atm, lwc, rh, &
756 ! co2_mixrat_in, photo_in, xprescribe_ph, &
757 ! iradical_in, idecomp_hmsa_hso5, iaq, &
758 ! yaq_beg, yaq_end, ph_cmuaq_cur )
760 call aqoperator1( &
761 istat_fatal, istat_warn, &
762 tstep_beg_sec_dp, tstep_end_sec_dp, &
763 gas_dp, aerosol_dp, gas_aqfrac_cmu_dp, &
764 temp_dp, p_atm_dp, lwc_dp, rh_dp, &
765 co2_mixrat_in_dp, photo_in_dp, xprescribe_ph_dp, &
766 iradical_in, idecomp_hmsa_hso5, iaq, &
767 yaq_beg_dp, yaq_end_dp, ph_cmuaq_cur_dp )
769 totsulf_end = ( aerosol_dp(na4)/96. &
770 + aerosol_dp(nahso5)/113. + aerosol_dp(nahmsa)/111. &
771 + (gas_dp(ngso2) + gas_dp(ng4))*cair_moleperm3 )*96.0
774 ! convert arguments back to single prec
775 tstep_beg_sec = tstep_beg_sec_dp
776 tstep_end_sec = tstep_end_sec_dp
777 temp = temp_dp
778 p_atm = p_atm_dp
779 lwc = lwc_dp
780 rh = rh_dp
781 ! co2_mixrat_in = co2_mixrat_in_dp ! this has intent(in)
782 ! photo_in = photo_in_dp ! this has intent(in)
783 ! xprescribe_ph = xprescribe_ph_dp ! this has intent(in)
784 ph_cmuaq_cur = ph_cmuaq_cur_dp
786 do i = 1, ngas
787 gas(i) = gas_dp(i)
788 end do
789 do i = 1, naers
790 aerosol(i) = aerosol_dp(i)
791 end do
792 do i = 1, ngas
793 gas_aqfrac_cmu(i) = gas_aqfrac_cmu_dp(i)
794 end do
795 do i = 1, meqn1max
796 yaq_beg(i) = yaq_beg_dp(i)
797 end do
798 do i = 1, meqn1max
799 yaq_end(i) = yaq_end_dp(i)
800 end do
803 !
804 ! warning message when status flags are non-zero
805 !
806 istat_aqop = 0
807 if (istat_fatal .ne. 0) then
808 write(6,*) &
809 '*** mosaic_cloudchem_driver, subr interface_to_aqoperator1'
810 write(6,'(a,4i5,2i10)') &
811 ' id,it,jt,kt, istat_fatal, warn =', &
812 id, it, jt, kt, istat_fatal, istat_warn
813 istat_aqop = -10
814 end if
816 !
817 ! warning message when sulfur mass balance error exceeds the greater
818 ! of (1.0e-3 ug/m3) OR (1.0e-3 X total sulfur mixing ratio)
819 !
820 dum = totsulf_end - totsulf_beg
821 dumb = max( totsulf_beg, totsulf_end )
822 if (abs(dum) .gt. max(1.0e-3,1.0e-3*dumb)) then
823 write(6,*) &
824 '*** mosaic_cloudchem_driver, sulfur balance warning'
825 write(6,'(a,4i5,1p,3e12.4)') &
826 ' id,it,jt,kt, total_sulfur_beg, _end, _error =', &
827 id, it, jt, kt, totsulf_beg, totsulf_end, dum
828 end if
830 !
831 ! map from gas,aerosol to rbox
832 !
833 rbox(p_nh3 ) = gas(nga ) / factgas
834 rbox(p_hno3 ) = gas(ngn ) / factgas
835 if(p_hcl .gt. param_first_scalar) rbox(p_hcl ) = gas(ngc ) / factgas
836 if(p_sulf .gt. param_first_scalar) rbox(p_sulf ) = gas(ng4) / factgas
837 if(p_h2so4 .gt. param_first_scalar) rbox(p_h2so4 ) = gas(ng4) / factgas
840 rbox(p_hcho ) = gas(nghcho ) / factgas
841 if(p_ora1 > param_first_scalar ) rbox(p_ora1 ) = gas(nghcooh ) / factgas
842 if(p_hcooh > param_first_scalar ) rbox(p_hcooh) = gas(nghcooh ) / factgas
843 rbox(p_so2 ) = gas(ngso2 ) / factgas
844 rbox(p_h2o2 ) = gas(ngh2o2 ) / factgas
845 rbox(p_o3 ) = gas(ngo3 ) / factgas
846 rbox(p_ho ) = gas(ngoh ) / factgas
847 rbox(p_ho2 ) = gas(ngho2 ) / factgas
848 rbox(p_no3 ) = gas(ngno3 ) / factgas
850 rbox(p_no ) = gas(ngno ) / factgas
851 rbox(p_no2 ) = gas(ngno2 ) / factgas
852 rbox(p_hono ) = gas(nghno2 ) / factgas
853 rbox(p_pan ) = gas(ngpan ) / factgas
854 rbox(p_ch3o2 ) = gas(ngch3o2 ) / factgas
855 rbox(p_ch3oh ) = gas(ngch3oh ) / factgas
856 if(p_op1 > param_first_scalar ) rbox(p_op1 ) = gas(ngch3o2h) / factgas
857 if(p_ch3ooh > param_first_scalar ) rbox(p_ch3ooh) = gas(ngch3o2h) / factgas
859 gas_aqfrac_box(:) = 0.0
861 if (p_nh3 .le. numgas_aqfrac) &
862 gas_aqfrac_box(p_nh3 ) = gas_aqfrac_cmu(nga )
863 if (p_hno3 .le. numgas_aqfrac) &
864 gas_aqfrac_box(p_hno3 ) = gas_aqfrac_cmu(ngn )
865 if (p_hcl .le. numgas_aqfrac .and. p_hcl .gt. param_first_scalar) &
866 gas_aqfrac_box(p_hcl ) = gas_aqfrac_cmu(ngc )
867 if (p_sulf .le. numgas_aqfrac .and. p_sulf .gt. param_first_scalar) &
868 gas_aqfrac_box(p_sulf ) = gas_aqfrac_cmu(ng4 )
869 if (p_h2so4 .le. numgas_aqfrac .and. p_h2so4 .gt. param_first_scalar) &
870 gas_aqfrac_box(p_h2so4 ) = gas_aqfrac_cmu(ng4 )
872 if (p_hcho .le. numgas_aqfrac) &
873 gas_aqfrac_box(p_hcho ) = gas_aqfrac_cmu(nghcho )
874 if (p_ora1 .le. numgas_aqfrac .and. p_ora1 .gt. param_first_scalar) &
875 gas_aqfrac_box(p_ora1 ) = gas_aqfrac_cmu(nghcooh )
876 if (p_hcooh .le. numgas_aqfrac .and. p_hcooh .gt. param_first_scalar) &
877 gas_aqfrac_box(p_hcooh ) = gas_aqfrac_cmu(nghcooh )
878 if (p_so2 .le. numgas_aqfrac) &
879 gas_aqfrac_box(p_so2 ) = gas_aqfrac_cmu(ngso2 )
880 if (p_h2o2 .le. numgas_aqfrac) &
881 gas_aqfrac_box(p_h2o2 ) = gas_aqfrac_cmu(ngh2o2 )
882 if (p_o3 .le. numgas_aqfrac) &
883 gas_aqfrac_box(p_o3 ) = gas_aqfrac_cmu(ngo3 )
884 if (p_ho .le. numgas_aqfrac) &
885 gas_aqfrac_box(p_ho ) = gas_aqfrac_cmu(ngoh )
886 if (p_ho2 .le. numgas_aqfrac) &
887 gas_aqfrac_box(p_ho2 ) = gas_aqfrac_cmu(ngho2 )
888 if (p_no3 .le. numgas_aqfrac) &
889 gas_aqfrac_box(p_no3 ) = gas_aqfrac_cmu(ngno3 )
891 if (p_no .le. numgas_aqfrac) &
892 gas_aqfrac_box(p_no ) = gas_aqfrac_cmu(ngno )
893 if (p_no2 .le. numgas_aqfrac) &
894 gas_aqfrac_box(p_no2 ) = gas_aqfrac_cmu(ngno2 )
895 if (p_hono .le. numgas_aqfrac) &
896 gas_aqfrac_box(p_hono ) = gas_aqfrac_cmu(nghno2 )
897 if (p_pan .le. numgas_aqfrac) &
898 gas_aqfrac_box(p_pan ) = gas_aqfrac_cmu(ngpan )
899 if (p_ch3o2 .le. numgas_aqfrac) &
900 gas_aqfrac_box(p_ch3o2 ) = gas_aqfrac_cmu(ngch3o2 )
901 if (p_ch3oh .le. numgas_aqfrac) &
902 gas_aqfrac_box(p_ch3oh ) = gas_aqfrac_cmu(ngch3oh )
903 if (p_op1 .le. numgas_aqfrac .and. p_op1 .gt. param_first_scalar) &
904 gas_aqfrac_box(p_op1 ) = gas_aqfrac_cmu(ngch3o2h)
905 if (p_ch3ooh.le. numgas_aqfrac .and. p_ch3ooh .gt. param_first_scalar) &
906 gas_aqfrac_box(p_ch3ooh) = gas_aqfrac_cmu(ngch3o2h)
908 rbulk_cwaer(l_so4_aqyy,2) = aerosol(na4) / factaerso4
909 rbulk_cwaer(l_no3_aqyy,2) = aerosol(nan) / factaerno3
910 rbulk_cwaer(l_cl_aqyy, 2) = aerosol(nac) / factaercl
911 rbulk_cwaer(l_nh4_aqyy,2) = aerosol(naa) / factaernh4
912 rbulk_cwaer(l_na_aqyy, 2) = aerosol(nas) / factaerna
913 rbulk_cwaer(l_oin_aqyy,2) = aerosol(nar) / factaeroin
914 rbulk_cwaer(l_bc_aqyy, 2) = aerosol(nae) / factaerbc
915 rbulk_cwaer(l_oc_aqyy, 2) = aerosol(nao) / factaeroc
918 !#if defined ( ccboxtest_box_testing_active)
919 #if defined ( temp_box_testing_active)
920 lunxx = 87
921 lunxx = -1
922 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
923 if (lunxx .gt. 0) then
924 write(lunxx,*)
925 write(lunxx,*) 'interface_to_aqoperator1 - after call'
926 write(lunxx,*) 'gas (1=nh3, 2=hno3, 3=hcl, 4=h2so4, 11=so2, 12=h2o2, 18=o3)'
927 write(lunxx,9875) gas
928 write(lunxx,*) 'rbox(nh3, hno3, hcl, h2so4, so2, h2o2, o3)'
929 write(lunxx,9875) rbox(p_nh3), rbox(p_hno3), rbox(p_hcl), &
930 rbox(p_sulf), rbox(p_h2so4), rbox(p_so2), rbox(p_h2o2), rbox(p_o3)
931 write(lunxx,*) 'aerosol (1=na, 3=nh4, 4=no3, 5=cl, 6=so4, 8=ec, 9=oc, 10=crus)'
932 write(lunxx,9875) aerosol
933 write(lunxx,*) 'rbulk_cwaer (1=so4, 2=no3, 3-cl, 4=nh4, 5=na, 6=oin, 7=bc, 8=oc)'
934 write(lunxx,9875) rbulk_cwaer(:,2)
935 write(lunxx,*) 'ph_cmuaq_cur'
936 write(lunxx,9875) ph_cmuaq_cur
937 if (icase .ge. 10) then
938 write(*,*) &
939 '*** stopping in interface_to_aqop1 at icase =', icase
940 ! stop
941 end if
942 end if
943 #endif
946 return
947 end subroutine interface_to_aqoperator1
951 !-----------------------------------------------------------------------
952 subroutine partition_cldwtr( &
953 rbox, fr_partit_cw, &
954 it, jt, kt, kpeg, mpeg, icase )
956 use module_state_description, only: &
957 param_first_scalar
959 use module_data_mosaic_asect, only: &
960 maxd_asize, maxd_atype, &
961 cw_phase, nsize_aer, ntype_aer, ncomp_aer, &
962 massptr_aer, numptr_aer, &
963 dens_aer, mw_aer, volumlo_sect, volumhi_sect
965 use module_data_mosaic_other, only: &
966 aboxtest_units_convert, cairclm, &
967 ktemp, l2maxd, ptotclm, rcldwtr_sub
970 implicit none
972 ! subr arguments
973 integer, intent(in) :: it, jt, kt, kpeg, mpeg, icase
975 real, intent(inout), dimension( 1:l2maxd ) :: rbox
977 real, intent(inout), dimension( maxd_asize, maxd_atype ) :: &
978 fr_partit_cw
980 ! local variables
981 integer :: iphase, isize, itype
982 integer :: jdone_mass, jdone_numb, jpos, jpos_mass, jpos_numb
983 integer :: l, ll, lunxx
984 integer :: p1st
986 real, parameter :: partit_wght_mass = 0.5
988 real :: dum, duma, dumb, dumc, dummass, dumnumb, dumvolu
989 real :: tmass, tnumb, umass, unumb, wmass, wnumb
990 real, dimension( maxd_asize, maxd_atype ) :: fmass, fnumb, xmass, xnumb
992 p1st = PARAM_FIRST_SCALAR
994 iphase = cw_phase
995 tmass = 0.0
996 tnumb = 0.0
997 umass = 0.0
998 unumb = 0.0
1000 ! compute
1001 ! xmass, xnumb = mass, number mixing ratio for a bin
1002 ! tmass, tnumb = sum over all bins of xmass, xnumb
1003 ! umass, unumb = max over all bins of xmass, xnumb
1004 ! set xmass, xnumb = 0.0 if bin mass, numb < 1.0e-37
1005 ! constrain xnumb so that mean particle volume is
1006 ! within bin boundaries
1007 do itype = 1, ntype_aer
1008 do isize = 1, nsize_aer(itype)
1009 dummass = 0.0
1010 dumvolu = 0.0
1011 do ll = 1, ncomp_aer(itype)
1012 l = massptr_aer(ll,isize,itype,iphase)
1013 if (l .ge. p1st) then
1014 dum = max( 0.0, rbox(l) )*mw_aer(ll,itype)
1015 dummass = dummass + dum
1016 dumvolu = dumvolu + dum/dens_aer(ll,itype)
1017 end if
1018 end do
1020 l = numptr_aer(isize,itype,iphase)
1021 dumnumb = max( 0.0, rbox(l) )
1022 if (dumnumb .gt. dumvolu/volumlo_sect(isize,itype)) then
1023 dumnumb = dumvolu/volumlo_sect(isize,itype)
1024 rbox(l) = dumnumb
1025 else if (dumnumb .lt. dumvolu/volumhi_sect(isize,itype)) then
1026 dumnumb = dumvolu/volumhi_sect(isize,itype)
1027 rbox(l) = dumnumb
1028 end if
1030 if (dummass .lt. 1.0e-37) dummass = 0.0
1031 xmass(isize,itype) = dummass
1032 if (dumnumb .lt. 1.0e-37) dumnumb = 0.0
1033 xnumb(isize,itype) = dumnumb
1035 tmass = tmass + xmass(isize,itype)
1036 tnumb = tnumb + xnumb(isize,itype)
1037 umass = max( umass, xmass(isize,itype) )
1038 unumb = max( unumb, xnumb(isize,itype) )
1039 end do
1040 end do
1042 ! compute
1043 ! fmass, fnumb = fraction of total mass, number that is in a bin
1044 ! if tmass<1e-35 and umass>0, set fmass=1 for bin with largest xmass
1045 ! if tmass<1e-35 and umass=0, set fmass=0 for all
1046 jdone_mass = 0
1047 jdone_numb = 0
1048 jpos_mass = 0
1049 jpos_numb = 0
1050 do itype = 1, ntype_aer
1051 do isize = 1, nsize_aer(itype)
1052 fmass(isize,itype) = 0.0
1053 if (tmass .ge. 1.0e-35) then
1054 fmass(isize,itype) = xmass(isize,itype)/tmass
1055 else if (umass .gt. 0.0) then
1056 if ( (jdone_mass .eq. 0) .and. &
1057 (xmass(isize,itype) .eq. umass) ) then
1058 jdone_mass = 1
1059 fmass(isize,itype) = 1.0
1060 end if
1061 end if
1062 if (fmass(isize,itype) .gt. 0) jpos_mass = jpos_mass + 1
1064 fnumb(isize,itype) = 0.0
1065 if (tnumb .ge. 1.0e-35) then
1066 fnumb(isize,itype) = xnumb(isize,itype)/tnumb
1067 else if (unumb .gt. 0.0) then
1068 if ( (jdone_numb .eq. 0) .and. &
1069 (xnumb(isize,itype) .eq. unumb) ) then
1070 jdone_numb = 1
1071 fnumb(isize,itype) = 1.0
1072 end if
1073 end if
1074 if (fnumb(isize,itype) .gt. 0) jpos_numb = jpos_numb + 1
1075 end do
1076 end do
1078 ! if only 1 bin has fmass or fnumb > 0, set value to 1.0 exactly
1079 if ((jpos_mass .eq. 1) .or. (jpos_numb .eq. 1)) then
1080 do itype = 1, ntype_aer
1081 do isize = 1, nsize_aer(itype)
1082 if (jpos_mass .eq. 1) then
1083 if (fmass(isize,itype) .gt. 0) fmass(isize,itype) = 1.0
1084 end if
1085 if (jpos_numb .eq. 1) then
1086 if (fnumb(isize,itype) .gt. 0) fnumb(isize,itype) = 1.0
1087 end if
1088 end do
1089 end do
1090 end if
1092 !
1093 ! compute fr_partit_cw as weighted average of fmass & fnumb, except
1094 ! if tmass<1e-35 and umass=0, use only fnumb
1095 ! if tnumb<1e-35 and unumb=0, use only fmass
1096 ! if tmass,tnumb<1e-35 and umass,unumb=0,
1097 ! set fr_partit_cw=1 for center bin of itype=1
1098 !
1099 fr_partit_cw(:,:) = 0.0
1100 if ((jpos_mass .eq. 0) .and. (jpos_numb .eq. 0)) then
1101 itype = 1
1102 isize = (nsize_aer(itype)+1)/2
1103 fr_partit_cw(isize,itype) = 1.0
1105 else if (jpos_mass .eq. 0) then
1106 fr_partit_cw(:,:) = fnumb(:,:)
1108 else if (jpos_numb .eq. 0) then
1109 fr_partit_cw(:,:) = fmass(:,:)
1111 else
1112 wmass = max( 0.0, min( 1.0, partit_wght_mass ) )
1113 wnumb = 1.0 - wmass
1114 fr_partit_cw(:,:) = wmass*fmass(:,:) + wnumb*fnumb(:,:)
1116 jpos = 0
1117 do itype = 1, ntype_aer
1118 do isize = 1, nsize_aer(itype)
1119 if (fr_partit_cw(isize,itype) .gt. 0.0) jpos = jpos + 1
1120 end do
1121 end do
1123 ! if only 1 bin has fr_partit_cw > 0, set value to 1.0 exactly
1124 if (jpos .eq. 1) then
1125 do itype = 1, ntype_aer
1126 do isize = 1, nsize_aer(itype)
1127 if (fr_partit_cw(isize,itype) .gt. 0.0) &
1128 fr_partit_cw(isize,itype) = 1.0
1129 end do
1130 end do
1131 end if
1132 end if
1135 !#if defined ( ccboxtest_box_testing_active)
1136 #if defined ( temp_box_testing_active)
1137 ! diagnostics when lunxx > 0
1138 lunxx = 86
1139 lunxx = -1
1140 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1141 if (lunxx .gt. 0) then
1142 if (icase .le. 9) then
1143 write(lunxx,9800)
1144 write(lunxx,9800)
1145 write(lunxx,9800) &
1146 'partition_cldwtr - icase, jpos, jpos_mass, jpos_numb'
1147 write(lunxx,9810) icase, jpos, jpos_mass, jpos_numb
1148 write(lunxx,9800) 'tmass, umass, wmass'
1149 write(lunxx,9820) tmass, umass, wmass
1150 write(lunxx,9800) 'tnumb, unumb, wnumb'
1151 write(lunxx,9820) tnumb, unumb, wnumb
1152 write(lunxx,9800) 'xmass, fmass, xnumb, fnumb, fr_partit_cw'
1153 duma = 0.0
1154 dumb = 0.0
1155 dumc = 0.0
1156 do itype = 1, ntype_aer
1157 do isize = 1, nsize_aer(itype)
1158 write(lunxx,9820) xmass(isize,itype), fmass(isize,itype), &
1159 xnumb(isize,itype), fnumb(isize,itype), &
1160 fr_partit_cw(isize,itype)
1161 duma = duma + fmass(isize,itype)
1162 dumb = dumb + fnumb(isize,itype)
1163 dumc = dumc + fr_partit_cw(isize,itype)
1164 end do
1165 end do
1166 write(lunxx,9800) &
1167 'sum_fmass-1.0, sum_fnumb-1.0, sum_fr_partit-1.0'
1168 write(lunxx,9820) (duma-1.0), (dumb-1.0), (dumc-1.0)
1169 if (icase .eq. -2) then
1170 write(*,*) '*** stopping in partition_cldwtr at icase =', icase
1171 ! stop
1172 end if
1173 9800 format( a )
1174 9810 format( 5i10 )
1175 9820 format( 5(1pe10.2) )
1176 end if
1177 end if
1178 #endif
1181 return
1182 end subroutine partition_cldwtr
1186 !-----------------------------------------------------------------------
1187 subroutine distribute_bulk_changes( &
1188 rbox, rbox_sv1, fr_partit_cw, &
1189 rbulk_cwaer, lptr_yyy_cwaer, &
1190 it, jt, kt, kpeg, mpeg, icase )
1192 use module_state_description, only: &
1193 param_first_scalar
1195 use module_data_mosaic_asect, only: &
1196 maxd_asize, maxd_atype, &
1197 cw_phase, nsize_aer, ntype_aer, &
1198 lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
1199 lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
1200 lptr_oin_aer, lptr_bc_aer, lptr_oc_aer
1202 use module_data_mosaic_other, only: l2maxd, lunout, name
1205 implicit none
1207 ! subr arguments
1208 integer, intent(in) :: it, jt, kt, kpeg, mpeg, icase
1210 integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
1212 real, intent(inout), dimension( 1:l2maxd ) :: rbox, rbox_sv1
1214 real, intent(in), dimension( maxd_asize, maxd_atype ) :: &
1215 fr_partit_cw
1217 real, intent(in), dimension( nyyy, 2 ) :: rbulk_cwaer
1220 ! local variables
1221 integer :: iphase, isize, itype
1222 integer :: idone, icount, ncount
1223 integer :: jpos, jpos_sv
1224 integer :: l, lunxx, lunxxaa, lunxxbb, lyyy
1225 integer :: p1st
1227 real :: duma, dumb, dumc
1228 real :: fr, frsum_cur
1229 real :: fr_cur(maxd_asize,maxd_atype)
1230 real :: del_r_current, del_r_remain
1231 real :: del_rbulk_cwaer(nyyy)
1234 p1st = param_first_scalar
1236 do lyyy = 1, nyyy
1237 del_rbulk_cwaer(lyyy) = rbulk_cwaer(lyyy,2) - rbulk_cwaer(lyyy,1)
1238 end do
1240 iphase = cw_phase
1243 jpos = 0
1244 do itype = 1, ntype_aer
1245 do isize = 1, nsize_aer(itype)
1246 if (fr_partit_cw(isize,itype) .gt. 0) jpos = jpos + 1
1247 end do
1248 end do
1249 jpos_sv = jpos
1251 !
1252 ! distribution is trivial when only 1 bin has fr_partit_cw > 0
1253 !
1254 if (jpos_sv .eq. 1) then
1255 do lyyy = 1, nyyy
1257 do itype = 1, ntype_aer
1258 do isize = 1, nsize_aer(itype)
1259 fr = fr_partit_cw(isize,itype)
1260 if (fr .eq. 1.0) then
1261 l = lptr_yyy_cwaer(isize,itype,lyyy)
1262 if (l .ge. p1st) rbox(l) = rbulk_cwaer(lyyy,2)
1263 end if
1264 end do
1265 end do
1267 end do
1268 goto 7900
1269 end if
1272 do 3900 lyyy = 1, nyyy
1274 !
1275 ! distribution is simple when del_rbulk_cwaer(lyyy) >= 0
1276 !
1277 if (del_rbulk_cwaer(lyyy) .eq. 0.0) then
1278 goto 3900
1279 else if (del_rbulk_cwaer(lyyy) .gt. 0.0) then
1280 do itype = 1, ntype_aer
1281 do isize = 1, nsize_aer(itype)
1282 fr = fr_partit_cw(isize,itype)
1283 if (fr .gt. 0.0) then
1284 l = lptr_yyy_cwaer(isize,itype,lyyy)
1285 if (l .ge. p1st) then
1286 rbox(l) = rbox(l) + fr*del_rbulk_cwaer(lyyy)
1287 end if
1288 end if
1289 end do
1290 end do
1292 goto 3900
1293 end if
1295 !
1296 ! distribution is complicated when del_rbulk_cwaer(lyyy) < 0,
1297 ! because you cannot produce any negative mixrats
1298 !
1299 del_r_remain = del_rbulk_cwaer(lyyy)
1300 fr_cur(:,:) = fr_partit_cw(:,:)
1302 ncount = max( 1, jpos_sv*2 )
1303 icount = 0
1305 ! iteration loop
1306 do while (icount .le. ncount)
1308 icount = icount + 1
1309 del_r_current = del_r_remain
1310 jpos = 0
1311 frsum_cur = 0.0
1313 do itype = 1, ntype_aer
1314 do isize = 1, nsize_aer(itype)
1316 fr = fr_cur(isize,itype)
1318 if (fr .gt. 0.0) then
1319 l = lptr_yyy_cwaer(isize,itype,lyyy)
1320 if (l .ge. p1st) then
1321 duma = fr*del_r_current
1322 dumb = rbox(l) + duma
1323 if (dumb .gt. 0.0) then
1324 jpos = jpos + 1
1325 else if (dumb .eq. 0.0) then
1326 fr_cur(isize,itype) = 0.0
1327 else
1328 duma = -rbox(l)
1329 dumb = 0.0
1330 fr_cur(isize,itype) = 0.0
1331 end if
1332 del_r_remain = del_r_remain - duma
1333 rbox(l) = dumb
1334 frsum_cur = frsum_cur + fr_cur(isize,itype)
1335 else
1336 fr_cur(isize,itype) = 0.0
1337 end if
1338 end if
1340 end do ! isize = 1, nsize_aer
1341 end do ! itype = 1, ntype_aer
1343 ! done if jpos = jpos_sv, because bins reached zero mixrat
1344 if (jpos .eq. jpos_sv) then
1345 idone = 1
1346 ! del_r_remain starts as negative, so done if non-negative
1347 else if (del_r_remain .ge. 0.0) then
1348 idone = 2
1349 ! del_r_remain starts as negative, so done if non-negative
1350 else if (abs(del_r_remain) .le. 1.0e-7*abs(del_rbulk_cwaer(lyyy))) then
1351 idone = 3
1352 ! done if all bins have fr_cur = 0
1353 else if (frsum_cur .le. 0.0) then
1354 idone = 4
1355 ! same thing basically
1356 else if (jpos .le. 0) then
1357 idone = 5
1358 else
1359 idone = 0
1360 end if
1362 ! check for done, and (conditionally) print message
1363 if (idone .gt. 0) then
1364 lunxxaa = 6
1365 !#if defined ( ccboxtest_box_testing_active)
1366 #if defined ( temp_box_testing_active)
1367 lunxxaa = 86
1368 #endif
1369 if ((lunxxaa .gt. 0) .and. (icount .gt. (1+jpos_sv)/2)) then
1370 write(lunxxaa,9800) &
1371 'distribute_bulk_changes - icount>jpos_sv/2 - i,j,k'
1372 write(lunxxaa,9810) it, jt, kt
1373 write(lunxxaa,9800) 'icase, lyyy, idone, icount, jpos, jpos_sv'
1374 write(lunxxaa,9810) icase, lyyy, idone, icount, jpos, jpos_sv
1375 end if
1376 goto 3900
1377 end if
1379 ! rescale fr_cur for next iteration
1380 fr_cur(:,:) = fr_cur(:,:)/frsum_cur
1382 end do ! while (icount .le. ncount)
1385 ! icount > ncount, so print message
1386 lunxxbb = 6
1387 !#if defined ( ccboxtest_box_testing_active)
1388 #if defined ( temp_box_testing_active)
1389 lunxxbb = 86
1390 #endif
1391 if (lunxxbb .gt. 0) then
1392 write(lunxxbb,9800)
1393 write(lunxxbb,9800) &
1394 'distribute_bulk_changes - icount>ncount - i,j,k'
1395 write(lunxxbb,9810) it, jt, kt
1396 write(lunxxbb,9800) 'icase, lyyy, icount, ncount, jpos_sv, jpos'
1397 write(lunxxbb,9810) icase, lyyy, icount, ncount, jpos_sv, jpos
1398 write(lunxxbb,9800) 'rbulk_cwaer(1), del_rbulk_cwaer, del_r_remain, frsum_cur, (frsum_cur-1.0)'
1399 write(lunxxbb,9820) rbulk_cwaer(lyyy,1), del_rbulk_cwaer(lyyy), &
1400 del_r_remain, frsum_cur, (frsum_cur-1.0)
1401 end if
1402 9800 format( a )
1403 9801 format( 3a )
1404 9810 format( 7i10 )
1405 9820 format( 7(1pe10.2) )
1406 9840 format( 2i3, 5(1pe14.6) )
1409 3900 continue
1411 7900 continue
1414 !#if defined ( ccboxtest_box_testing_active)
1415 #if defined ( temp_box_testing_active)
1416 ! diagnostics for testing
1417 lunxx = 88
1418 lunxx = -1
1419 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1420 if (lunxx .gt. 0) then
1421 icount = 0
1422 do lyyy = 1, nyyy
1423 duma = del_rbulk_cwaer(lyyy)
1424 if ( abs(duma) .gt. &
1425 max( 1.0e-35, 1.0e-5*abs(rbulk_cwaer(lyyy,1)) ) ) then
1426 icount = icount + 1
1427 if (icount .eq. 1) write(lunxx,9800)
1428 if (icount .eq. 1) write(lunxx,9800)
1429 write(lunxx,9800)
1430 write(lunxx,9801) 'distribute_bulk_changes - ',name(lptr_yyy_cwaer(1,1,lyyy)), ' - icase, lyyy'
1431 write(lunxx,9810) icase, lyyy
1432 write(lunxx,9800) ' tp sz rbox_sv1, rbox, del_rbox, del_rbox/del_rbulk_cwaer, (...-fr_partit_cw)'
1433 write(lunxx,9840) 0, 0, rbulk_cwaer(lyyy,1:2), duma
1434 do itype = 1, ntype_aer
1435 do isize = 1, nsize_aer(itype)
1436 l = lptr_yyy_cwaer(isize,itype,lyyy)
1437 dumb = rbox(l) - rbox_sv1(l)
1438 dumc = dumb/max( abs(duma), 1.0e-35 )
1439 if (duma .lt. 0.0) dumc = -dumc
1440 write(lunxx,9840) itype, isize, rbox_sv1(l), rbox(l), &
1441 dumb, dumc, (dumc-fr_partit_cw(isize,itype))
1442 end do
1443 end do
1444 end if
1445 end do
1446 if (icase .ge. 5) then
1447 write(*,*) &
1448 '*** stop in distribute_bulk_changes diags, icase =', icase
1449 ! stop
1450 end if
1451 end if
1452 #endif
1455 return
1456 end subroutine distribute_bulk_changes
1460 !-----------------------------------------------------------------------
1461 subroutine cloudchem_apply_move_sections( &
1462 rbox, rbox_sv1, &
1463 it, jt, kt, kpeg, mpeg, icase )
1465 use module_state_description, only: &
1466 param_first_scalar
1468 use module_data_mosaic_asect, only: &
1469 msectional, &
1470 maxd_asize, maxd_atype, &
1471 cw_phase, nsize_aer, ntype_aer, ncomp_aer, &
1472 massptr_aer, numptr_aer, mw_aer, dens_aer, &
1473 lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
1474 lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
1475 lptr_oin_aer, lptr_bc_aer, lptr_oc_aer, &
1476 drymass_aftgrow, drymass_pregrow, &
1477 drydens_aftgrow, drydens_pregrow
1479 use module_data_mosaic_other, only: l2maxd, name, rsub
1481 use module_mosaic_movesect, only: move_sections
1484 implicit none
1486 ! subr arguments
1487 integer, intent(in) :: it, jt, kt, kpeg, mpeg, icase
1489 real, intent(inout), dimension( 1:l2maxd ) :: rbox, rbox_sv1
1492 ! local variables
1493 integer :: idum_msect
1494 integer :: iphase, isize, itype
1495 integer :: l, ll, lunxx
1496 integer :: p1st
1497 integer :: lptr_dum(maxd_asize,maxd_atype)
1499 real :: densdefault
1500 real :: dmaft, dmpre, dvaft, dvpre
1501 real :: duma, dumb, dumc
1502 real :: smallmassbb
1505 p1st = param_first_scalar
1506 iphase = cw_phase
1508 !
1509 ! compute drymass before and after growth
1510 ! set drydens = -1.0, so it will be calculated
1511 !
1512 densdefault = 2.0
1513 smallmassbb = 1.0e-30
1515 do 1800 itype = 1, ntype_aer
1516 do 1800 isize = 1, nsize_aer(itype)
1517 dmaft = 0.0
1518 dmpre = 0.0
1519 dvaft = 0.0
1520 dvpre = 0.0
1522 do ll = 1, ncomp_aer(itype)
1523 l = massptr_aer(ll,isize,itype,iphase)
1524 if (l .ge. p1st) then
1525 duma = mw_aer(ll,itype)
1526 dmaft = dmaft + duma*rbox(l)
1527 dmpre = dmpre + duma*rbox_sv1(l)
1529 duma = duma/dens_aer(ll,itype)
1530 dvaft = dvaft + duma*rbox(l)
1531 dvpre = dvpre + duma*rbox_sv1(l)
1532 end if
1533 end do
1535 drymass_aftgrow(isize,itype) = dmaft
1536 drymass_pregrow(isize,itype) = dmpre
1538 if (min(dmaft,dvaft) .le. smallmassbb) then
1539 drydens_aftgrow(isize,itype) = densdefault
1540 else
1541 drydens_aftgrow(isize,itype) = dmaft/dvaft
1542 end if
1543 if (min(dmpre,dvpre) .le. smallmassbb) then
1544 drydens_pregrow(isize,itype) = densdefault
1545 else
1546 drydens_pregrow(isize,itype) = dmpre/dvpre
1547 end if
1549 1800 continue
1552 ! apply move sections routine
1553 ! (and conditionally turn on movesect diagnostics)
1554 idum_msect = msectional
1556 !#if defined ( ccboxtest_box_testing_active )
1557 #if defined ( temp_box_testing_active)
1558 lunxx = 88
1559 lunxx = -1
1560 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1561 if (lunxx .gt. 0) then
1562 if (msectional .lt. 7000) msectional = msectional + 7000
1563 end if
1564 #endif
1566 call move_sections( 2, it, jt, kpeg, mpeg )
1568 msectional = idum_msect
1571 !#if defined ( ccboxtest_box_testing_active)
1572 #if defined ( temp_box_testing_active)
1573 ! diagnostics for testing
1574 lunxx = 88
1575 lunxx = -1
1576 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1577 if (lunxx .gt. 0) then
1578 do ll = 1, 4
1579 if (ll .eq. 1) then
1580 lptr_dum(:,:) = lptr_so4_aer(:,:,iphase)
1581 else if (ll .eq. 2) then
1582 lptr_dum(:,:) = lptr_nh4_aer(:,:,iphase)
1583 else if (ll .eq. 3) then
1584 lptr_dum(:,:) = lptr_oin_aer(:,:,iphase)
1585 else if (ll .eq. 4) then
1586 lptr_dum(:,:) = numptr_aer(:,:,iphase)
1587 end if
1589 if (ll .eq. 1) write(lunxx,9800)
1590 if (ll .eq. 1) write(lunxx,9800)
1591 write(lunxx,9800)
1592 write(lunxx,9800)
1593 write(lunxx,9801) 'cloudchem_apply_move_sections - ', &
1594 name(lptr_dum(1,1)), ' - icase, ll'
1595 write(lunxx,9810) icase, ll
1596 write(lunxx,9800) ' tp sz rbox_sv1, rbox, rsub'
1598 do itype = 1, ntype_aer
1599 do isize = 1, nsize_aer(itype)
1600 l = lptr_dum(isize,itype)
1601 dumb = rbox(l) - rbox_sv1(l)
1602 dumc = dumb/max( abs(duma), 1.0e-35 )
1603 if (duma .lt. 0.0) dumc = -dumc
1604 write(lunxx,9840) itype, isize, rbox_sv1(l), rbox(l), &
1605 rsub(l,kpeg,mpeg)
1606 end do
1607 end do
1608 end do
1610 if (icase .ge. 5) then
1611 write(*,*) &
1612 '*** stop in cloudchem_apply_move_sections diags, icase =', &
1613 icase
1614 ! stop
1615 end if
1616 end if
1617 9800 format( a )
1618 9801 format( 3a )
1619 9810 format( 7i10 )
1620 9820 format( 7(1pe10.2) )
1621 9840 format( 2i3, 5(1pe14.6) )
1622 #endif
1625 return
1626 end subroutine cloudchem_apply_move_sections
1630 !-----------------------------------------------------------------------
1631 subroutine mosaic_cloudchem_dumpaa( &
1632 id, ktau, ktauc, dtstepc, config_flags, &
1633 p_phy, t_phy, rho_phy, alt, &
1634 cldfra, ph_no2, &
1635 moist, chem, &
1636 gas_aqfrac, numgas_aqfrac, &
1637 ids,ide, jds,jde, kds,kde, &
1638 ims,ime, jms,jme, kms,kme, &
1639 its,ite, jts,jte, kts,kte, &
1640 qcldwtr_cutoff, &
1641 itcur, jtcur, ktcur )
1643 use module_state_description, only: &
1644 num_moist, num_chem, p_qc
1646 use module_configure, only: grid_config_rec_type
1648 use module_data_mosaic_asect
1650 use module_data_mosaic_other, only: k_pegbegin, name
1652 use module_mosaic_driver, only: mapaer_tofrom_host
1655 implicit none
1657 ! subr arguments
1658 integer, intent(in) :: &
1659 id, ktau, ktauc, &
1660 numgas_aqfrac, &
1661 ids, ide, jds, jde, kds, kde, &
1662 ims, ime, jms, jme, kms, kme, &
1663 its, ite, jts, jte, kts, kte, &
1664 itcur, jtcur, ktcur
1665 ! id - domain index
1666 ! ktau - time step number
1667 ! ktauc - gas and aerosol chemistry time step number
1668 ! numgas_aqfrac - last dimension of gas_aqfrac
1670 ! [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for 'domain'
1671 ! [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for 'memory'
1672 ! Most arrays that are arguments to chem_driver
1673 ! are dimensioned with these spatial indices.
1674 ! [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for 'tile'
1675 ! chem_driver and routines under it do calculations
1676 ! over these spatial indices.
1678 type(grid_config_rec_type), intent(in) :: config_flags
1679 ! config_flags - configuration and control parameters
1681 real, intent(in) :: &
1682 dtstepc, qcldwtr_cutoff
1683 ! dtstepc - time step for gas and aerosol chemistry(s)
1685 real, intent(in), &
1686 dimension( ims:ime, kms:kme, jms:jme ) :: &
1687 p_phy, t_phy, rho_phy, alt, cldfra, ph_no2
1688 ! p_phy - air pressure (Pa)
1689 ! t_phy - temperature (K)
1690 ! rho_phy - moist air density (kg/m^3)
1691 ! alt - dry air specific volume (m^3/kg)
1692 ! cldfra - cloud fractional area (0-1)
1693 ! ph_no2 - no2 photolysis rate (1/min)
1695 real, intent(in), &
1696 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
1697 moist
1698 ! moist - mixing ratios of moisture species (water vapor,
1699 ! cloud water, ...) (kg/kg for mass species, #/kg for number species)
1701 real, intent(inout), &
1702 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
1703 chem
1704 ! chem - mixing ratios of trace gas and aerosol species (ppm for gases,
1705 ! ug/kg for aerosol mass species, #/kg for aerosol number species)
1707 real, intent(inout), &
1708 dimension( ims:ime, kms:kme, jms:jme, numgas_aqfrac ) :: &
1709 gas_aqfrac
1710 ! gas_aqfrac - fraction (0-1) of gas that is dissolved in cloud water
1713 ! local variables
1714 integer :: it, jt, kt, l, ll, n
1715 integer :: isize, itype
1717 real :: dumai, dumcw
1718 real :: qcldwtr
1721 it = itcur
1722 jt = jtcur
1723 kt = ktcur
1725 write(*,*)
1726 write(*,*)
1727 write(*,*)
1728 write(*,9100)
1729 write(*,9102) ktau, it, jt, kt
1730 9100 format( 7('----------') )
1731 9102 format( &
1732 'mosaic_cloudchem_dumpaa - ktau, i, j, k =', 4i5 )
1734 itype = 1
1735 do 2900 isize = 1, nsize_aer(itype)
1737 write(*,9110) isize
1738 9110 format( / 'isize =', i3 / &
1739 ' k cldwtr mass-ai numb-ai mass-cw numb-cw' )
1741 do 2800 kt = kte, kts, -1
1743 dumai = 0.0
1744 dumcw = 0.0
1745 do ll = 1, ncomp_aer(itype)
1746 l = massptr_aer(ll,isize,itype,1)
1747 dumai = dumai + chem(it,kt,jt,l)
1748 l = massptr_aer(ll,isize,itype,2)
1749 dumcw = dumcw + chem(it,kt,jt,l)
1750 end do
1751 write(*,9120) kt, &
1752 moist(it,kt,jt,p_qc), &
1753 dumai, chem(it,kt,jt,numptr_aer(isize,itype,1)), &
1754 dumcw, chem(it,kt,jt,numptr_aer(isize,itype,2))
1755 9120 format( i3, 1p, e10.2, 2(3x, 2e10.2) )
1757 2800 continue
1758 2900 continue
1760 write(*,*)
1761 write(*,9100)
1762 write(*,*)
1764 ! map from wrf-chem 3d arrays to pegasus clm & sub arrays
1765 kt = ktcur
1766 if ((ktau .eq. 30) .and. (it .eq. 23) .and. &
1767 (jt .eq. 1) .and. (kt .eq. 11)) then
1768 qcldwtr = moist(it,kt,jt,p_qc)
1769 write(*,*)
1770 write(*,*)
1771 write(*,9102) ktau, it, jt, kt
1772 write(*,*)
1773 write( *, '(3(1pe10.2,3x,a))' ) &
1774 (chem(it,kt,jt,l), name(l)(1:10), l=1,num_chem)
1775 write(*,*)
1776 write( *, '(3(1pe10.2,3x,a))' ) &
1777 p_phy(it,kt,jt), 'p_phy ', &
1778 t_phy(it,kt,jt), 't_phy ', &
1779 rho_phy(it,kt,jt), 'rho_phy ', &
1780 alt(it,kt,jt), 'alt ', &
1781 qcldwtr, 'qcldwtr ', &
1782 qcldwtr_cutoff, 'qcldwtrcut'
1783 write(*,*)
1784 write(*,9100)
1785 write(*,*)
1786 end if
1789 return
1790 end subroutine mosaic_cloudchem_dumpaa
1794 !-----------------------------------------------------------------------
1795 end module module_mosaic_cloudchem