!**********************************************************************************
! This computer software was prepared by Battelle Memorial Institute, hereinafter
! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! MOSAIC module: see module_mosaic_driver.F for references and terms of use
!**********************************************************************************
module module_mosaic_cloudchem
integer, parameter :: l_so4_aqyy = 1
integer, parameter :: l_no3_aqyy = 2
integer, parameter :: l_cl_aqyy = 3
integer, parameter :: l_nh4_aqyy = 4
integer, parameter :: l_na_aqyy = 5
integer, parameter :: l_oin_aqyy = 6
integer, parameter :: l_bc_aqyy = 7
integer, parameter :: l_oc_aqyy = 8
integer, parameter :: nyyy = 8
contains
!-----------------------------------------------------------------------
subroutine mosaic_cloudchem_driver( &
id, ktau, ktauc, dtstepc, config_flags, &
p_phy, t_phy, rho_phy, alt, &
cldfra, ph_no2, &
moist, chem, &
gas_aqfrac, numgas_aqfrac, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
use module_state_description, only: &
num_moist, num_chem, p_qc
use module_configure, only: grid_config_rec_type
use module_data_mosaic_asect, only: cw_phase, nphase_aer
use module_data_mosaic_other, only: k_pegbegin, name
use module_mosaic_driver, only: mapaer_tofrom_host
implicit none
! subr arguments
integer, intent(in) :: &
id, ktau, ktauc, &
numgas_aqfrac, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte
! id - domain index
! ktau - time step number
! ktauc - gas and aerosol chemistry time step number
! numgas_aqfrac - last dimension of gas_aqfrac
! [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for 'domain'
! [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for 'memory'
! Most arrays that are arguments to chem_driver
! are dimensioned with these spatial indices.
! [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for 'tile'
! chem_driver and routines under it do calculations
! over these spatial indices.
type(grid_config_rec_type), intent(in) :: config_flags
! config_flags - configuration and control parameters
real, intent(in) :: &
dtstepc
! dtstepc - time step for gas and aerosol chemistry(s)
real, intent(in), &
dimension( ims:ime, kms:kme, jms:jme ) :: &
p_phy, t_phy, rho_phy, alt, cldfra, ph_no2
! p_phy - air pressure (Pa)
! t_phy - temperature (K)
! rho_phy - moist air density (kg/m^3)
! alt - dry air specific volume (m^3/kg)
! cldfra - cloud fractional area (0-1)
! ph_no2 - no2 photolysis rate (1/min)
real, intent(in), &
dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
moist
! moist - mixing ratios of moisture species (water vapor,
! cloud water, ...) (kg/kg for mass species, #/kg for number species)
real, intent(inout), &
dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
chem
! chem - mixing ratios of trace gas and aerosol species (ppm for gases,
! ug/kg for aerosol mass species, #/kg for aerosol number species)
real, intent(inout), &
dimension( ims:ime, kms:kme, jms:jme, numgas_aqfrac ) :: &
gas_aqfrac
! gas_aqfrac - fraction (0-1) of gas that is dissolved in cloud water
! local variables
integer :: it, jt, kt, kpeg, k_pegshift, l, mpeg
integer :: icase
integer :: igaschem_onoff, iphotol_onoff, iradical_onoff
real :: gas_aqfrac_box(numgas_aqfrac)
real :: ph_aq_box
real, parameter :: qcldwtr_cutoff = 1.0e-6
real :: qcldwtr
! check that cw_phase is active
if ((cw_phase .le. 0) .or. (cw_phase .gt. nphase_aer)) then
print *, '*** mosaic_cloudchem_driver - cw_phase not active'
return
end if
print 93010, 'entering mosaic_cloudchem_driver - ktau =', ktau
icase = 0
! iphotol_onoff = 1 if photolysis rate calcs are on; 0 if off
iphotol_onoff = 0
if (config_flags%phot_opt .gt. 0) iphotol_onoff = 1
! igaschem_onoff = 1 if gas-phase chemistry is on; 0 if off
igaschem_onoff = 0
if (config_flags%gaschem_onoff .gt. 0) igaschem_onoff = 1
! iradical_onoff turns aqueous radical chemistry on/off
! set iradical_onoff=0 if either photolysis or gas-phase chem are off
if ((igaschem_onoff .le. 0) .or. (iphotol_onoff .le. 0)) then
iradical_onoff = 0
else
iradical_onoff = 1
end if
! following line turns aqueous radical chem off unconditionally
iradical_onoff = 0
do 3920 jt = jts, jte
do 3910 it = its, ite
do 3800 kt = kts, kte
qcldwtr = moist(it,kt,jt,p_qc)
if (qcldwtr .le. qcldwtr_cutoff) goto 3800
k_pegshift = k_pegbegin - kts
kpeg = kt + k_pegshift
mpeg = 1
icase = icase + 1
! detailed dump for debugging
if (ktau .eq. -13579) then
! if ((ktau .eq. 30) .and. (it .eq. 23) .and. &
! (jt .eq. 1) .and. (kt .eq. 11)) then
call mosaic_cloudchem_dumpaa( &
id, ktau, ktauc, dtstepc, config_flags, &
p_phy, t_phy, rho_phy, alt, &
cldfra, ph_no2, &
moist, chem, &
gas_aqfrac, numgas_aqfrac, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
qcldwtr_cutoff, &
it, jt, kt )
end if
! map from wrf-chem 3d arrays to pegasus clm & sub arrays
call mapaer_tofrom_host( 0, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
it, jt, kt, kt, &
num_moist, num_chem, moist, chem, &
t_phy, p_phy, rho_phy )
! make call '1box' cloudchem routine
! print 93010, 'calling mosaic_cloudchem_1 at ijk =', it, jt, kt
call mosaic_cloudchem_1box( &
id, ktau, ktauc, dtstepc, &
iphotol_onoff, iradical_onoff, &
ph_no2(it,kt,jt), &
ph_aq_box, gas_aqfrac_box, &
numgas_aqfrac, it, jt, kt, kpeg, mpeg, icase )
! map back to wrf-chem 3d arrays
call mapaer_tofrom_host( 1, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
it, jt, kt, kt, &
num_moist, num_chem, moist, chem, &
t_phy, p_phy, rho_phy )
gas_aqfrac(it,kt,jt,:) = gas_aqfrac_box(:)
3800 continue
3910 continue
3920 continue
print 93010, 'leaving mosaic_cloudchem_driver - ktau =', ktau, icase
93010 format( a, 8(1x,i6) )
return
end subroutine mosaic_cloudchem_driver
!-----------------------------------------------------------------------
subroutine mosaic_cloudchem_1box( &
id, ktau, ktauc, dtstepc, &
iphotol_onoff, iradical_onoff, &
photol_no2_box, &
ph_aq_box, gas_aqfrac_box, &
numgas_aqfrac, it, jt, kt, kpeg, mpeg, icase )
use module_state_description, only: &
num_moist, num_chem
use module_data_mosaic_asect, only: &
msectional, &
maxd_asize, maxd_atype, &
cw_phase, nsize_aer, ntype_aer, &
lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
lptr_oin_aer, lptr_bc_aer, lptr_oc_aer
use module_data_mosaic_other, only: &
l2maxd, ltot2, rsub
use module_data_cmu_bulkaqchem, only: &
meqn1max
implicit none
! subr arguments
integer, intent(in) :: &
id, ktau, ktauc, &
numgas_aqfrac, it, jt, kt, kpeg, mpeg, &
icase, iphotol_onoff, iradical_onoff
real, intent(in) :: &
dtstepc, photol_no2_box
real, intent(inout) :: ph_aq_box
real, intent(inout), dimension( numgas_aqfrac ) :: gas_aqfrac_box
! local variables
integer :: iphase
integer :: icase_in, idecomp_hmsa_hso5, &
iradical_in, istat_aqop
integer :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
real :: co2_mixrat_in
real :: ph_cmuaq_cur
real :: photol_no2_in
real :: xprescribe_ph
real :: yaq_beg(meqn1max), yaq_end(meqn1max)
real :: rbox(l2maxd), rbox_sv1(l2maxd)
real :: rbulk_cwaer(nyyy,2)
real, dimension( maxd_asize, maxd_atype ) :: fr_partit_cw
!
! set the lptr_yyy_cwaer
!
iphase = cw_phase
lptr_yyy_cwaer(:,:,l_so4_aqyy) = lptr_so4_aer(:,:,iphase)
lptr_yyy_cwaer(:,:,l_no3_aqyy) = lptr_no3_aer(:,:,iphase)
lptr_yyy_cwaer(:,:,l_cl_aqyy ) = lptr_cl_aer( :,:,iphase)
lptr_yyy_cwaer(:,:,l_nh4_aqyy) = lptr_nh4_aer(:,:,iphase)
lptr_yyy_cwaer(:,:,l_na_aqyy ) = lptr_na_aer( :,:,iphase)
lptr_yyy_cwaer(:,:,l_oin_aqyy) = lptr_oin_aer(:,:,iphase)
lptr_yyy_cwaer(:,:,l_bc_aqyy ) = lptr_bc_aer( :,:,iphase)
lptr_yyy_cwaer(:,:,l_oc_aqyy ) = lptr_oc_aer( :,:,iphase)
!
! xfer from rsub to rbox
!
rbox(1:ltot2) = max( 0.0, rsub(1:ltot2,kpeg,mpeg) )
rbox_sv1(1:ltot2) = rbox(1:ltot2)
!
!
! do bulk cloud-water chemistry
!
!
icase_in = icase
iradical_in = 1
idecomp_hmsa_hso5 = 1
co2_mixrat_in = 350.0
photol_no2_in = photol_no2_box
! when xprescribe_ph >= 0, ph is forced to its value
xprescribe_ph = -1.0e31
! turn off aqueous phase photolytic and radical chemistry
! if either of the iphotol_onoff and iradical_onoff flags are 0
if ((iphotol_onoff .le. 0) .or. (iradical_onoff .le. 0)) then
photol_no2_in = 0.0
iradical_in = 0
end if
#if defined ( ccboxtest_box_testing_active)
! following is for off-line box testing only
call ccboxtest_extra_args_aa( 'get', &
co2_mixrat_in, iradical_in, &
idecomp_hmsa_hso5, icase_in, &
xprescribe_ph )
#endif
gas_aqfrac_box(:) = 0.0
! make call to interface_to_aqoperator1
call interface_to_aqoperator1( &
istat_aqop, &
dtstepc, &
rbox, gas_aqfrac_box, &
rbulk_cwaer, lptr_yyy_cwaer, &
co2_mixrat_in, photol_no2_in, xprescribe_ph, &
iradical_in, idecomp_hmsa_hso5, &
yaq_beg, yaq_end, ph_cmuaq_cur, &
numgas_aqfrac, id, it, jt, kt, kpeg, mpeg, ktau, icase_in )
ph_aq_box = ph_cmuaq_cur
#if defined ( ccboxtest_box_testing_active)
! following is for off-line box testing only
call ccboxtest_extra_args_bb( 'put', &
yaq_beg, yaq_end, ph_cmuaq_cur )
#endif
!
!
! calculate fraction of cloud-water associated with each activated aerosol bin
!
!
call partition_cldwtr( &
rbox, fr_partit_cw, &
it, jt, kt, kpeg, mpeg, icase_in )
!
!
! distribute changes in bulk cloud-water composition among size bins
!
!
call distribute_bulk_changes( &
rbox, rbox_sv1, fr_partit_cw, &
rbulk_cwaer, lptr_yyy_cwaer, &
it, jt, kt, kpeg, mpeg, icase_in )
!
! xfer back to rsub
!
rsub(1:ltot2,kpeg,mpeg) = max( 0.0, rbox(1:ltot2) )
!
! do move-sections
!
if (msectional .lt. 1000000000) then
call cloudchem_apply_move_sections( &
rbox, rbox_sv1, &
it, jt, kt, kpeg, mpeg, icase_in )
end if
return
end subroutine mosaic_cloudchem_1box
!-----------------------------------------------------------------------
subroutine interface_to_aqoperator1( &
istat_aqop, &
dtstepc, &
rbox, gas_aqfrac_box, &
rbulk_cwaer, lptr_yyy_cwaer, &
co2_mixrat_in, photol_no2_in, xprescribe_ph, &
iradical_in, idecomp_hmsa_hso5, &
yaq_beg, yaq_end, ph_cmuaq_cur, &
numgas_aqfrac, id, it, jt, kt, kpeg, mpeg, ktau, icase )
use module_state_description, only: &
num_chem, param_first_scalar, p_qc, &
p_nh3, p_hno3, p_hcl, p_sulf, p_h2so4, p_hcho, &
p_ora1, p_so2, p_h2o2, p_o3, p_ho, &
p_ho2, p_no3, p_no, p_no2, p_hono, &
! p_pan, p_ch3o2, p_ch3oh, p_op1
p_pan, p_ch3o2, p_ch3oh, p_op1, &
p_hcooh, p_ch3ooh
use module_data_cmu_bulkaqchem, only: &
meqn1max, naers, ngas, &
na4, naa, nac, nae, nah, nahmsa, nahso5, &
nan, nao, nar, nas, naw, &
ng4, nga, ngc, ngch3co3h, ngch3o2, ngch3o2h, ngch3oh, &
ngh2o2, nghcho, nghcooh, nghno2, ngho2, &
ngn, ngno, ngno2, ngno3, ngo3, ngoh, ngpan, ngso2
use module_cmu_bulkaqchem, only: aqoperator1
use module_data_mosaic_asect, only: &
maxd_asize, maxd_atype, &
cw_phase, nsize_aer, ntype_aer, &
lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
lptr_oin_aer, lptr_bc_aer, lptr_oc_aer, &
mw_cl_aer, mw_na_aer, mw_nh4_aer, mw_no3_aer, mw_so4_aer
use module_data_mosaic_other, only: &
aboxtest_units_convert, cairclm, &
ktemp, l2maxd, ptotclm, rcldwtr_sub
implicit none
! subr arguments
integer, intent(in) :: &
iradical_in, idecomp_hmsa_hso5, &
numgas_aqfrac, id, it, jt, kt, kpeg, mpeg, ktau, icase
integer, intent(inout) :: &
istat_aqop
integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
real, intent(in) :: &
dtstepc, co2_mixrat_in, &
photol_no2_in, xprescribe_ph
real, intent(inout) :: ph_cmuaq_cur
real, intent(inout), dimension( 1:l2maxd ) :: rbox
real, intent(inout), dimension( nyyy, 2 ) :: rbulk_cwaer
real, intent(inout), dimension( 1:numgas_aqfrac ) :: gas_aqfrac_box
real, intent(inout), dimension( meqn1max ) :: yaq_beg, yaq_end
! local variables
integer :: i, iphase, isize, itype
integer :: iaq, istat_fatal, istat_warn
integer :: l, lunxx, lyyy
integer :: p1st
real, parameter :: eps=0.622 ! (mw h2o)/(mw air)
real :: cair_moleperm3
real :: dum, dumb
real :: factgas, factlwc, factpatm, factphoto
real :: factaerbc, factaercl, factaerna, factaernh4, &
factaerno3, factaeroc, factaeroin, factaerso4
real :: lwc
real :: p_atm, photo_in
real :: rh
real :: temp, tstep_beg_sec, tstep_end_sec
real :: totsulf_beg, totsulf_end
real :: gas(ngas), aerosol(naers)
real :: gas_aqfrac_cmu(ngas)
double precision tstep_beg_sec_dp, tstep_end_sec_dp, &
temp_dp, p_atm_dp, lwc_dp, rh_dp, &
co2_mixrat_in_dp, photo_in_dp, ph_cmuaq_cur_dp, &
xprescribe_ph_dp
double precision gas_dp(ngas), gas_aqfrac_cmu_dp(ngas), &
aerosol_dp(naers), yaq_beg_dp(meqn1max), yaq_end_dp(meqn1max)
p1st = param_first_scalar
!
! units conversion factors
! 'cmuaq-bulk' value = pegasus value X factor
!
! [pres in atmospheres] = [pres in dynes/cm2] * factpatm
factpatm = 1.0/1.01325e6
! [cldwtr in g-h2o/m3-air] = [cldwtr in mole-h2o/mole-air] * factlwc
factlwc = 28.966*eps*1.0e6*cairclm(kpeg)
! [aq photolysis rate scaling factor in --] = [jno2 in 1/min] * factphoto
factphoto = 1.6
! [gas in ppm] = [gas in mole/mole-air] * factgas
factgas = 1.0e6
! [aerosol in ug/m3-air] = [aerosol in mole/mole-air] * factaer
dum = cairclm(kpeg)*1.0e12
factaerso4 = dum*mw_so4_aer
factaerno3 = dum*mw_no3_aer
factaercl = dum*mw_cl_aer
factaernh4 = dum*mw_nh4_aer
factaerna = dum*mw_na_aer
factaeroin = dum
factaeroc = dum
factaerbc = dum
! rce 2005-jul-11 - use same molecular weights here as in cmu code
! factaerso4 = dum*96.0
! factaerno3 = dum*62.0
! factaercl = dum*35.5
! factaernh4 = dum*18.0
! factaerna = dum*23.0
! If aboxtest_units_convert=10, turn off units conversions both here
! and in module_mosaic. This is for testing, to allow exact agreements.
if (aboxtest_units_convert .eq. 10) then
factpatm = 1.0
factlwc = 1.0
factphoto = 1.0
factgas = 1.0
factaerso4 = 1.0
factaerno3 = 1.0
factaercl = 1.0
factaernh4 = 1.0
factaerna = 1.0
factaeroin = 1.0
factaeroc = 1.0
factaerbc = 1.0
end if
!
! map from rbox to gas,aerosol
!
temp = rbox(ktemp)
lwc = rcldwtr_sub(kpeg,mpeg) * factlwc
p_atm = ptotclm(kpeg) * factpatm
! rce 2005-jul-11 - set p_atm so that cmu code's cair will match cairclm
p_atm = cairclm(kpeg)*1.0e3*0.082058e0*temp
photo_in = photol_no2_in * factphoto
rh = 1.0
iaq = 1
tstep_beg_sec = 0.0
tstep_end_sec = dtstepc
! map gases and convert to ppm
gas(:) = 0.0
gas(nga ) = rbox(p_nh3 ) * factgas
gas(ngn ) = rbox(p_hno3 ) * factgas
if(p_hcl > param_first_scalar ) gas(ngc ) = rbox(p_hcl ) * factgas
if(p_sulf > param_first_scalar ) gas(ng4 ) = rbox(p_sulf ) * factgas
if(p_h2so4 > param_first_scalar ) gas(ng4 ) = rbox(p_h2so4 ) * factgas
gas(nghcho ) = rbox(p_hcho ) * factgas
if(p_ora1 > param_first_scalar ) gas(nghcooh ) = rbox(p_ora1 ) * factgas
if(p_hcooh > param_first_scalar ) gas(nghcooh ) = rbox(p_hcooh ) * factgas
gas(ngso2 ) = rbox(p_so2 ) * factgas
gas(ngh2o2 ) = rbox(p_h2o2 ) * factgas
gas(ngo3 ) = rbox(p_o3 ) * factgas
gas(ngoh ) = rbox(p_ho ) * factgas
gas(ngho2 ) = rbox(p_ho2 ) * factgas
gas(ngno3 ) = rbox(p_no3 ) * factgas
gas(ngno ) = rbox(p_no ) * factgas
gas(ngno2 ) = rbox(p_no2 ) * factgas
gas(nghno2 ) = rbox(p_hono ) * factgas
gas(ngpan ) = rbox(p_pan ) * factgas
gas(ngch3o2 ) = rbox(p_ch3o2 ) * factgas
gas(ngch3oh ) = rbox(p_ch3oh ) * factgas
if(p_op1 > param_first_scalar ) gas(ngch3o2h) = rbox(p_op1 ) * factgas
if(p_ch3ooh > param_first_scalar ) gas(ngch3o2h) = rbox(p_ch3ooh) * factgas
! compute bulk activated-aerosol mixing ratios
aerosol(:) = 0.0
rbulk_cwaer(:,:) = 0.0
iphase = cw_phase
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
do lyyy = 1, nyyy
l = lptr_yyy_cwaer(isize,itype,lyyy)
if (l .ge. p1st) rbulk_cwaer(lyyy,1) = rbulk_cwaer(lyyy,1) + rbox(l)
end do
end do
end do
! map them to 'aerosol' array and convert to ug/m3
aerosol(na4) = rbulk_cwaer(l_so4_aqyy,1) * factaerso4
aerosol(nan) = rbulk_cwaer(l_no3_aqyy,1) * factaerno3
aerosol(nac) = rbulk_cwaer(l_cl_aqyy, 1) * factaercl
aerosol(naa) = rbulk_cwaer(l_nh4_aqyy,1) * factaernh4
aerosol(nas) = rbulk_cwaer(l_na_aqyy, 1) * factaerna
aerosol(nar) = rbulk_cwaer(l_oin_aqyy,1) * factaeroin
aerosol(nae) = rbulk_cwaer(l_bc_aqyy, 1) * factaerbc
aerosol(nao) = rbulk_cwaer(l_oc_aqyy, 1) * factaeroc
!
! make call to aqoperator1
!
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
lunxx = 87
lunxx = -1
if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
if (lunxx .gt. 0) then
write(lunxx,*)
write(lunxx,*)
write(lunxx,*) 'interface_to_aqoperator1 - icase, irad, idecomp'
write(lunxx,9870) icase, iradical_in, idecomp_hmsa_hso5
write(lunxx,*) 'it, jt, kt, kpeg, mpeg, ktau'
write(lunxx,9870) it, jt, kt, kpeg, mpeg, ktau
write(lunxx,*) 'temp, p_atm, lwc, photo, co2, xprescribe_ph'
write(lunxx,9875) temp, p_atm, lwc, photo_in, co2_mixrat_in, xprescribe_ph
write(lunxx,*) 'ptot, cair, rcldwtr_clm, dt_sec'
write(lunxx,9875) ptotclm(kpeg), cairclm(kpeg), &
rcldwtr_sub(kpeg,mpeg), (tstep_end_sec-tstep_beg_sec)
write(lunxx,*) 'gas (1=nh3, 2=hno3, 3=hcl, 4=h2so4, 11=so2, 12=h2o2, 18=o3)'
write(lunxx,9875) gas
write(lunxx,*) 'rbox(nh3, hno3, hcl, h2so4, so2, h2o2, o3)'
write(lunxx,9875) rbox(p_nh3), rbox(p_hno3), rbox(p_hcl), &
rbox(p_sulf),rbox(p_h2so4), rbox(p_so2), rbox(p_h2o2), rbox(p_o3)
write(lunxx,*) 'aerosol (1=na, 3=nh4, 4=no3, 5=cl, 6=so4, 8=ec, 9=oc, 10=crus)'
write(lunxx,9875) aerosol
write(lunxx,*) 'rbulk_cwaer (1=so4, 2=no3, 3-cl, 4=nh4, 5=na, 6=oin, 7=bc, 8=oc)'
write(lunxx,9875) rbulk_cwaer(:,1)
! if (icase .ge. 3) then
! write(*,*) &
! '*** stopping in interface_to_aqop1 at icase =', icase
! stop
! end if
end if
9870 format( 8i5 )
9875 format( 5(1pe14.6) )
#endif
#if 0
! Print outs for debugging of aqoperator1... wig, 26-Oct-2005
!!$ if( (id == 1 .and. ktau >= 207 ) .or. &
!!$ (id == 2 .and. ktau >= 610 ) .or. &
!!$ (id == 3 .and. ktau >= 1830 ) ) then
write(6,'(a)') '---Begin input for aqoperator1---'
write(6,'(a,4i)') 'id, it, jt, kt =', id, it, jt, kt
write(6,'(a,1p,2e20.12)') 'tstep_beg_sec, tstep_end_sec = ',tstep_beg_sec, tstep_end_sec
do l=1,ngas
write(6,'("gas(",i2,") = ",1p,1e20.12)') l, gas(l)
end do
do l=1,naers
write(6,'("aerosol(",i2,") = ",1p,1e20.12)') l, aerosol(l)
end do
write(6,'(a,1p,4e20.12)') "temp, p_atm, lwc, rh = ", temp, p_atm, lwc, rh
write(6,'(a,1p,3e20.12)') "co2_mixrat_in, photo_in, xprescribe_ph = ", &
co2_mixrat_in, photo_in, xprescribe_ph
write(6,'(a,3i)') " iradical_in, idecomp_hmsa_hso5, iaq = ", iradical_in, idecomp_hmsa_hso5, iaq
write(6,'(a)') "---End input for aqoperator1---"
!!$ end if
#endif
! convert arguments to double prec
tstep_beg_sec_dp = 0.0d0
if (tstep_beg_sec .ne. 0.0) tstep_beg_sec_dp = tstep_beg_sec
tstep_end_sec_dp = 0.0d0
if (tstep_end_sec .ne. 0.0) tstep_end_sec_dp = tstep_end_sec
temp_dp = 0.0d0
if (temp .ne. 0.0) temp_dp = temp
p_atm_dp = 0.0d0
if (p_atm .ne. 0.0) p_atm_dp = p_atm
lwc_dp = 0.0d0
if (lwc .ne. 0.0) lwc_dp = lwc
rh_dp = 0.0d0
if (rh .ne. 0.0) rh_dp = rh
co2_mixrat_in_dp = 0.0d0
if (co2_mixrat_in .ne. 0.0) co2_mixrat_in_dp = co2_mixrat_in
photo_in_dp = 0.0d0
if (photo_in .ne. 0.0) photo_in_dp = photo_in
xprescribe_ph_dp = 0.0d0
if (xprescribe_ph .ne. 0.0) xprescribe_ph_dp = xprescribe_ph
ph_cmuaq_cur_dp = 0.0d0
if (ph_cmuaq_cur .ne. 0.0) ph_cmuaq_cur_dp = ph_cmuaq_cur
do i = 1, ngas
gas_dp(i) = 0.0d0
if (gas(i) .ne. 0.0) gas_dp(i) = gas(i)
end do
do i = 1, naers
aerosol_dp(i) = 0.0d0
if (aerosol(i) .ne. 0.0) aerosol_dp(i) = aerosol(i)
end do
do i = 1, ngas
gas_aqfrac_cmu_dp(i) = 0.0d0
if (gas_aqfrac_cmu(i) .ne. 0.0) gas_aqfrac_cmu_dp(i) = gas_aqfrac_cmu(i)
end do
do i = 1, meqn1max
yaq_beg_dp(i) = 0.0d0
if (yaq_beg(i) .ne. 0.0) yaq_beg_dp(i) = yaq_beg(i)
end do
do i = 1, meqn1max
yaq_end_dp(i) = 0.0d0
if (yaq_end(i) .ne. 0.0) yaq_end_dp(i) = yaq_end(i)
end do
! total sulfur species conc as sulfate (ug/m3)
cair_moleperm3 = 1.0e3*p_atm_dp/(0.082058e0*temp_dp)
totsulf_beg = ( aerosol_dp(na4)/96. &
+ aerosol_dp(nahso5)/113. + aerosol_dp(nahmsa)/111. &
+ (gas_dp(ngso2) + gas_dp(ng4))*cair_moleperm3 )*96.0
! call aqoperator1( &
! istat_fatal, istat_warn, &
! tstep_beg_sec, tstep_end_sec, &
! gas, aerosol, gas_aqfrac_cmu, &
! temp, p_atm, lwc, rh, &
! co2_mixrat_in, photo_in, xprescribe_ph, &
! iradical_in, idecomp_hmsa_hso5, iaq, &
! yaq_beg, yaq_end, ph_cmuaq_cur )
call aqoperator1( &
istat_fatal, istat_warn, &
tstep_beg_sec_dp, tstep_end_sec_dp, &
gas_dp, aerosol_dp, gas_aqfrac_cmu_dp, &
temp_dp, p_atm_dp, lwc_dp, rh_dp, &
co2_mixrat_in_dp, photo_in_dp, xprescribe_ph_dp, &
iradical_in, idecomp_hmsa_hso5, iaq, &
yaq_beg_dp, yaq_end_dp, ph_cmuaq_cur_dp )
totsulf_end = ( aerosol_dp(na4)/96. &
+ aerosol_dp(nahso5)/113. + aerosol_dp(nahmsa)/111. &
+ (gas_dp(ngso2) + gas_dp(ng4))*cair_moleperm3 )*96.0
! convert arguments back to single prec
tstep_beg_sec = tstep_beg_sec_dp
tstep_end_sec = tstep_end_sec_dp
temp = temp_dp
p_atm = p_atm_dp
lwc = lwc_dp
rh = rh_dp
! co2_mixrat_in = co2_mixrat_in_dp ! this has intent(in)
! photo_in = photo_in_dp ! this has intent(in)
! xprescribe_ph = xprescribe_ph_dp ! this has intent(in)
ph_cmuaq_cur = ph_cmuaq_cur_dp
do i = 1, ngas
gas(i) = gas_dp(i)
end do
do i = 1, naers
aerosol(i) = aerosol_dp(i)
end do
do i = 1, ngas
gas_aqfrac_cmu(i) = gas_aqfrac_cmu_dp(i)
end do
do i = 1, meqn1max
yaq_beg(i) = yaq_beg_dp(i)
end do
do i = 1, meqn1max
yaq_end(i) = yaq_end_dp(i)
end do
!
! warning message when status flags are non-zero
!
istat_aqop = 0
if (istat_fatal .ne. 0) then
write(6,*) &
'*** mosaic_cloudchem_driver, subr interface_to_aqoperator1'
write(6,'(a,4i5,2i10)') &
' id,it,jt,kt, istat_fatal, warn =', &
id, it, jt, kt, istat_fatal, istat_warn
istat_aqop = -10
end if
!
! warning message when sulfur mass balance error exceeds the greater
! of (1.0e-3 ug/m3) OR (1.0e-3 X total sulfur mixing ratio)
!
dum = totsulf_end - totsulf_beg
dumb = max( totsulf_beg, totsulf_end )
if (abs(dum) .gt. max(1.0e-3,1.0e-3*dumb)) then
write(6,*) &
'*** mosaic_cloudchem_driver, sulfur balance warning'
write(6,'(a,4i5,1p,3e12.4)') &
' id,it,jt,kt, total_sulfur_beg, _end, _error =', &
id, it, jt, kt, totsulf_beg, totsulf_end, dum
end if
!
! map from gas,aerosol to rbox
!
rbox(p_nh3 ) = gas(nga ) / factgas
rbox(p_hno3 ) = gas(ngn ) / factgas
if(p_hcl .gt. param_first_scalar) rbox(p_hcl ) = gas(ngc ) / factgas
if(p_sulf .gt. param_first_scalar) rbox(p_sulf ) = gas(ng4) / factgas
if(p_h2so4 .gt. param_first_scalar) rbox(p_h2so4 ) = gas(ng4) / factgas
rbox(p_hcho ) = gas(nghcho ) / factgas
if(p_ora1 > param_first_scalar ) rbox(p_ora1 ) = gas(nghcooh ) / factgas
if(p_hcooh > param_first_scalar ) rbox(p_hcooh) = gas(nghcooh ) / factgas
rbox(p_so2 ) = gas(ngso2 ) / factgas
rbox(p_h2o2 ) = gas(ngh2o2 ) / factgas
rbox(p_o3 ) = gas(ngo3 ) / factgas
rbox(p_ho ) = gas(ngoh ) / factgas
rbox(p_ho2 ) = gas(ngho2 ) / factgas
rbox(p_no3 ) = gas(ngno3 ) / factgas
rbox(p_no ) = gas(ngno ) / factgas
rbox(p_no2 ) = gas(ngno2 ) / factgas
rbox(p_hono ) = gas(nghno2 ) / factgas
rbox(p_pan ) = gas(ngpan ) / factgas
rbox(p_ch3o2 ) = gas(ngch3o2 ) / factgas
rbox(p_ch3oh ) = gas(ngch3oh ) / factgas
if(p_op1 > param_first_scalar ) rbox(p_op1 ) = gas(ngch3o2h) / factgas
if(p_ch3ooh > param_first_scalar ) rbox(p_ch3ooh) = gas(ngch3o2h) / factgas
gas_aqfrac_box(:) = 0.0
if (p_nh3 .le. numgas_aqfrac) &
gas_aqfrac_box(p_nh3 ) = gas_aqfrac_cmu(nga )
if (p_hno3 .le. numgas_aqfrac) &
gas_aqfrac_box(p_hno3 ) = gas_aqfrac_cmu(ngn )
if (p_hcl .le. numgas_aqfrac .and. p_hcl .gt. param_first_scalar) &
gas_aqfrac_box(p_hcl ) = gas_aqfrac_cmu(ngc )
if (p_sulf .le. numgas_aqfrac .and. p_sulf .gt. param_first_scalar) &
gas_aqfrac_box(p_sulf ) = gas_aqfrac_cmu(ng4 )
if (p_h2so4 .le. numgas_aqfrac .and. p_h2so4 .gt. param_first_scalar) &
gas_aqfrac_box(p_h2so4 ) = gas_aqfrac_cmu(ng4 )
if (p_hcho .le. numgas_aqfrac) &
gas_aqfrac_box(p_hcho ) = gas_aqfrac_cmu(nghcho )
if (p_ora1 .le. numgas_aqfrac .and. p_ora1 .gt. param_first_scalar) &
gas_aqfrac_box(p_ora1 ) = gas_aqfrac_cmu(nghcooh )
if (p_hcooh .le. numgas_aqfrac .and. p_hcooh .gt. param_first_scalar) &
gas_aqfrac_box(p_hcooh ) = gas_aqfrac_cmu(nghcooh )
if (p_so2 .le. numgas_aqfrac) &
gas_aqfrac_box(p_so2 ) = gas_aqfrac_cmu(ngso2 )
if (p_h2o2 .le. numgas_aqfrac) &
gas_aqfrac_box(p_h2o2 ) = gas_aqfrac_cmu(ngh2o2 )
if (p_o3 .le. numgas_aqfrac) &
gas_aqfrac_box(p_o3 ) = gas_aqfrac_cmu(ngo3 )
if (p_ho .le. numgas_aqfrac) &
gas_aqfrac_box(p_ho ) = gas_aqfrac_cmu(ngoh )
if (p_ho2 .le. numgas_aqfrac) &
gas_aqfrac_box(p_ho2 ) = gas_aqfrac_cmu(ngho2 )
if (p_no3 .le. numgas_aqfrac) &
gas_aqfrac_box(p_no3 ) = gas_aqfrac_cmu(ngno3 )
if (p_no .le. numgas_aqfrac) &
gas_aqfrac_box(p_no ) = gas_aqfrac_cmu(ngno )
if (p_no2 .le. numgas_aqfrac) &
gas_aqfrac_box(p_no2 ) = gas_aqfrac_cmu(ngno2 )
if (p_hono .le. numgas_aqfrac) &
gas_aqfrac_box(p_hono ) = gas_aqfrac_cmu(nghno2 )
if (p_pan .le. numgas_aqfrac) &
gas_aqfrac_box(p_pan ) = gas_aqfrac_cmu(ngpan )
if (p_ch3o2 .le. numgas_aqfrac) &
gas_aqfrac_box(p_ch3o2 ) = gas_aqfrac_cmu(ngch3o2 )
if (p_ch3oh .le. numgas_aqfrac) &
gas_aqfrac_box(p_ch3oh ) = gas_aqfrac_cmu(ngch3oh )
if (p_op1 .le. numgas_aqfrac .and. p_op1 .gt. param_first_scalar) &
gas_aqfrac_box(p_op1 ) = gas_aqfrac_cmu(ngch3o2h)
if (p_ch3ooh.le. numgas_aqfrac .and. p_ch3ooh .gt. param_first_scalar) &
gas_aqfrac_box(p_ch3ooh) = gas_aqfrac_cmu(ngch3o2h)
rbulk_cwaer(l_so4_aqyy,2) = aerosol(na4) / factaerso4
rbulk_cwaer(l_no3_aqyy,2) = aerosol(nan) / factaerno3
rbulk_cwaer(l_cl_aqyy, 2) = aerosol(nac) / factaercl
rbulk_cwaer(l_nh4_aqyy,2) = aerosol(naa) / factaernh4
rbulk_cwaer(l_na_aqyy, 2) = aerosol(nas) / factaerna
rbulk_cwaer(l_oin_aqyy,2) = aerosol(nar) / factaeroin
rbulk_cwaer(l_bc_aqyy, 2) = aerosol(nae) / factaerbc
rbulk_cwaer(l_oc_aqyy, 2) = aerosol(nao) / factaeroc
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
lunxx = 87
lunxx = -1
if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
if (lunxx .gt. 0) then
write(lunxx,*)
write(lunxx,*) 'interface_to_aqoperator1 - after call'
write(lunxx,*) 'gas (1=nh3, 2=hno3, 3=hcl, 4=h2so4, 11=so2, 12=h2o2, 18=o3)'
write(lunxx,9875) gas
write(lunxx,*) 'rbox(nh3, hno3, hcl, h2so4, so2, h2o2, o3)'
write(lunxx,9875) rbox(p_nh3), rbox(p_hno3), rbox(p_hcl), &
rbox(p_sulf), rbox(p_h2so4), rbox(p_so2), rbox(p_h2o2), rbox(p_o3)
write(lunxx,*) 'aerosol (1=na, 3=nh4, 4=no3, 5=cl, 6=so4, 8=ec, 9=oc, 10=crus)'
write(lunxx,9875) aerosol
write(lunxx,*) 'rbulk_cwaer (1=so4, 2=no3, 3-cl, 4=nh4, 5=na, 6=oin, 7=bc, 8=oc)'
write(lunxx,9875) rbulk_cwaer(:,2)
write(lunxx,*) 'ph_cmuaq_cur'
write(lunxx,9875) ph_cmuaq_cur
if (icase .ge. 10) then
write(*,*) &
'*** stopping in interface_to_aqop1 at icase =', icase
! stop
end if
end if
#endif
return
end subroutine interface_to_aqoperator1
!-----------------------------------------------------------------------
subroutine partition_cldwtr( &
rbox, fr_partit_cw, &
it, jt, kt, kpeg, mpeg, icase )
use module_state_description, only: &
param_first_scalar
use module_data_mosaic_asect, only: &
maxd_asize, maxd_atype, &
cw_phase, nsize_aer, ntype_aer, ncomp_aer, &
massptr_aer, numptr_aer, &
dens_aer, mw_aer, volumlo_sect, volumhi_sect
use module_data_mosaic_other, only: &
aboxtest_units_convert, cairclm, &
ktemp, l2maxd, ptotclm, rcldwtr_sub
implicit none
! subr arguments
integer, intent(in) :: it, jt, kt, kpeg, mpeg, icase
real, intent(inout), dimension( 1:l2maxd ) :: rbox
real, intent(inout), dimension( maxd_asize, maxd_atype ) :: &
fr_partit_cw
! local variables
integer :: iphase, isize, itype
integer :: jdone_mass, jdone_numb, jpos, jpos_mass, jpos_numb
integer :: l, ll, lunxx
integer :: p1st
real, parameter :: partit_wght_mass = 0.5
real :: dum, duma, dumb, dumc, dummass, dumnumb, dumvolu
real :: tmass, tnumb, umass, unumb, wmass, wnumb
real, dimension( maxd_asize, maxd_atype ) :: fmass, fnumb, xmass, xnumb
p1st = PARAM_FIRST_SCALAR
iphase = cw_phase
tmass = 0.0
tnumb = 0.0
umass = 0.0
unumb = 0.0
! compute
! xmass, xnumb = mass, number mixing ratio for a bin
! tmass, tnumb = sum over all bins of xmass, xnumb
! umass, unumb = max over all bins of xmass, xnumb
! set xmass, xnumb = 0.0 if bin mass, numb < 1.0e-37
! constrain xnumb so that mean particle volume is
! within bin boundaries
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
dummass = 0.0
dumvolu = 0.0
do ll = 1, ncomp_aer(itype)
l = massptr_aer(ll,isize,itype,iphase)
if (l .ge. p1st) then
dum = max( 0.0, rbox(l) )*mw_aer(ll,itype)
dummass = dummass + dum
dumvolu = dumvolu + dum/dens_aer(ll,itype)
end if
end do
l = numptr_aer(isize,itype,iphase)
dumnumb = max( 0.0, rbox(l) )
if (dumnumb .gt. dumvolu/volumlo_sect(isize,itype)) then
dumnumb = dumvolu/volumlo_sect(isize,itype)
rbox(l) = dumnumb
else if (dumnumb .lt. dumvolu/volumhi_sect(isize,itype)) then
dumnumb = dumvolu/volumhi_sect(isize,itype)
rbox(l) = dumnumb
end if
if (dummass .lt. 1.0e-37) dummass = 0.0
xmass(isize,itype) = dummass
if (dumnumb .lt. 1.0e-37) dumnumb = 0.0
xnumb(isize,itype) = dumnumb
tmass = tmass + xmass(isize,itype)
tnumb = tnumb + xnumb(isize,itype)
umass = max( umass, xmass(isize,itype) )
unumb = max( unumb, xnumb(isize,itype) )
end do
end do
! compute
! fmass, fnumb = fraction of total mass, number that is in a bin
! if tmass<1e-35 and umass>0, set fmass=1 for bin with largest xmass
! if tmass<1e-35 and umass=0, set fmass=0 for all
jdone_mass = 0
jdone_numb = 0
jpos_mass = 0
jpos_numb = 0
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
fmass(isize,itype) = 0.0
if (tmass .ge. 1.0e-35) then
fmass(isize,itype) = xmass(isize,itype)/tmass
else if (umass .gt. 0.0) then
if ( (jdone_mass .eq. 0) .and. &
(xmass(isize,itype) .eq. umass) ) then
jdone_mass = 1
fmass(isize,itype) = 1.0
end if
end if
if (fmass(isize,itype) .gt. 0) jpos_mass = jpos_mass + 1
fnumb(isize,itype) = 0.0
if (tnumb .ge. 1.0e-35) then
fnumb(isize,itype) = xnumb(isize,itype)/tnumb
else if (unumb .gt. 0.0) then
if ( (jdone_numb .eq. 0) .and. &
(xnumb(isize,itype) .eq. unumb) ) then
jdone_numb = 1
fnumb(isize,itype) = 1.0
end if
end if
if (fnumb(isize,itype) .gt. 0) jpos_numb = jpos_numb + 1
end do
end do
! if only 1 bin has fmass or fnumb > 0, set value to 1.0 exactly
if ((jpos_mass .eq. 1) .or. (jpos_numb .eq. 1)) then
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
if (jpos_mass .eq. 1) then
if (fmass(isize,itype) .gt. 0) fmass(isize,itype) = 1.0
end if
if (jpos_numb .eq. 1) then
if (fnumb(isize,itype) .gt. 0) fnumb(isize,itype) = 1.0
end if
end do
end do
end if
!
! compute fr_partit_cw as weighted average of fmass & fnumb, except
! if tmass<1e-35 and umass=0, use only fnumb
! if tnumb<1e-35 and unumb=0, use only fmass
! if tmass,tnumb<1e-35 and umass,unumb=0,
! set fr_partit_cw=1 for center bin of itype=1
!
fr_partit_cw(:,:) = 0.0
if ((jpos_mass .eq. 0) .and. (jpos_numb .eq. 0)) then
itype = 1
isize = (nsize_aer(itype)+1)/2
fr_partit_cw(isize,itype) = 1.0
else if (jpos_mass .eq. 0) then
fr_partit_cw(:,:) = fnumb(:,:)
else if (jpos_numb .eq. 0) then
fr_partit_cw(:,:) = fmass(:,:)
else
wmass = max( 0.0, min( 1.0, partit_wght_mass ) )
wnumb = 1.0 - wmass
fr_partit_cw(:,:) = wmass*fmass(:,:) + wnumb*fnumb(:,:)
jpos = 0
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
if (fr_partit_cw(isize,itype) .gt. 0.0) jpos = jpos + 1
end do
end do
! if only 1 bin has fr_partit_cw > 0, set value to 1.0 exactly
if (jpos .eq. 1) then
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
if (fr_partit_cw(isize,itype) .gt. 0.0) &
fr_partit_cw(isize,itype) = 1.0
end do
end do
end if
end if
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
! diagnostics when lunxx > 0
lunxx = 86
lunxx = -1
if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
if (lunxx .gt. 0) then
if (icase .le. 9) then
write(lunxx,9800)
write(lunxx,9800)
write(lunxx,9800) &
'partition_cldwtr - icase, jpos, jpos_mass, jpos_numb'
write(lunxx,9810) icase, jpos, jpos_mass, jpos_numb
write(lunxx,9800) 'tmass, umass, wmass'
write(lunxx,9820) tmass, umass, wmass
write(lunxx,9800) 'tnumb, unumb, wnumb'
write(lunxx,9820) tnumb, unumb, wnumb
write(lunxx,9800) 'xmass, fmass, xnumb, fnumb, fr_partit_cw'
duma = 0.0
dumb = 0.0
dumc = 0.0
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
write(lunxx,9820) xmass(isize,itype), fmass(isize,itype), &
xnumb(isize,itype), fnumb(isize,itype), &
fr_partit_cw(isize,itype)
duma = duma + fmass(isize,itype)
dumb = dumb + fnumb(isize,itype)
dumc = dumc + fr_partit_cw(isize,itype)
end do
end do
write(lunxx,9800) &
'sum_fmass-1.0, sum_fnumb-1.0, sum_fr_partit-1.0'
write(lunxx,9820) (duma-1.0), (dumb-1.0), (dumc-1.0)
if (icase .eq. -2) then
write(*,*) '*** stopping in partition_cldwtr at icase =', icase
! stop
end if
9800 format( a )
9810 format( 5i10 )
9820 format( 5(1pe10.2) )
end if
end if
#endif
return
end subroutine partition_cldwtr
!-----------------------------------------------------------------------
subroutine distribute_bulk_changes( &
rbox, rbox_sv1, fr_partit_cw, &
rbulk_cwaer, lptr_yyy_cwaer, &
it, jt, kt, kpeg, mpeg, icase )
use module_state_description, only: &
param_first_scalar
use module_data_mosaic_asect, only: &
maxd_asize, maxd_atype, &
cw_phase, nsize_aer, ntype_aer, &
lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
lptr_oin_aer, lptr_bc_aer, lptr_oc_aer
use module_data_mosaic_other, only: l2maxd, lunout, name
implicit none
! subr arguments
integer, intent(in) :: it, jt, kt, kpeg, mpeg, icase
integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
real, intent(inout), dimension( 1:l2maxd ) :: rbox, rbox_sv1
real, intent(in), dimension( maxd_asize, maxd_atype ) :: &
fr_partit_cw
real, intent(in), dimension( nyyy, 2 ) :: rbulk_cwaer
! local variables
integer :: iphase, isize, itype
integer :: idone, icount, ncount
integer :: jpos, jpos_sv
integer :: l, lunxx, lunxxaa, lunxxbb, lyyy
integer :: p1st
real :: duma, dumb, dumc
real :: fr, frsum_cur
real :: fr_cur(maxd_asize,maxd_atype)
real :: del_r_current, del_r_remain
real :: del_rbulk_cwaer(nyyy)
p1st = param_first_scalar
do lyyy = 1, nyyy
del_rbulk_cwaer(lyyy) = rbulk_cwaer(lyyy,2) - rbulk_cwaer(lyyy,1)
end do
iphase = cw_phase
jpos = 0
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
if (fr_partit_cw(isize,itype) .gt. 0) jpos = jpos + 1
end do
end do
jpos_sv = jpos
!
! distribution is trivial when only 1 bin has fr_partit_cw > 0
!
if (jpos_sv .eq. 1) then
do lyyy = 1, nyyy
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
fr = fr_partit_cw(isize,itype)
if (fr .eq. 1.0) then
l = lptr_yyy_cwaer(isize,itype,lyyy)
if (l .ge. p1st) rbox(l) = rbulk_cwaer(lyyy,2)
end if
end do
end do
end do
goto 7900
end if
do 3900 lyyy = 1, nyyy
!
! distribution is simple when del_rbulk_cwaer(lyyy) >= 0
!
if (del_rbulk_cwaer(lyyy) .eq. 0.0) then
goto 3900
else if (del_rbulk_cwaer(lyyy) .gt. 0.0) then
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
fr = fr_partit_cw(isize,itype)
if (fr .gt. 0.0) then
l = lptr_yyy_cwaer(isize,itype,lyyy)
if (l .ge. p1st) then
rbox(l) = rbox(l) + fr*del_rbulk_cwaer(lyyy)
end if
end if
end do
end do
goto 3900
end if
!
! distribution is complicated when del_rbulk_cwaer(lyyy) < 0,
! because you cannot produce any negative mixrats
!
del_r_remain = del_rbulk_cwaer(lyyy)
fr_cur(:,:) = fr_partit_cw(:,:)
ncount = max( 1, jpos_sv*2 )
icount = 0
! iteration loop
do while (icount .le. ncount)
icount = icount + 1
del_r_current = del_r_remain
jpos = 0
frsum_cur = 0.0
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
fr = fr_cur(isize,itype)
if (fr .gt. 0.0) then
l = lptr_yyy_cwaer(isize,itype,lyyy)
if (l .ge. p1st) then
duma = fr*del_r_current
dumb = rbox(l) + duma
if (dumb .gt. 0.0) then
jpos = jpos + 1
else if (dumb .eq. 0.0) then
fr_cur(isize,itype) = 0.0
else
duma = -rbox(l)
dumb = 0.0
fr_cur(isize,itype) = 0.0
end if
del_r_remain = del_r_remain - duma
rbox(l) = dumb
frsum_cur = frsum_cur + fr_cur(isize,itype)
else
fr_cur(isize,itype) = 0.0
end if
end if
end do ! isize = 1, nsize_aer
end do ! itype = 1, ntype_aer
! done if jpos = jpos_sv, because bins reached zero mixrat
if (jpos .eq. jpos_sv) then
idone = 1
! del_r_remain starts as negative, so done if non-negative
else if (del_r_remain .ge. 0.0) then
idone = 2
! del_r_remain starts as negative, so done if non-negative
else if (abs(del_r_remain) .le. 1.0e-7*abs(del_rbulk_cwaer(lyyy))) then
idone = 3
! done if all bins have fr_cur = 0
else if (frsum_cur .le. 0.0) then
idone = 4
! same thing basically
else if (jpos .le. 0) then
idone = 5
else
idone = 0
end if
! check for done, and (conditionally) print message
if (idone .gt. 0) then
lunxxaa = 6
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
lunxxaa = 86
#endif
if ((lunxxaa .gt. 0) .and. (icount .gt. (1+jpos_sv)/2)) then
write(lunxxaa,9800) &
'distribute_bulk_changes - icount>jpos_sv/2 - i,j,k'
write(lunxxaa,9810) it, jt, kt
write(lunxxaa,9800) 'icase, lyyy, idone, icount, jpos, jpos_sv'
write(lunxxaa,9810) icase, lyyy, idone, icount, jpos, jpos_sv
end if
goto 3900
end if
! rescale fr_cur for next iteration
fr_cur(:,:) = fr_cur(:,:)/frsum_cur
end do ! while (icount .le. ncount)
! icount > ncount, so print message
lunxxbb = 6
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
lunxxbb = 86
#endif
if (lunxxbb .gt. 0) then
write(lunxxbb,9800)
write(lunxxbb,9800) &
'distribute_bulk_changes - icount>ncount - i,j,k'
write(lunxxbb,9810) it, jt, kt
write(lunxxbb,9800) 'icase, lyyy, icount, ncount, jpos_sv, jpos'
write(lunxxbb,9810) icase, lyyy, icount, ncount, jpos_sv, jpos
write(lunxxbb,9800) 'rbulk_cwaer(1), del_rbulk_cwaer, del_r_remain, frsum_cur, (frsum_cur-1.0)'
write(lunxxbb,9820) rbulk_cwaer(lyyy,1), del_rbulk_cwaer(lyyy), &
del_r_remain, frsum_cur, (frsum_cur-1.0)
end if
9800 format( a )
9801 format( 3a )
9810 format( 7i10 )
9820 format( 7(1pe10.2) )
9840 format( 2i3, 5(1pe14.6) )
3900 continue
7900 continue
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
! diagnostics for testing
lunxx = 88
lunxx = -1
if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
if (lunxx .gt. 0) then
icount = 0
do lyyy = 1, nyyy
duma = del_rbulk_cwaer(lyyy)
if ( abs(duma) .gt. &
max( 1.0e-35, 1.0e-5*abs(rbulk_cwaer(lyyy,1)) ) ) then
icount = icount + 1
if (icount .eq. 1) write(lunxx,9800)
if (icount .eq. 1) write(lunxx,9800)
write(lunxx,9800)
write(lunxx,9801) 'distribute_bulk_changes - ',name(lptr_yyy_cwaer(1,1,lyyy)), ' - icase, lyyy'
write(lunxx,9810) icase, lyyy
write(lunxx,9800) ' tp sz rbox_sv1, rbox, del_rbox, del_rbox/del_rbulk_cwaer, (...-fr_partit_cw)'
write(lunxx,9840) 0, 0, rbulk_cwaer(lyyy,1:2), duma
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
l = lptr_yyy_cwaer(isize,itype,lyyy)
dumb = rbox(l) - rbox_sv1(l)
dumc = dumb/max( abs(duma), 1.0e-35 )
if (duma .lt. 0.0) dumc = -dumc
write(lunxx,9840) itype, isize, rbox_sv1(l), rbox(l), &
dumb, dumc, (dumc-fr_partit_cw(isize,itype))
end do
end do
end if
end do
if (icase .ge. 5) then
write(*,*) &
'*** stop in distribute_bulk_changes diags, icase =', icase
! stop
end if
end if
#endif
return
end subroutine distribute_bulk_changes
!-----------------------------------------------------------------------
subroutine cloudchem_apply_move_sections( &
rbox, rbox_sv1, &
it, jt, kt, kpeg, mpeg, icase )
use module_state_description, only: &
param_first_scalar
use module_data_mosaic_asect, only: &
msectional, &
maxd_asize, maxd_atype, &
cw_phase, nsize_aer, ntype_aer, ncomp_aer, &
massptr_aer, numptr_aer, mw_aer, dens_aer, &
lptr_so4_aer, lptr_no3_aer, lptr_cl_aer, lptr_co3_aer, &
lptr_msa_aer, lptr_nh4_aer, lptr_na_aer, lptr_ca_aer, &
lptr_oin_aer, lptr_bc_aer, lptr_oc_aer, &
drymass_aftgrow, drymass_pregrow, &
drydens_aftgrow, drydens_pregrow
use module_data_mosaic_other, only: l2maxd, name, rsub
use module_mosaic_movesect, only: move_sections
implicit none
! subr arguments
integer, intent(in) :: it, jt, kt, kpeg, mpeg, icase
real, intent(inout), dimension( 1:l2maxd ) :: rbox, rbox_sv1
! local variables
integer :: idum_msect
integer :: iphase, isize, itype
integer :: l, ll, lunxx
integer :: p1st
integer :: lptr_dum(maxd_asize,maxd_atype)
real :: densdefault
real :: dmaft, dmpre, dvaft, dvpre
real :: duma, dumb, dumc
real :: smallmassbb
p1st = param_first_scalar
iphase = cw_phase
!
! compute drymass before and after growth
! set drydens = -1.0, so it will be calculated
!
densdefault = 2.0
smallmassbb = 1.0e-30
do 1800 itype = 1, ntype_aer
do 1800 isize = 1, nsize_aer(itype)
dmaft = 0.0
dmpre = 0.0
dvaft = 0.0
dvpre = 0.0
do ll = 1, ncomp_aer(itype)
l = massptr_aer(ll,isize,itype,iphase)
if (l .ge. p1st) then
duma = mw_aer(ll,itype)
dmaft = dmaft + duma*rbox(l)
dmpre = dmpre + duma*rbox_sv1(l)
duma = duma/dens_aer(ll,itype)
dvaft = dvaft + duma*rbox(l)
dvpre = dvpre + duma*rbox_sv1(l)
end if
end do
drymass_aftgrow(isize,itype) = dmaft
drymass_pregrow(isize,itype) = dmpre
if (min(dmaft,dvaft) .le. smallmassbb) then
drydens_aftgrow(isize,itype) = densdefault
else
drydens_aftgrow(isize,itype) = dmaft/dvaft
end if
if (min(dmpre,dvpre) .le. smallmassbb) then
drydens_pregrow(isize,itype) = densdefault
else
drydens_pregrow(isize,itype) = dmpre/dvpre
end if
1800 continue
! apply move sections routine
! (and conditionally turn on movesect diagnostics)
idum_msect = msectional
!#if defined ( ccboxtest_box_testing_active )
#if defined ( temp_box_testing_active)
lunxx = 88
lunxx = -1
if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
if (lunxx .gt. 0) then
if (msectional .lt. 7000) msectional = msectional + 7000
end if
#endif
call move_sections( 2, it, jt, kpeg, mpeg )
msectional = idum_msect
!#if defined ( ccboxtest_box_testing_active)
#if defined ( temp_box_testing_active)
! diagnostics for testing
lunxx = 88
lunxx = -1
if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
if (lunxx .gt. 0) then
do ll = 1, 4
if (ll .eq. 1) then
lptr_dum(:,:) = lptr_so4_aer(:,:,iphase)
else if (ll .eq. 2) then
lptr_dum(:,:) = lptr_nh4_aer(:,:,iphase)
else if (ll .eq. 3) then
lptr_dum(:,:) = lptr_oin_aer(:,:,iphase)
else if (ll .eq. 4) then
lptr_dum(:,:) = numptr_aer(:,:,iphase)
end if
if (ll .eq. 1) write(lunxx,9800)
if (ll .eq. 1) write(lunxx,9800)
write(lunxx,9800)
write(lunxx,9800)
write(lunxx,9801) 'cloudchem_apply_move_sections - ', &
name(lptr_dum(1,1)), ' - icase, ll'
write(lunxx,9810) icase, ll
write(lunxx,9800) ' tp sz rbox_sv1, rbox, rsub'
do itype = 1, ntype_aer
do isize = 1, nsize_aer(itype)
l = lptr_dum(isize,itype)
dumb = rbox(l) - rbox_sv1(l)
dumc = dumb/max( abs(duma), 1.0e-35 )
if (duma .lt. 0.0) dumc = -dumc
write(lunxx,9840) itype, isize, rbox_sv1(l), rbox(l), &
rsub(l,kpeg,mpeg)
end do
end do
end do
if (icase .ge. 5) then
write(*,*) &
'*** stop in cloudchem_apply_move_sections diags, icase =', &
icase
! stop
end if
end if
9800 format( a )
9801 format( 3a )
9810 format( 7i10 )
9820 format( 7(1pe10.2) )
9840 format( 2i3, 5(1pe14.6) )
#endif
return
end subroutine cloudchem_apply_move_sections
!-----------------------------------------------------------------------
subroutine mosaic_cloudchem_dumpaa( &
id, ktau, ktauc, dtstepc, config_flags, &
p_phy, t_phy, rho_phy, alt, &
cldfra, ph_no2, &
moist, chem, &
gas_aqfrac, numgas_aqfrac, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
qcldwtr_cutoff, &
itcur, jtcur, ktcur )
use module_state_description, only: &
num_moist, num_chem, p_qc
use module_configure, only: grid_config_rec_type
use module_data_mosaic_asect
use module_data_mosaic_other, only: k_pegbegin, name
use module_mosaic_driver, only: mapaer_tofrom_host
implicit none
! subr arguments
integer, intent(in) :: &
id, ktau, ktauc, &
numgas_aqfrac, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
itcur, jtcur, ktcur
! id - domain index
! ktau - time step number
! ktauc - gas and aerosol chemistry time step number
! numgas_aqfrac - last dimension of gas_aqfrac
! [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for 'domain'
! [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for 'memory'
! Most arrays that are arguments to chem_driver
! are dimensioned with these spatial indices.
! [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for 'tile'
! chem_driver and routines under it do calculations
! over these spatial indices.
type(grid_config_rec_type), intent(in) :: config_flags
! config_flags - configuration and control parameters
real, intent(in) :: &
dtstepc, qcldwtr_cutoff
! dtstepc - time step for gas and aerosol chemistry(s)
real, intent(in), &
dimension( ims:ime, kms:kme, jms:jme ) :: &
p_phy, t_phy, rho_phy, alt, cldfra, ph_no2
! p_phy - air pressure (Pa)
! t_phy - temperature (K)
! rho_phy - moist air density (kg/m^3)
! alt - dry air specific volume (m^3/kg)
! cldfra - cloud fractional area (0-1)
! ph_no2 - no2 photolysis rate (1/min)
real, intent(in), &
dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
moist
! moist - mixing ratios of moisture species (water vapor,
! cloud water, ...) (kg/kg for mass species, #/kg for number species)
real, intent(inout), &
dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
chem
! chem - mixing ratios of trace gas and aerosol species (ppm for gases,
! ug/kg for aerosol mass species, #/kg for aerosol number species)
real, intent(inout), &
dimension( ims:ime, kms:kme, jms:jme, numgas_aqfrac ) :: &
gas_aqfrac
! gas_aqfrac - fraction (0-1) of gas that is dissolved in cloud water
! local variables
integer :: it, jt, kt, l, ll, n
integer :: isize, itype
real :: dumai, dumcw
real :: qcldwtr
it = itcur
jt = jtcur
kt = ktcur
write(*,*)
write(*,*)
write(*,*)
write(*,9100)
write(*,9102) ktau, it, jt, kt
9100 format( 7('----------') )
9102 format( &
'mosaic_cloudchem_dumpaa - ktau, i, j, k =', 4i5 )
itype = 1
do 2900 isize = 1, nsize_aer(itype)
write(*,9110) isize
9110 format( / 'isize =', i3 / &
' k cldwtr mass-ai numb-ai mass-cw numb-cw' )
do 2800 kt = kte, kts, -1
dumai = 0.0
dumcw = 0.0
do ll = 1, ncomp_aer(itype)
l = massptr_aer(ll,isize,itype,1)
dumai = dumai + chem(it,kt,jt,l)
l = massptr_aer(ll,isize,itype,2)
dumcw = dumcw + chem(it,kt,jt,l)
end do
write(*,9120) kt, &
moist(it,kt,jt,p_qc), &
dumai, chem(it,kt,jt,numptr_aer(isize,itype,1)), &
dumcw, chem(it,kt,jt,numptr_aer(isize,itype,2))
9120 format( i3, 1p, e10.2, 2(3x, 2e10.2) )
2800 continue
2900 continue
write(*,*)
write(*,9100)
write(*,*)
! map from wrf-chem 3d arrays to pegasus clm & sub arrays
kt = ktcur
if ((ktau .eq. 30) .and. (it .eq. 23) .and. &
(jt .eq. 1) .and. (kt .eq. 11)) then
qcldwtr = moist(it,kt,jt,p_qc)
write(*,*)
write(*,*)
write(*,9102) ktau, it, jt, kt
write(*,*)
write( *, '(3(1pe10.2,3x,a))' ) &
(chem(it,kt,jt,l), name(l)(1:10), l=1,num_chem)
write(*,*)
write( *, '(3(1pe10.2,3x,a))' ) &
p_phy(it,kt,jt), 'p_phy ', &
t_phy(it,kt,jt), 't_phy ', &
rho_phy(it,kt,jt), 'rho_phy ', &
alt(it,kt,jt), 'alt ', &
qcldwtr, 'qcldwtr ', &
qcldwtr_cutoff, 'qcldwtrcut'
write(*,*)
write(*,9100)
write(*,*)
end if
return
end subroutine mosaic_cloudchem_dumpaa
!-----------------------------------------------------------------------
end module module_mosaic_cloudchem