!********************************************************************************** ! 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_sorgam_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 ! "negligible volume" = (1 #/kg ~= 1e-6 #/cm3) * ((dp=1e-6 cm)**3 * pi/6) real, parameter :: smallvolaa = 0.5e-18 contains !----------------------------------------------------------------------- subroutine sorgam_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_sorgam, only: cw_phase, nphase_aer 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, l integer :: icase integer :: igaschem_onoff, iphotol_onoff, iradical_onoff real :: rbox(num_chem) 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 *, '*** sorgam_cloudchem_driver - cw_phase not active' return end if print 93010, 'entering sorgam_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 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 sorgam_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 rbox(1:num_chem) = chem(it,kt,jt,1:num_chem) ! make call '1box' cloudchem routine ! print 93010, 'calling sorgam_cloudchem_1 at ijk =', it, jt, kt call sorgam_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, icase, & rbox, qcldwtr, & t_phy(it,kt,jt), p_phy(it,kt,jt), rho_phy(it,kt,jt), & config_flags) ! map back to wrf-chem 3d arrays chem(it,kt,jt,1:num_chem) = rbox(1:num_chem) gas_aqfrac(it,kt,jt,:) = gas_aqfrac_box(:) 3800 continue 3910 continue 3920 continue print 93010, 'leaving sorgam_cloudchem_driver - ktau =', ktau, icase 93010 format( a, 8(1x,i6) ) return end subroutine sorgam_cloudchem_driver !----------------------------------------------------------------------- subroutine sorgam_cloudchem_1box( & id, ktau, ktauc, dtstepc, & iphotol_onoff, iradical_onoff, & photol_no2_box, & ph_aq_box, gas_aqfrac_box, & numgas_aqfrac, it, jt, kt, icase, & rbox, qcw_box, temp_box, pres_box, rho_box, & config_flags ) use module_configure, only: grid_config_rec_type use module_state_description, only: & num_moist, num_chem use module_data_sorgam, only: & msectional, maxd_asize, maxd_atype, & cw_phase, nsize_aer, ntype_aer, do_cloudchem_aer, & lptr_so4_aer, lptr_no3_aer, lptr_nh4_aer, & lptr_orgpa_aer, lptr_ec_aer, lptr_p25_aer, & lptr_cl_aer, lptr_na_aer use module_data_cmu_bulkaqchem, only: & meqn1max implicit none ! subr arguments type(grid_config_rec_type), intent(in) :: config_flags integer, intent(in) :: & id, ktau, ktauc, & numgas_aqfrac, it, jt, kt, & icase, iphotol_onoff, iradical_onoff real, intent(in) :: & dtstepc, photol_no2_box, & qcw_box, & ! cloud water (kg/kg) temp_box, & ! air temp (K) pres_box, & ! air pres (Pa) rho_box ! air dens (kg/m3) real, intent(inout) :: ph_aq_box real, intent(inout), dimension( num_chem ) :: rbox ! rbox has same units as chem [gas = ppmv, AP mass = ug/kg, AP number = #/kg] 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_sv1(num_chem) real :: rbulk_cwaer(nyyy,2) real, dimension( maxd_asize, maxd_atype ) :: fr_partit_cw real, dimension( 2, 3 ) :: xvol_old ! ! 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_nh4_aqyy) = lptr_nh4_aer(:,:,iphase) lptr_yyy_cwaer(:,:,l_oin_aqyy) = lptr_p25_aer(:,:,iphase) lptr_yyy_cwaer(:,:,l_bc_aqyy ) = lptr_ec_aer( :,:,iphase) lptr_yyy_cwaer(:,:,l_oc_aqyy ) = lptr_orgpa_aer( :,:,iphase) ! lptr_yyy_cwaer(:,:,l_cl_aqyy ) = -999888777 ! lptr_yyy_cwaer(:,:,l_na_aqyy ) = -999888777 lptr_yyy_cwaer(:,:,l_cl_aqyy ) = lptr_cl_aer(:,:,iphase) lptr_yyy_cwaer(:,:,l_na_aqyy ) = lptr_na_aer(:,:,iphase) ! ! ! 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 rbox_sv1(:) = rbox(:) gas_aqfrac_box(:) = 0.0 ! make call to interface_to_aqoperator1 call sorgam_interface_to_aqoperator1( & istat_aqop, & dtstepc, & rbox, gas_aqfrac_box, & qcw_box, temp_box, pres_box, rho_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, ktau, icase_in, & config_flags) 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 sorgam_partition_cldwtr( & rbox, fr_partit_cw, xvol_old, & id, it, jt, kt, icase_in ) ! ! ! distribute changes in bulk cloud-water composition among size bins ! ! call sorgam_distribute_bulk_changes( & rbox, rbox_sv1, fr_partit_cw, & rbulk_cwaer, lptr_yyy_cwaer, & id, it, jt, kt, icase_in ) ! ! do move-sections ! if (msectional .lt. 1000000000) then call sorgam_cloudchem_apply_mode_transfer( & rbox, rbox_sv1, xvol_old, & id, it, jt, kt, icase_in ) end if return end subroutine sorgam_cloudchem_1box !----------------------------------------------------------------------- subroutine sorgam_interface_to_aqoperator1( & istat_aqop, & dtstepc, & rbox, gas_aqfrac_box, & qcw_box, temp_box, pres_box, rho_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, ktau, icase, & config_flags ) use module_configure, only: grid_config_rec_type use module_state_description, only: & num_chem, param_first_scalar, p_qc, & p_nh3, p_hno3, p_hcl, p_sulf, 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_form, p_facd, p_oh, p_meo2, p_meoh, p_mepx, & CB05_SORG_AQ_KPP 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_sorgam, only: & maxd_asize, maxd_atype, & cw_phase, nsize_aer, ntype_aer, do_cloudchem_aer, & lptr_so4_aer, lptr_no3_aer, lptr_nh4_aer, & lptr_orgpa_aer, lptr_ec_aer, lptr_p25_aer, & lptr_cl_aer, lptr_na_aer implicit none ! subr arguments type(grid_config_rec_type), intent(in) :: config_flags integer, intent(in) :: & iradical_in, idecomp_hmsa_hso5, & numgas_aqfrac, id, it, jt, kt, 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, & qcw_box, temp_box, pres_box, rho_box real, intent(inout) :: ph_cmuaq_cur real, intent(inout), dimension( num_chem ) :: rbox ! ppm or ug/kg real, intent(inout), dimension( numgas_aqfrac ) :: gas_aqfrac_box real, intent(inout), dimension( nyyy, 2 ) :: rbulk_cwaer real, intent(inout), dimension( meqn1max ) :: yaq_beg, yaq_end ! local variables integer :: i, iphase, isize, itype integer :: iaq, istat_fatal, istat_warn integer :: l, lyyy integer :: p1st #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) integer :: lunxx #endif 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 pascals] * factpatm factpatm = 1.0/1.01325e5 ! [cldwtr in g-h2o/m3-air] = [cldwtr in kg-h2o/kg-air] * factlwc factlwc = 1.0e3*rho_box ! [aq photolysis rate scaling factor in --] = [jno2 in 1/min] * factphoto factphoto = 1.6 ! [gas in ppm] = [gas in ppm] * factgas factgas = 1.0 ! [aerosol in ug/m3-air] = [aerosol in ug/kg-air] * factaer dum = rho_box factaerso4 = dum factaerno3 = dum factaercl = dum factaernh4 = dum factaerna = dum factaeroin = dum factaeroc = dum factaerbc = dum #if defined ( ccboxtest_box_testing_active) ! 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 #endif ! ! map from rbox to gas,aerosol ! temp = temp_box lwc = qcw_box * factlwc p_atm = pres_box * 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 ! for made-sorgam, set p_atm so that cmu code's (cair*28.966e3) ! will match rho_box p_atm = (rho_box/28.966)*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 if (p_nh3 >= p1st) gas(nga ) = rbox(p_nh3 )*factgas if (p_hno3 >= p1st) gas(ngn ) = rbox(p_hno3 )*factgas if (p_hcl >= p1st) gas(ngc ) = rbox(p_hcl )*factgas if (p_sulf >= p1st) gas(ng4 ) = rbox(p_sulf )*factgas ! if (p_hcho >= p1st) gas(nghcho ) = rbox(p_hcho )*factgas ! if (p_ora1 >= p1st) gas(nghcooh ) = rbox(p_ora1 )*factgas if (p_so2 >= p1st) gas(ngso2 ) = rbox(p_so2 )*factgas if (p_h2o2 >= p1st) gas(ngh2o2 ) = rbox(p_h2o2 )*factgas if (p_o3 >= p1st) gas(ngo3 ) = rbox(p_o3 )*factgas ! if (p_ho >= p1st) gas(ngoh ) = rbox(p_ho )*factgas if (p_ho2 >= p1st) gas(ngho2 ) = rbox(p_ho2 )*factgas if (p_no3 >= p1st) gas(ngno3 ) = rbox(p_no3 )*factgas if (p_no >= p1st) gas(ngno ) = rbox(p_no )*factgas if (p_no2 >= p1st) gas(ngno2 ) = rbox(p_no2 )*factgas if (p_hono >= p1st) gas(nghno2 ) = rbox(p_hono )*factgas if (p_pan >= p1st) gas(ngpan ) = rbox(p_pan )*factgas ! if (p_ch3o2 >= p1st) gas(ngch3o2 ) = rbox(p_ch3o2)*factgas ! if (p_ch3oh >= p1st) gas(ngch3oh ) = rbox(p_ch3oh)*factgas ! if (p_op1 >= p1st) gas(ngch3o2h) = rbox(p_op1 )*factgas if ((config_flags%chem_opt == CB05_SORG_AQ_KPP)) then if (p_form >= p1st) gas(nghcho ) = rbox(p_form )*factgas if (p_facd >= p1st) gas(nghcooh ) = rbox(p_facd )*factgas if (p_oh >= p1st) gas(ngoh ) = rbox(p_oh )*factgas if (p_meo2 >= p1st) gas(ngch3o2 ) = rbox(p_meo2 )*factgas if (p_meoh >= p1st) gas(ngch3oh ) = rbox(p_meoh )*factgas if (p_mepx >= p1st) gas(ngch3o2h) = rbox(p_mepx )*factgas else if (p_hcho >= p1st) gas(nghcho ) = rbox(p_hcho )*factgas if (p_ora1 >= p1st) gas(nghcooh ) = rbox(p_ora1 )*factgas if (p_ho >= p1st) gas(ngoh ) = rbox(p_ho )*factgas if (p_ch3o2 >= p1st) gas(ngch3o2 ) = rbox(p_ch3o2)*factgas if (p_ch3oh >= p1st) gas(ngch3oh ) = rbox(p_ch3oh)*factgas if (p_op1 >= p1st) gas(ngch3o2h) = rbox(p_op1 )*factgas endif ! 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) if ( .not. do_cloudchem_aer(isize,itype) ) cycle 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 ) || defined ( cctemp_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, ktau' write(lunxx,9870) it, jt, kt, 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,*) 'pres_box, rho_box, qcw_box, dt_sec' write(lunxx,9875) pres_box, rho_box, qcw_box, & (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_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 .le. -5) then write(*,*) & '*** stopping in interface_to_aqop1 at icase =', icase call wrf_error_fatal('*** stopping in interface_to_aqop1') 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,*) & '*** sorgam_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,*) & '*** sorgam_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,gas_aqfrac_box] to [rbox,gas_aqfrac_box] ! gas_aqfrac_box(:) = 0.0 if (p_nh3 >= p1st) then rbox(p_nh3 ) = gas(nga )/factgas if (p_nh3 <= numgas_aqfrac) & gas_aqfrac_box(p_nh3 ) = gas_aqfrac_cmu(nga ) end if if (p_hno3 >= p1st) then rbox(p_hno3 ) = gas(ngn )/factgas if (p_hno3 <= numgas_aqfrac) & gas_aqfrac_box(p_hno3 ) = gas_aqfrac_cmu(ngn ) end if if (p_hcl >= p1st) then rbox(p_hcl ) = gas(ngc )/factgas if (p_hcl <= numgas_aqfrac) & gas_aqfrac_box(p_hcl ) = gas_aqfrac_cmu(ngc ) end if if (p_sulf >= p1st) then rbox(p_sulf ) = gas(ng4 )/factgas if (p_sulf <= numgas_aqfrac) & gas_aqfrac_box(p_sulf ) = gas_aqfrac_cmu(ng4 ) end if ! if (p_hcho >= p1st) then ! rbox(p_hcho ) = gas(nghcho )/factgas ! if (p_hcho <= numgas_aqfrac) & ! gas_aqfrac_box(p_hcho ) = gas_aqfrac_cmu(nghcho ) ! end if ! if (p_ora1 >= p1st) then ! rbox(p_ora1 ) = gas(nghcooh )/factgas ! if (p_ora1 <= numgas_aqfrac) & ! gas_aqfrac_box(p_ora1 ) = gas_aqfrac_cmu(nghcooh ) ! end if if (p_so2 >= p1st) then rbox(p_so2 ) = gas(ngso2 )/factgas if (p_so2 <= numgas_aqfrac) & gas_aqfrac_box(p_so2 ) = gas_aqfrac_cmu(ngso2 ) end if if (p_h2o2 >= p1st) then rbox(p_h2o2 ) = gas(ngh2o2 )/factgas if (p_h2o2 <= numgas_aqfrac) & gas_aqfrac_box(p_h2o2 ) = gas_aqfrac_cmu(ngh2o2 ) end if if (p_o3 >= p1st) then rbox(p_o3 ) = gas(ngo3 )/factgas if (p_o3 <= numgas_aqfrac) & gas_aqfrac_box(p_o3 ) = gas_aqfrac_cmu(ngo3 ) end if ! if (p_ho >= p1st) then ! rbox(p_ho ) = gas(ngoh )/factgas ! if (p_ho <= numgas_aqfrac) & ! gas_aqfrac_box(p_ho ) = gas_aqfrac_cmu(ngoh ) ! end if if (p_ho2 >= p1st) then rbox(p_ho2 ) = gas(ngho2 )/factgas if (p_ho2 <= numgas_aqfrac) & gas_aqfrac_box(p_ho2 ) = gas_aqfrac_cmu(ngho2 ) end if if (p_no3 >= p1st) then rbox(p_no3 ) = gas(ngno3 )/factgas if (p_no3 <= numgas_aqfrac) & gas_aqfrac_box(p_no3 ) = gas_aqfrac_cmu(ngno3 ) end if if (p_no >= p1st) then rbox(p_no ) = gas(ngno )/factgas if (p_no <= numgas_aqfrac) & gas_aqfrac_box(p_no ) = gas_aqfrac_cmu(ngno ) end if if (p_no2 >= p1st) then rbox(p_no2 ) = gas(ngno2 )/factgas if (p_no2 <= numgas_aqfrac) & gas_aqfrac_box(p_no2 ) = gas_aqfrac_cmu(ngno2 ) end if if (p_hono >= p1st) then rbox(p_hono ) = gas(nghno2 )/factgas if (p_hono <= numgas_aqfrac) & gas_aqfrac_box(p_hono ) = gas_aqfrac_cmu(nghno2 ) end if if (p_pan >= p1st) then rbox(p_pan ) = gas(ngpan )/factgas if (p_pan <= numgas_aqfrac) & gas_aqfrac_box(p_pan ) = gas_aqfrac_cmu(ngpan ) end if ! if (p_ch3o2 >= p1st) then ! rbox(p_ch3o2) = gas(ngch3o2 )/factgas ! if (p_ch3o2 <= numgas_aqfrac) & ! gas_aqfrac_box(p_ch3o2 ) = gas_aqfrac_cmu(ngch3o2 ) ! end if ! if (p_ch3oh >= p1st) then ! rbox(p_ch3oh) = gas(ngch3oh )/factgas ! if (p_ch3oh <= numgas_aqfrac) & ! gas_aqfrac_box(p_ch3oh ) = gas_aqfrac_cmu(ngch3oh ) ! end if ! if (p_op1 >= p1st) then ! rbox(p_op1 ) = gas(ngch3o2h)/factgas ! if (p_op1 <= numgas_aqfrac) & ! gas_aqfrac_box(p_op1 ) = gas_aqfrac_cmu(ngch3o2h) ! end if ! if ( (config_flags%chem_opt == CB05_SORG_AQ_KPP) ) then if (p_form >= p1st) then rbox(p_form ) = gas(nghcho )/factgas if (p_form <= numgas_aqfrac) & gas_aqfrac_box(p_form ) = gas_aqfrac_cmu(nghcho ) end if if (p_facd >= p1st) then rbox(p_facd ) = gas(nghcooh )/factgas if (p_facd <= numgas_aqfrac) & gas_aqfrac_box(p_facd ) = gas_aqfrac_cmu(nghcooh ) end if if (p_oh >= p1st) then rbox(p_oh ) = gas(ngoh )/factgas if (p_oh <= numgas_aqfrac) & gas_aqfrac_box(p_oh ) = gas_aqfrac_cmu(ngoh ) end if if (p_meo2 >= p1st) then rbox(p_meo2 ) = gas(ngch3o2 )/factgas if (p_meo2 <= numgas_aqfrac) & gas_aqfrac_box(p_meo2 ) = gas_aqfrac_cmu(ngch3o2 ) end if if (p_meoh >= p1st) then rbox(p_meoh ) = gas(ngch3oh )/factgas if (p_meoh <= numgas_aqfrac) & gas_aqfrac_box(p_meoh ) = gas_aqfrac_cmu(ngch3oh ) end if if (p_mepx >= p1st) then rbox(p_mepx ) = gas(ngch3o2h)/factgas if (p_mepx <= numgas_aqfrac) & gas_aqfrac_box(p_mepx ) = gas_aqfrac_cmu(ngch3o2h) end if else if (p_hcho >= p1st) then rbox(p_hcho ) = gas(nghcho )/factgas if (p_hcho <= numgas_aqfrac) & gas_aqfrac_box(p_hcho ) = gas_aqfrac_cmu(nghcho ) end if if (p_ora1 >= p1st) then rbox(p_ora1 ) = gas(nghcooh )/factgas if (p_ora1 <= numgas_aqfrac) & gas_aqfrac_box(p_ora1 ) = gas_aqfrac_cmu(nghcooh ) end if if (p_ho >= p1st) then rbox(p_ho ) = gas(ngoh )/factgas if (p_ho <= numgas_aqfrac) & gas_aqfrac_box(p_ho ) = gas_aqfrac_cmu(ngoh ) end if if (p_ch3o2 >= p1st) then rbox(p_ch3o2) = gas(ngch3o2 )/factgas if (p_ch3o2 <= numgas_aqfrac) & gas_aqfrac_box(p_ch3o2 ) = gas_aqfrac_cmu(ngch3o2 ) end if if (p_ch3oh >= p1st) then rbox(p_ch3oh) = gas(ngch3oh )/factgas if (p_ch3oh <= numgas_aqfrac) & gas_aqfrac_box(p_ch3oh ) = gas_aqfrac_cmu(ngch3oh ) end if if (p_op1 >= p1st) then rbox(p_op1 ) = gas(ngch3o2h)/factgas if (p_op1 <= numgas_aqfrac) & gas_aqfrac_box(p_op1 ) = gas_aqfrac_cmu(ngch3o2h) end if end if 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 ) || defined ( cctemp_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_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 .le. -5) then write(*,*) & '*** stopping in interface_to_aqop1 at icase =', icase call wrf_error_fatal('*** stopping in interface_to_aqop1') end if end if #endif return end subroutine sorgam_interface_to_aqoperator1 !----------------------------------------------------------------------- subroutine sorgam_partition_cldwtr( & rbox, fr_partit_cw, xvol_old, & id, it, jt, kt, icase ) use module_state_description, only: & param_first_scalar, num_chem use module_data_sorgam, only: & maxd_asize, maxd_atype, & ai_phase, cw_phase, nsize_aer, ntype_aer, ncomp_aer, & do_cloudchem_aer, massptr_aer, numptr_aer, & dens_aer, sigmag_aer, & dcen_sect, dlo_sect, dhi_sect, & volumcen_sect, volumlo_sect, volumhi_sect implicit none ! subr arguments integer, intent(in) :: id, it, jt, kt, icase real, intent(inout), dimension( 1:num_chem ) :: rbox real, intent(inout), dimension( maxd_asize, maxd_atype ) :: & fr_partit_cw real, intent(inout), dimension( 2, 3 ) :: xvol_old ! local variables integer :: isize, itype integer :: jdone_mass, jdone_numb, jpos, jpos_mass, jpos_numb integer :: l, ll integer :: p1st #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) integer :: lunxx #endif real, parameter :: partit_wght_mass = 0.5 real :: tmpa, tmpb, tmpc real :: tmp_cwvolfrac, tmp_lnsg real :: tmass, tnumb, umass, unumb, wmass, wnumb real :: xmass_c, xmass_a, xmass_t, xvolu_c, xvolu_a, xvolu_t real :: xnumb_c1, xnumb_a1, xnumb_t1, xnumb_c2, xnumb_a2, xnumb_t2 real, dimension( maxd_asize, maxd_atype ) :: fmass, fnumb, xmass, xnumb, xnumbsv p1st = PARAM_FIRST_SCALAR 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 ! for made-sorgam, x/t/umass are g/kg-air, x/t/unumb are #/kg-air do itype = 1, ntype_aer do isize = 1, nsize_aer(itype) if ( .not. do_cloudchem_aer(isize,itype) ) cycle xmass_c = 0.0 xvolu_c = 0.0 xvolu_a = 0.0 do ll = 1, ncomp_aer(itype) l = massptr_aer(ll,isize,itype,cw_phase) if (l .ge. p1st) then tmpa = max( 0.0, rbox(l) )*1.0e-6 xmass_c = xmass_c + tmpa xvolu_c = xvolu_c + tmpa/dens_aer(ll,itype) end if l = massptr_aer(ll,isize,itype,ai_phase) if (l .ge. p1st) then tmpa = max( 0.0, rbox(l) )*1.0e-6 xvolu_a = xvolu_a + tmpa/dens_aer(ll,itype) end if end do xnumb_c1 = max( 0.0, rbox(numptr_aer(isize,itype,cw_phase)) ) xnumb_a1 = max( 0.0, rbox(numptr_aer(isize,itype,ai_phase)) ) xnumbsv(isize,itype) = xnumb_c1 xnumb_t1 = xnumb_a1 + xnumb_c1 xvolu_t = xvolu_a + xvolu_c ! do "bounding" activated+interstitial combined number ! and calculate dgnum for activated+interstitial combined if (xvolu_t < smallvolaa) then xnumb_t2 = xvolu_t/volumcen_sect(isize,itype) else if (xnumb_t1 < xvolu_t/volumhi_sect(isize,itype)) then xnumb_t2 = xvolu_t/volumhi_sect(isize,itype) else if (xnumb_t1 > xvolu_t/volumlo_sect(isize,itype)) then xnumb_t2 = xvolu_t/volumlo_sect(isize,itype) else xnumb_t2 = xnumb_t1 end if ! do "bounding" of activated number ! tmp_cwvolfrac = (cw volume)/(ai + cw volume) tmp_cwvolfrac = xvolu_c/max(xvolu_t,1.e-30) tmp_lnsg = log(sigmag_aer(isize,itype)) if ((xvolu_c < smallvolaa) .or. (tmp_cwvolfrac < 1.0e-10)) then ! for very small cw volume or volume fraction, ! use (ai+cw number)*(cw volume fraction) xnumb_c2 = xnumb_t2*tmp_cwvolfrac tmpa = -7.0 ; tmpb = -7.0 ; tmpc = -7.0 else ! tmpa is value of (ln(dpcut)-ln(dgvol))/ln(sigmag) for which ! "norm01 upper tail" cummulative pdf is equal to tmp_cwvolfrac tmpa = norm01_uptail_inv( tmp_cwvolfrac ) ! tmpb is corresponding value of (ln(dp)-ln(dgnum))/ln(sigmag) tmpb = tmpa + 3.0*tmp_lnsg ! tmpc is corresponding "norm01 upper tail" cummulative pdf tmpc = norm01_uptail( tmpb ) ! minimum number of activated particles occurs when ! activated particles are dp>dpcut and interstitial are dp 0 lunxx = 86 lunxx = -1 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171 if (lunxx > 0) then if ((isize == 1) .and. (itype == 1)) write(lunxx,'(a)') write(lunxx,'(a,3i4,4x,2i4,2x,l2)') 'partition-cw i,j,k, is,it', & it, jt, kt, isize, itype write(lunxx,'(a,1p,5e12.4)') 'cw vol, num old/adj ', & xvolu_c, xnumb_c1, xnumb_c2 write(lunxx,'(a,1p,5e12.4)') 'ai vol, num old/adj ', & xvolu_a, xnumb_a1, xnumb_a2 write(lunxx,'(a,1p,5e12.4)') 'a+c vol, num old/adj ', & xvolu_t, xnumb_t1, xnumb_t2 write(lunxx,'(a,1p,5e12.4)') 'lnsg, cwvolfr, cwnumfr 1/2', & tmp_lnsg, tmp_cwvolfrac, & xnumb_c1/max(xnumb_t1,1.0e-30), & xnumb_c2/max(xnumb_t2,1.0e-30) write(lunxx,'(a,1p,5e12.4)') 'tmpa/b/c ', & tmpa, tmpb, tmpc write(lunxx,'(a,1p,5e12.4)') 'dlo, dcen, dhi_sect ', & dlo_sect(isize,itype), dcen_sect(isize,itype), & dhi_sect(isize,itype) write(lunxx,'(a,1p,5e12.4)') 'vlo, vcen, vhi_sect ', & volumlo_sect(isize,itype), volumcen_sect(isize,itype), & volumhi_sect(isize,itype) end if #endif if (xmass_c .lt. 1.0e-37) xmass_c = 0.0 xmass(isize,itype) = xmass_c if (xnumb_c2 .lt. 1.0e-37) xnumb_c2 = 0.0 xnumb(isize,itype) = xnumb_c2 xnumbsv(isize,itype) = xnumb_c1 tmass = tmass + xmass(isize,itype) tnumb = tnumb + xnumb(isize,itype) umass = max( umass, xmass(isize,itype) ) unumb = max( unumb, xnumb(isize,itype) ) if ((itype == 1) .and. (isize <= 2)) then xvol_old(isize,1) = xvolu_c xvol_old(isize,2) = xvolu_a xvol_old(isize,3) = xvolu_t end if 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) if ( .not. do_cloudchem_aer(isize,itype) ) cycle 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 ( .not. do_cloudchem_aer(isize,itype) ) cycle 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 ( .not. do_cloudchem_aer(isize,itype) ) cycle 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 ( .not. do_cloudchem_aer(isize,itype) ) cycle 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 ) || defined ( cctemp_testing_active ) ! diagnostics when lunxx > 0 lunxx = 86 lunxx = -1 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171 ! if (icase .gt. 9) lunxx = -1 if (lunxx .gt. 0) then 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_orig/adj, fnumb, fr_partit_cw' tmpa = 0.0 tmpb = 0.0 tmpc = 0.0 do itype = 1, ntype_aer do isize = 1, nsize_aer(itype) if ( .not. do_cloudchem_aer(isize,itype) ) cycle write(lunxx,9820) xmass(isize,itype), fmass(isize,itype), & xnumbsv(isize,itype), xnumb(isize,itype), & fnumb(isize,itype), fr_partit_cw(isize,itype) tmpa = tmpa + fmass(isize,itype) tmpb = tmpb + fnumb(isize,itype) tmpc = tmpc + 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) (tmpa-1.0), (tmpb-1.0), (tmpc-1.0) if (icase .le. -5) then write(*,*) '*** stopping in partition_cldwtr at icase =', icase call wrf_error_fatal('*** stopping in partition_cldwtr') end if 9800 format( a ) 9810 format( 5i10 ) 9820 format( 6(1pe10.2) ) end if #endif return end subroutine sorgam_partition_cldwtr !----------------------------------------------------------------------- subroutine sorgam_distribute_bulk_changes( & rbox, rbox_sv1, fr_partit_cw, & rbulk_cwaer, lptr_yyy_cwaer, & id, it, jt, kt, icase ) use module_state_description, only: & param_first_scalar, num_chem use module_scalar_tables, only: chem_dname_table use module_data_sorgam, only: & maxd_asize, maxd_atype, & cw_phase, nsize_aer, ntype_aer, do_cloudchem_aer, & lptr_so4_aer, lptr_no3_aer, lptr_nh4_aer, & lptr_orgpa_aer, lptr_ec_aer, lptr_p25_aer implicit none ! subr arguments integer, intent(in) :: id, it, jt, kt, icase integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy) real, intent(inout), dimension( 1:num_chem ) :: 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, lunxxaa, lunxxbb, lyyy integer :: p1st #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) integer :: lunxx #endif 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 ( .not. do_cloudchem_aer(isize,itype) ) cycle 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) if ( .not. do_cloudchem_aer(isize,itype) ) cycle 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) if ( .not. do_cloudchem_aer(isize,itype) ) cycle 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) if ( .not. do_cloudchem_aer(isize,itype) ) cycle 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 ) || defined ( cctemp_testing_active ) ! lunxxaa = 86 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxxaa = 171 #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 ) || defined ( cctemp_testing_active ) ! lunxxbb = 86 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxxbb = 171 #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 ) || defined ( cctemp_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) l = lptr_yyy_cwaer(1,1,lyyy) if (l .ge. p1st) then write(lunxx,9801) 'distribute_bulk_changes - ', & chem_dname_table(id,l)(1:12), ' - icase, lyyy, l11' else write(lunxx,9801) 'distribute_bulk_changes - ', & 'name = ?????', ' - icase, lyyy, l11' end if write(lunxx,9810) icase, lyyy, l 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) if ( .not. do_cloudchem_aer(isize,itype) ) cycle l = lptr_yyy_cwaer(isize,itype,lyyy) if (l .lt. p1st) cycle 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 .le. -5) then write(*,*) & '*** stop in distribute_bulk_changes diags, icase =', icase call wrf_error_fatal('*** stop in distribute_bulk_changes diags') end if end if #endif return end subroutine sorgam_distribute_bulk_changes !----------------------------------------------------------------------- subroutine sorgam_cloudchem_apply_mode_transfer( & rbox, rbox_sv1, xvol_old, & id, it, jt, kt, icase ) use module_state_description, only: & param_first_scalar, num_chem use module_scalar_tables, only: chem_dname_table use module_data_sorgam, only: & pirs, & msectional, & maxd_asize, maxd_atype, & ai_phase, cw_phase, nsize_aer, ntype_aer, ncomp_aer, & do_cloudchem_aer, massptr_aer, numptr_aer, dens_aer, & sigmag_aer, dcen_sect, dlo_sect, dhi_sect, & volumcen_sect, volumlo_sect, volumhi_sect, & lptr_so4_aer, lptr_nh4_aer, lptr_p25_aer use module_aerosols_sorgam, only: getaf implicit none ! subr arguments integer, intent(in) :: id, it, jt, kt, icase real, intent(inout), dimension( 1:num_chem ) :: rbox, rbox_sv1 real, intent(in), dimension( 2, 3 ) :: xvol_old ! local variables integer :: idum_msect integer :: ii, isize, isize_ait, isize_acc, itype integer :: jj integer :: l, lfrm, ltoo, ll integer :: lptr_dum(maxd_asize,maxd_atype) integer :: p1st #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) integer :: lunxx #endif logical :: skip_xfer real :: delvol(2) real :: fracrem_num, fracrem_vol, fracxfr_num, fracxfr_vol real :: rbox_sv2(1:num_chem) real :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf real :: tmp_cwnumfrac, tmp_cwvolfrac real :: tmp_dpmeanvol, tmp_lnsg real :: xcut_num, xcut_vol real :: xdgnum_aaa(2), xlnsg(2) real :: xnum_aaa(2,3) real :: xvol_aaa(2,3) #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) ! diagnostics for testing lunxx = -1 if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171 #endif p1st = param_first_scalar ! ! initial calculations for aitken (ii=1) and accum (ii=2) modes ! skip_xfer = .false. do ii = 1, 2 itype = 1 isize = ii ! calculate new volumes for activated (jj=1), interstitial (jj=2), and combined (jj=3) tmpa = 0.0 do ll = 1, ncomp_aer(itype) l = massptr_aer(ll,isize,itype,cw_phase) if (l >= p1st) tmpa = tmpa + rbox(l)/dens_aer(ll,itype) end do xvol_aaa(ii,1) = tmpa*1.0e-6 xvol_aaa(ii,2) = xvol_old(ii,2) xnum_aaa(ii,1) = rbox(numptr_aer(isize,itype,cw_phase)) xnum_aaa(ii,2) = rbox(numptr_aer(isize,itype,ai_phase)) xvol_aaa(ii,3) = xvol_aaa(ii,1) + xvol_aaa(ii,2) xnum_aaa(ii,3) = xnum_aaa(ii,1) + xnum_aaa(ii,2) delvol(ii) = xvol_aaa(ii,1) - xvol_old(ii,1) ! check for negligible number or volume in aitken mode if (ii == 1) then if (xvol_aaa(ii,3) < smallvolaa) then skip_xfer = .true. exit end if end if ! calculate dgnum for activated+interstitial combined using new volume if (xvol_aaa(ii,3) < smallvolaa) then tmp_dpmeanvol = dcen_sect(isize,itype) else tmp_dpmeanvol = xvol_aaa(ii,3)/xnum_aaa(ii,3) tmp_dpmeanvol = (tmp_dpmeanvol*6.0/pirs)**0.33333333 end if xlnsg(ii) = log(sigmag_aer(isize,itype)) xdgnum_aaa(ii) = tmp_dpmeanvol*exp(-1.5*xlnsg(ii)*xlnsg(ii)) ! tmp_cwvolfrac = (cw volume)/(ai + cw volume) tmp_cwvolfrac = xvol_aaa(ii,1)/max(xvol_aaa(ii,3),1.e-30) #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) ! diagnostics for testing if (lunxx > 0) then if (ii == 1) write(lunxx,'(a)') write(lunxx,'(a,3i4,i8,2x,l2)') 'merge i,j,k, ii,skip', & it, jt, kt, ii, skip_xfer write(lunxx,'(a,1p,5e12.4)') 'cw vol old/aaa, num aaa ', & xvol_old(ii,1), xvol_aaa(ii,1), xnum_aaa(ii,1) write(lunxx,'(a,1p,5e12.4)') 'ai vol old/aaa, num aaa ', & xvol_old(ii,2), xvol_aaa(ii,2), xnum_aaa(ii,2) write(lunxx,'(a,1p,5e12.4)') 'a+c vol old/aaa, num aaa ', & xvol_old(ii,3), xvol_aaa(ii,3), xnum_aaa(ii,3) write(lunxx,'(a,1p,5e12.4)') 'cwnum/volfrac, dpmeanvol, dgnum ', & (xnum_aaa(ii,1)/max(xnum_aaa(ii,3),1.e-30)), & tmp_cwvolfrac, tmp_dpmeanvol, xdgnum_aaa(ii) end if #endif end do ! ii = 1, 2 ! check for mode merging if ( skip_xfer ) return if (delvol(1) > delvol(2)) then continue else if ( (xdgnum_aaa(1) > 0.03e-4) .and. (xnum_aaa(1,3) > xnum_aaa(2,3)) ) then continue else return end if #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) if (lunxx > 0) then if (delvol(1) > delvol(2)) then write(lunxx,'(a)') 'do merging - criterion 1' else if ( (xdgnum_aaa(1) > 0.03e-4) .and. (xnum_aaa(1,3) > xnum_aaa(2,3)) ) then write(lunxx,'(a)') 'do merging - criterion 2' end if end if #endif ! ! calc transfer fractions for volume/mass and number ! approach follows that in module_aerosols_sorgam (subr aerostep) except ! >> the first steps of the calculation are done using the total ! (interstitial+activated) size distributions ! >> the number and volume (moment-3) transfer amounts are then limited ! to the number/volume of the aitken activated distribution ! ! xcut_num = [ln(dintsect/dgnuc)/xxlsgn], where dintsect is the diameter ! at which the aitken-mode and accum-mode number distribs intersect (overlap). tmpa = sqrt(2.0) ! aaa = getaf( nu0, ac0, dgnuc, dgacc, xxlsgn, xxlsga, sqrt2 ) xcut_num = tmpa * getaf( xnum_aaa(1,3), xnum_aaa(2,3), & xdgnum_aaa(1), xdgnum_aaa(2), xlnsg(1), xlnsg(2), tmpa ) ! forcing xcut_vol>0 means that no more than half of the aitken volume ! will be transferred tmpd = xcut_num tmpc = 3.0*xlnsg(1) xcut_vol = max( xcut_num-tmpc, 0.0 ) xcut_num = xcut_vol + tmpc fracxfr_vol = norm01_uptail( xcut_vol ) fracxfr_num = norm01_uptail( xcut_num ) tmpe = fracxfr_vol ; tmpf = fracxfr_num tmp_cwvolfrac = xvol_aaa(1,1)/max(xvol_aaa(1,3),1.e-30) tmp_cwnumfrac = xnum_aaa(1,1)/max(xnum_aaa(1,3),1.e-30) if ( (fracxfr_vol >= tmp_cwvolfrac) .or. & (fracxfr_num >= tmp_cwnumfrac) ) then ! limit volume fraction transferred to tmp_cwvolfrac ! limit number fraction transferred to tmp_cwnumfrac fracxfr_num = 1.0 fracxfr_vol = 1.0 else ! at this point, fracxfr_num/vol are fraction of ! interstitial+activated num/vol to be transferred ! convert them to fraction of activated num/vol to be transferred fracxfr_vol = fracxfr_vol/max(1.0e-10,tmp_cwvolfrac) fracxfr_num = fracxfr_num/max(1.0e-10,tmp_cwnumfrac) ! number fraction transferred cannot exceed volume fraction fracxfr_num = min( fracxfr_num, fracxfr_vol ) end if fracrem_vol = 1.0 - fracxfr_vol fracrem_num = 1.0 - fracxfr_num #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) if (lunxx > 0) then write(lunxx,'(a,1p,5e12.4)') 'xcut_num1/num2/vol ', & tmpd, xcut_num, xcut_vol write(lunxx,'(a,1p,5e12.4)') 'fracxfr_num3/vol3/num1,vol1 ', & tmpf, tmpe, fracxfr_num, fracxfr_vol end if #endif if ( skip_xfer ) return ! do the transfer rbox_sv2(:) = rbox(:) itype = 1 isize_ait = 1 isize_acc = 2 lfrm = numptr_aer(isize_ait,itype,cw_phase) ltoo = numptr_aer(isize_acc,itype,cw_phase) rbox(ltoo) = rbox(ltoo) + rbox(lfrm)*fracxfr_num rbox(lfrm) = rbox(lfrm)*fracrem_num do ll = 1, ncomp_aer(itype) lfrm = massptr_aer(ll,isize_ait,itype,cw_phase) ltoo = massptr_aer(ll,isize_acc,itype,cw_phase) if (lfrm >= p1st) then if (ltoo >= p1st) rbox(ltoo) = rbox(ltoo) & + rbox(lfrm)*fracxfr_vol rbox(lfrm) = rbox(lfrm)*fracrem_vol end if end do #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active ) ! more diagnostics for testing if (lunxx .gt. 0) then do ll = 1, 4 if (ll .eq. 1) then lptr_dum(:,:) = lptr_so4_aer(:,:,cw_phase) else if (ll .eq. 2) then lptr_dum(:,:) = lptr_nh4_aer(:,:,cw_phase) else if (ll .eq. 3) then lptr_dum(:,:) = lptr_p25_aer(:,:,cw_phase) else if (ll .eq. 4) then lptr_dum(:,:) = numptr_aer(:,:,cw_phase) end if if (ll .eq. 1) write(lunxx,'(a)') write(lunxx,'(2a,i6,i3,2x,a)') 'sorgam_cloudchem_apply_mode_transfer', & ' - icase, ll', icase, ll, & chem_dname_table(id,lptr_dum(1,1))(1:12) write(lunxx,'(a)') ' ty sz rbox_sv1, rbox, rsub' itype = 1 do ii = 1, 2 if (ii == 1) then isize = isize_ait else isize = isize_acc end if l = lptr_dum(isize,itype) write(lunxx,'(2i3,1p,5e14.6)') & itype, isize, rbox_sv1(l), rbox_sv2(l), rbox(l) end do end do if (icase .le. -5) then write(*,*) & '*** stop in sorgam_cloudchem_apply_mode_transfer diags, icase =', & icase call wrf_error_fatal('*** stop in sorgam_cloudchem_apply_mode_transfer diags') end if end if #endif return end subroutine sorgam_cloudchem_apply_mode_transfer !----------------------------------------------------------------------- real function norm01_uptail( x ) ! ! norm01_uptail = cummulative pdf complement of normal(0,1) pdf ! = integral from x to +infinity of [normal(0,1) pdf] ! ! erfc_num_recipes is from press et al, numerical recipes, 1990, page 164 ! implicit none real x, xabs real*8 erfc_approx, tmpa, t, z xabs = abs(x) if (xabs >= 12.962359) then if (x > 0.0) then norm01_uptail = 0.0 else norm01_uptail = 1.0 end if return end if z = xabs / sqrt(2.0_8) t = 1.0_8/(1.0_8 + 0.5_8*z) ! erfc_approx = ! & t*exp( -z*z - 1.26551223 + t*(1.00002368 + t*(0.37409196 + ! & t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 + ! & t*(-1.13520398 + ! & t*(1.48851587 + t*(-0.82215223 + t*0.17087277 ))))))))) tmpa = ( -z*z - 1.26551223_8 + t*(1.00002368_8 + t*(0.37409196_8 + & t*(0.09678418_8 + t*(-0.18628806_8 + t*(0.27886807_8 + & t*(-1.13520398_8 + & t*(1.48851587_8 + t*(-0.82215223_8 + t*0.17087277_8 ))))))))) erfc_approx = t * exp(tmpa) if (x .lt. 0.0) erfc_approx = 2.0_8 - erfc_approx norm01_uptail = 0.5_8 * erfc_approx return end function norm01_uptail !----------------------------------------------------------------------- real function norm01_uptail_inv( x ) ! ! norm01_uptail_inv = inverse of norm01_uptail ! if y = norm01_uptail_inv( x ), then ! {integral from y to +infinity of [normal(0,1) pdf]} = x ! y is computed using newton's method ! implicit none ! fn parameters real x ! local variables integer niter real dfdyinv, f, pi, sqrt2pi, tmpa, y, ynew parameter (pi = 3.1415926535897932384626434) if (x .le. 1.0e-38) then norm01_uptail_inv = 12.962359 return else if (x .ge. 1.0) then norm01_uptail_inv = -12.962359 return end if sqrt2pi = sqrt( 2.0*pi ) ! initial guess ! crude ! y = 3.0*(0.5 - x) ! better tmpa = x tmpa = max( 0.0, min( 1.0, tmpa ) ) tmpa = 4.0*tmpa*(1.0 - tmpa) tmpa = max( 1.0e-38, min( 1.0, tmpa ) ) y = sqrt( -(pi/2.0)*log(tmpa) ) if (x > 0.5) y = -y f = norm01_uptail(y) - x do niter = 1, 100 ! iterate - dfdy is computed analytically dfdyinv = -sqrt2pi * exp( 0.5*y*y ) ynew = y - f*dfdyinv f = norm01_uptail(ynew) - x if ( (ynew == y) .or. & (abs(f) <= abs(x)*1.0e-5) ) then exit end if y = ynew end do 9100 format( 'niter/x/f/y/ynew', i5, 4(1pe16.8) ) norm01_uptail_inv = ynew return end function norm01_uptail_inv !----------------------------------------------------------------------- subroutine sorgam_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_scalar_tables, only: chem_dname_table use module_configure, only: grid_config_rec_type use module_data_sorgam 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( & 'sorgam_cloudchem_dumpaa - ktau, i, j, k =', 4i5 ) do 2900 itype = 1, ntype_aer do 2900 isize = 1, nsize_aer(itype) if ( .not. do_cloudchem_aer(isize,itype) ) goto 2900 write(*,9110) isize 9110 format( / 'isize, itype =', 2i3 / & ' 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), chem_dname_table(id,l)(1:12), 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 sorgam_cloudchem_dumpaa !----------------------------------------------------------------------- end module module_sorgam_cloudchem