module mo_gas_phase_chemdr use shr_kind_mod, only : r8 => shr_kind_r8 use shr_const_mod, only : pi => shr_const_pi use constituents, only : pcnst use cam_history, only : fieldname_len use chem_mods, only : phtcnt, rxntot, gas_pcnst use chem_mods, only : rxt_tag_cnt, rxt_tag_lst, rxt_tag_map, extcnt use dust_model, only : dust_names, ndust => dust_nbin use ppgrid, only : pcols, pver use phys_control, only : phys_getopts use carma_flags_mod, only : carma_hetchem_feedback use chem_prod_loss_diags, only: chem_prod_loss_diags_init, chem_prod_loss_diags_out implicit none save private public :: gas_phase_chemdr, gas_phase_chemdr_inti public :: map2chm integer :: map2chm(pcnst) = 0 ! index map to/from chemistry/constituents list integer :: synoz_ndx, so4_ndx, h2o_ndx, o2_ndx, o_ndx, hno3_ndx, hcl_ndx, dst_ndx, cldice_ndx, snow_ndx integer :: o3_ndx, o3s_ndx integer :: het1_ndx integer :: ndx_cldfr, ndx_cmfdqr, ndx_nevapr, ndx_cldtop, ndx_prain integer :: ndx_h2so4 ! ! CCMI ! integer :: st80_25_ndx integer :: st80_25_tau_ndx integer :: aoa_nh_ndx integer :: aoa_nh_ext_ndx integer :: nh_5_ndx integer :: nh_50_ndx integer :: nh_50w_ndx integer :: sad_pbf_ndx integer :: cb1_ndx,cb2_ndx,oc1_ndx,oc2_ndx,dst1_ndx,dst2_ndx,sslt1_ndx,sslt2_ndx integer :: soa_ndx,soai_ndx,soam_ndx,soat_ndx,soab_ndx,soax_ndx character(len=fieldname_len),dimension(rxntot-phtcnt) :: rxn_names character(len=fieldname_len),dimension(phtcnt) :: pht_names character(len=fieldname_len),dimension(rxt_tag_cnt) :: tag_names character(len=fieldname_len),dimension(extcnt) :: extfrc_name logical :: pm25_srf_diag logical :: pm25_srf_diag_soa logical :: convproc_do_aer integer :: ele_temp_ndx, ion_temp_ndx contains subroutine gas_phase_chemdr_inti() use mo_chem_utls, only : get_spc_ndx, get_extfrc_ndx, get_rxt_ndx use cam_history, only : addfld,add_default,horiz_only use mo_chm_diags, only : chm_diags_inti use constituents, only : cnst_get_ind use physics_buffer, only : pbuf_get_index use rate_diags, only : rate_diags_init implicit none character(len=3) :: string integer :: n, m, err logical :: history_cesm_forcing !----------------------------------------------------------------------- call phys_getopts( convproc_do_aer_out = convproc_do_aer, history_cesm_forcing_out=history_cesm_forcing ) ndx_h2so4 = get_spc_ndx('H2SO4') ! ! CCMI ! st80_25_ndx = get_spc_ndx ('ST80_25') st80_25_tau_ndx = get_rxt_ndx ('ST80_25_tau') aoa_nh_ndx = get_spc_ndx ('AOA_NH') aoa_nh_ext_ndx = get_extfrc_ndx('AOA_NH') nh_5_ndx = get_spc_ndx('NH_5') nh_50_ndx = get_spc_ndx('NH_50') nh_50w_ndx = get_spc_ndx('NH_50W') ! cb1_ndx = get_spc_ndx('CB1') cb2_ndx = get_spc_ndx('CB2') oc1_ndx = get_spc_ndx('OC1') oc2_ndx = get_spc_ndx('OC2') dst1_ndx = get_spc_ndx('DST01') dst2_ndx = get_spc_ndx('DST02') sslt1_ndx = get_spc_ndx('SSLT01') sslt2_ndx = get_spc_ndx('SSLT02') soa_ndx = get_spc_ndx('SOA') soam_ndx = get_spc_ndx('SOAM') soai_ndx = get_spc_ndx('SOAI') soat_ndx = get_spc_ndx('SOAT') soab_ndx = get_spc_ndx('SOAB') soax_ndx = get_spc_ndx('SOAX') pm25_srf_diag = cb1_ndx>0 .and. cb2_ndx>0 .and. oc1_ndx>0 .and. oc2_ndx>0 & .and. dst1_ndx>0 .and. dst2_ndx>0 .and. sslt1_ndx>0 .and. sslt2_ndx>0 & .and. soa_ndx>0 pm25_srf_diag_soa = cb1_ndx>0 .and. cb2_ndx>0 .and. oc1_ndx>0 .and. oc2_ndx>0 & .and. dst1_ndx>0 .and. dst2_ndx>0 .and. sslt1_ndx>0 .and. sslt2_ndx>0 & .and. soam_ndx>0 .and. soai_ndx>0 .and. soat_ndx>0 .and. soab_ndx>0 .and. soax_ndx>0 if ( pm25_srf_diag .or. pm25_srf_diag_soa) then call addfld('PM25_SRF',horiz_only,'I','kg/kg','bottom layer PM2.5 mixing ratio' ) endif call addfld('U_SRF',horiz_only,'I','m/s','bottom layer wind velocity' ) call addfld('V_SRF',horiz_only,'I','m/s','bottom layer wind velocity' ) call addfld('Q_SRF',horiz_only,'I','kg/kg','bottom layer specific humidity' ) ! het1_ndx= get_rxt_ndx('het1') o3_ndx = get_spc_ndx('O3') o3s_ndx = get_spc_ndx('O3S') o_ndx = get_spc_ndx('O') o2_ndx = get_spc_ndx('O2') so4_ndx = get_spc_ndx('SO4') h2o_ndx = get_spc_ndx('H2O') hno3_ndx = get_spc_ndx('HNO3') hcl_ndx = get_spc_ndx('HCL') dst_ndx = get_spc_ndx( dust_names(1) ) synoz_ndx = get_extfrc_ndx( 'SYNOZ' ) call cnst_get_ind( 'CLDICE', cldice_ndx ) call cnst_get_ind( 'SNOWQM', snow_ndx, abort=.false. ) do m = 1,extcnt WRITE(UNIT=string, FMT='(I2.2)') m extfrc_name(m) = 'extfrc_'// trim(string) call addfld( extfrc_name(m), (/ 'lev' /), 'I', ' ', 'ext frcing' ) end do do n = 1,rxt_tag_cnt tag_names(n) = trim(rxt_tag_lst(n)) if (n<=phtcnt) then call addfld( tag_names(n), (/ 'lev' /), 'I', '/s', 'photolysis rate' ) else call addfld( tag_names(n), (/ 'lev' /), 'I', '/cm3/s', 'reaction rate' ) endif call add_default(tag_names(n), 2,' ') enddo do n = 1,phtcnt WRITE(UNIT=string, FMT='(I3.3)') n pht_names(n) = 'J_' // trim(string) call addfld( pht_names(n), (/ 'lev' /), 'I', '/s', 'photolysis rate' ) enddo do n = 1,rxntot-phtcnt WRITE(UNIT=string, FMT='(I3.3)') n rxn_names(n) = 'R_' // trim(string) call addfld( rxn_names(n), (/ 'lev' /), 'I', '/cm3/s', 'reaction rate' ) enddo call addfld( 'DTCBS', horiz_only, 'I', ' ','photolysis diagnostic black carbon OD' ) call addfld( 'DTOCS', horiz_only, 'I', ' ','photolysis diagnostic organic carbon OD' ) call addfld( 'DTSO4', horiz_only, 'I', ' ','photolysis diagnostic SO4 OD' ) call addfld( 'DTSOA', horiz_only, 'I', ' ','photolysis diagnostic SOA OD' ) call addfld( 'DTANT', horiz_only, 'I', ' ','photolysis diagnostic NH4SO4 OD' ) call addfld( 'DTSAL', horiz_only, 'I', ' ','photolysis diagnostic salt OD' ) call addfld( 'DTDUST', horiz_only, 'I', ' ','photolysis diagnostic dust OD' ) call addfld( 'DTTOTAL', horiz_only, 'I', ' ','photolysis diagnostic total aerosol OD' ) call addfld( 'FRACDAY', horiz_only, 'I', ' ','photolysis diagnostic fraction of day' ) call addfld( 'QDSAD', (/ 'lev' /), 'I', '/s', 'water vapor sad delta' ) call addfld( 'SAD_STRAT', (/ 'lev' /), 'I', 'cm2/cm3', 'stratospheric aerosol SAD' ) call addfld( 'SAD_SULFC', (/ 'lev' /), 'I', 'cm2/cm3', 'chemical sulfate aerosol SAD' ) call addfld( 'SAD_SAGE', (/ 'lev' /), 'I', 'cm2/cm3', 'SAGE sulfate aerosol SAD' ) call addfld( 'SAD_LNAT', (/ 'lev' /), 'I', 'cm2/cm3', 'large-mode NAT aerosol SAD' ) call addfld( 'SAD_ICE', (/ 'lev' /), 'I', 'cm2/cm3', 'water-ice aerosol SAD' ) call addfld( 'RAD_SULFC', (/ 'lev' /), 'I', 'cm', 'chemical sad sulfate' ) call addfld( 'RAD_LNAT', (/ 'lev' /), 'I', 'cm', 'large nat radius' ) call addfld( 'RAD_ICE', (/ 'lev' /), 'I', 'cm', 'sad ice' ) call addfld( 'SAD_TROP', (/ 'lev' /), 'I', 'cm2/cm3', 'tropospheric aerosol SAD' ) call addfld( 'SAD_AERO', (/ 'lev' /), 'I', 'cm2/cm3', 'aerosol surface area density' ) if (history_cesm_forcing) then call add_default ('SAD_AERO',8,' ') endif call addfld( 'REFF_AERO', (/ 'lev' /), 'I', 'cm', 'aerosol effective radius' ) call addfld( 'SULF_TROP', (/ 'lev' /), 'I', 'mol/mol', 'tropospheric aerosol SAD' ) call addfld( 'QDSETT', (/ 'lev' /), 'I', '/s', 'water vapor settling delta' ) call addfld( 'QDCHEM', (/ 'lev' /), 'I', '/s', 'water vapor chemistry delta') call addfld( 'HNO3_TOTAL', (/ 'lev' /), 'I', 'mol/mol', 'total HNO3' ) call addfld( 'HNO3_STS', (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HNO3' ) call addfld( 'HNO3_NAT', (/ 'lev' /), 'I', 'mol/mol', 'NAT condensed HNO3' ) call addfld( 'HNO3_GAS', (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hno3' ) call addfld( 'H2O_GAS', (/ 'lev' /), 'I', 'mol/mol', 'gas-phase h2o' ) call addfld( 'HCL_TOTAL', (/ 'lev' /), 'I', 'mol/mol', 'total hcl' ) call addfld( 'HCL_GAS', (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hcl' ) call addfld( 'HCL_STS', (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HCL' ) if (het1_ndx>0) then call addfld( 'het1_total', (/ 'lev' /), 'I', '/s', 'total N2O5 + H2O het rate constant' ) endif call addfld( 'SZA', horiz_only, 'I', 'degrees', 'solar zenith angle' ) call chm_diags_inti() call rate_diags_init() !----------------------------------------------------------------------- ! get pbuf indicies !----------------------------------------------------------------------- ndx_cldfr = pbuf_get_index('CLD') ndx_cmfdqr = pbuf_get_index('RPRDTOT') ndx_nevapr = pbuf_get_index('NEVAPR') ndx_prain = pbuf_get_index('PRAIN') ndx_cldtop = pbuf_get_index('CLDTOP') sad_pbf_ndx= pbuf_get_index('VOLC_SAD',errcode=err) ! prescribed strat aerosols (volcanic) if (.not.sad_pbf_ndx>0) sad_pbf_ndx = pbuf_get_index('SADSULF',errcode=err) ! CARMA's version of strat aerosols ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index ! diagnostics for stratospheric heterogeneous reactions call addfld( 'GAMMA_HET1', (/ 'lev' /), 'I', '1', 'Reaction Probability' ) call addfld( 'GAMMA_HET2', (/ 'lev' /), 'I', '1', 'Reaction Probability' ) call addfld( 'GAMMA_HET3', (/ 'lev' /), 'I', '1', 'Reaction Probability' ) call addfld( 'GAMMA_HET4', (/ 'lev' /), 'I', '1', 'Reaction Probability' ) call addfld( 'GAMMA_HET5', (/ 'lev' /), 'I', '1', 'Reaction Probability' ) call addfld( 'GAMMA_HET6', (/ 'lev' /), 'I', '1', 'Reaction Probability' ) call addfld( 'WTPER', (/ 'lev' /), 'I', '%', 'H2SO4 Weight Percent' ) call chem_prod_loss_diags_init end subroutine gas_phase_chemdr_inti !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, & phis, zm, zi, calday, & tfld, pmid, pdel, pint, & cldw, troplev, troplevchem, & ncldwtr, ufld, vfld, & delt, ps, xactive_prates, & fsds, ts, asdir, ocnfrac, icefrac, & precc, precl, snowhland, ghg_chem, latmapback, & drydepflx, wetdepflx, cflx, fire_sflx, fire_ztop, nhx_nitrogen_flx, noy_nitrogen_flx, qtend, pbuf) !----------------------------------------------------------------------- ! ... Chem_solver advances the volumetric mixing ratio ! forward one time step via a combination of explicit, ! ebi, hov, fully implicit, and/or rodas algorithms. !----------------------------------------------------------------------- use chem_mods, only : nabscol, nfs, indexm, clscnt4 use physconst, only : rga use mo_photo, only : set_ub_col, setcol, table_photo, xactive_photo use mo_exp_sol, only : exp_sol use mo_imp_sol, only : imp_sol use mo_setrxt, only : setrxt use mo_adjrxt, only : adjrxt use mo_phtadj, only : phtadj use llnl_O1D_to_2OH_adj,only : O1D_to_2OH_adj use mo_usrrxt, only : usrrxt use mo_setinv, only : setinv use mo_negtrc, only : negtrc use mo_sulf, only : sulf_interp use mo_setext, only : setext use fire_emissions, only : fire_emissions_vrt use mo_sethet, only : sethet use mo_drydep, only : drydep, set_soilw use seq_drydep_mod, only : DD_XLND, DD_XATM, DD_TABL, drydep_method use mo_fstrat, only : set_fstrat_vals, set_fstrat_h2o use noy_ubc, only : noy_ubc_set use mo_flbc, only : flbc_set use phys_grid, only : get_rlat_all_p, get_rlon_all_p, get_lat_all_p, get_lon_all_p use mo_mean_mass, only : set_mean_mass use cam_history, only : outfld use wv_saturation, only : qsat use constituents, only : cnst_mw use mo_drydep, only : has_drydep use time_manager, only : get_ref_date use mo_ghg_chem, only : ghg_chem_set_rates, ghg_chem_set_flbc use mo_sad, only : sad_strat_calc use charge_neutrality, only : charge_balance use mo_strato_rates, only : ratecon_sfstrat use mo_aero_settling, only : strat_aer_settling use shr_orb_mod, only : shr_orb_decl use cam_control_mod, only : lambm0, eccen, mvelpp, obliqr use mo_strato_rates, only : has_strato_chem use short_lived_species,only: set_short_lived_species,get_short_lived_species use mo_chm_diags, only : chm_diags, het_diags use perf_mod, only : t_startf, t_stopf use gas_wetdep_opts, only : gas_wetdep_method use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx use infnan, only : nan, assignment(=) use rate_diags, only : rate_diags_calc use mo_mass_xforms, only : mmr2vmr, vmr2mmr, h2o_to_vmr, mmr2vmri use orbit, only : zenith ! ! LINOZ ! use lin_strat_chem, only : do_lin_strat_chem, lin_strat_chem_solve use linoz_data, only : has_linoz_data ! ! for aqueous chemistry and aerosol growth ! use aero_model, only : aero_model_gasaerexch use aero_model, only : aero_model_strat_surfarea use time_manager, only : is_first_step implicit none !----------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: lchnk ! chunk index integer, intent(in) :: ncol ! number columns in chunk integer, intent(in) :: imozart ! gas phase start index in q real(r8), intent(in) :: delt ! timestep (s) real(r8), intent(in) :: calday ! day of year real(r8), intent(in) :: ps(pcols) ! surface pressure real(r8), intent(in) :: phis(pcols) ! surface geopotential real(r8),target,intent(in) :: tfld(pcols,pver) ! midpoint temperature (K) real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures (Pa) real(r8), intent(in) :: pdel(pcols,pver) ! pressure delta about midpoints (Pa) real(r8), intent(in) :: ufld(pcols,pver) ! zonal velocity (m/s) real(r8), intent(in) :: vfld(pcols,pver) ! meridional velocity (m/s) real(r8), intent(in) :: cldw(pcols,pver) ! cloud water (kg/kg) real(r8), intent(in) :: ncldwtr(pcols,pver) ! droplet number concentration (#/kg) real(r8), intent(in) :: zm(pcols,pver) ! midpoint geopotential height above the surface (m) real(r8), intent(in) :: zi(pcols,pver+1) ! interface geopotential height above the surface (m) real(r8), intent(in) :: pint(pcols,pver+1) ! interface pressures (Pa) real(r8), intent(in) :: q(pcols,pver,pcnst) ! species concentrations (kg/kg) real(r8),pointer, intent(in) :: fire_sflx(:,:) ! fire emssions surface flux (kg/m^2/s) real(r8),pointer, intent(in) :: fire_ztop(:) ! top of vertical distribution of fire emssions (m) logical, intent(in) :: xactive_prates real(r8), intent(in) :: fsds(pcols) ! longwave down at sfc real(r8), intent(in) :: icefrac(pcols) ! sea-ice areal fraction real(r8), intent(in) :: ocnfrac(pcols) ! ocean areal fraction real(r8), intent(in) :: asdir(pcols) ! albedo: shortwave, direct real(r8), intent(in) :: ts(pcols) ! sfc temp (merged w/ocean if coupled) real(r8), intent(in) :: precc(pcols) ! real(r8), intent(in) :: precl(pcols) ! real(r8), intent(in) :: snowhland(pcols) ! logical, intent(in) :: ghg_chem integer, intent(in) :: latmapback(pcols) integer, intent(in) :: troplev(pcols) ! trop/strat separation vertical index integer, intent(in) :: troplevchem(pcols) ! trop/strat chemistry separation vertical index real(r8), intent(inout) :: qtend(pcols,pver,pcnst) ! species tendencies (kg/kg/s) real(r8), intent(inout) :: cflx(pcols,pcnst) ! constituent surface flux (kg/m^2/s) real(r8), intent(out) :: drydepflx(pcols,pcnst) ! dry deposition flux (kg/m^2/s) real(r8), intent(in) :: wetdepflx(pcols,pcnst) ! wet deposition flux (kg/m^2/s) real(r8), intent(out) :: nhx_nitrogen_flx(pcols) real(r8), intent(out) :: noy_nitrogen_flx(pcols) type(physics_buffer_desc), pointer :: pbuf(:) !----------------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------------- real(r8), parameter :: m2km = 1.e-3_r8 real(r8), parameter :: Pa2mb = 1.e-2_r8 real(r8), pointer :: prain(:,:) real(r8), pointer :: nevapr(:,:) real(r8), pointer :: cmfdqr(:,:) real(r8), pointer :: cldfr(:,:) real(r8), pointer :: cldtop(:) integer :: i, k, m, n integer :: tim_ndx real(r8) :: delt_inverse real(r8) :: esfact integer :: latndx(pcols) ! chunk lat indicies integer :: lonndx(pcols) ! chunk lon indicies real(r8) :: invariants(ncol,pver,nfs) real(r8) :: col_dens(ncol,pver,max(1,nabscol)) ! column densities (molecules/cm^2) real(r8) :: col_delta(ncol,0:pver,max(1,nabscol)) ! layer column densities (molecules/cm^2) real(r8) :: extfrc(ncol,pver,max(1,extcnt)) real(r8) :: vmr(ncol,pver,gas_pcnst) ! xported species (vmr) real(r8) :: reaction_rates(ncol,pver,max(1,rxntot)) ! reaction rates real(r8) :: depvel(ncol,gas_pcnst) ! dry deposition velocity (cm/s) real(r8) :: het_rates(ncol,pver,max(1,gas_pcnst)) ! washout rate (1/s) real(r8), dimension(ncol,pver) :: & h2ovmr, & ! water vapor volume mixing ratio mbar, & ! mean wet atmospheric mass ( amu ) zmid, & ! midpoint geopotential in km zmidr, & ! midpoint geopotential in km realitive to surf sulfate, & ! trop sulfate aerosols pmb ! pressure at midpoints ( hPa ) real(r8), dimension(ncol,pver) :: & cwat, & ! cloud water mass mixing ratio (kg/kg) wrk real(r8), dimension(ncol,pver+1) :: & zintr ! interface geopotential in km realitive to surf real(r8), dimension(ncol,pver+1) :: & zint ! interface geopotential in km real(r8), dimension(ncol) :: & zen_angle, & ! solar zenith angles zsurf, & ! surface height (m) rlats, rlons ! chunk latitudes and longitudes (radians) real(r8) :: sza(ncol) ! solar zenith angles (degrees) real(r8), parameter :: rad2deg = 180._r8/pi ! radians to degrees conversion factor real(r8) :: relhum(ncol,pver) ! relative humidity real(r8) :: satv(ncol,pver) ! wrk array for relative humidity real(r8) :: satq(ncol,pver) ! wrk array for relative humidity integer :: j integer :: ltrop_sol(pcols) ! tropopause vertical index used in chem solvers real(r8), pointer :: strato_sad(:,:) ! stratospheric sad (1/cm) real(r8) :: sad_trop(pcols,pver) ! total tropospheric sad (cm^2/cm^3) real(r8) :: reff(pcols,pver) ! aerosol effective radius (cm) real(r8) :: reff_strat(pcols,pver) ! stratospheric aerosol effective radius (cm) real(r8) :: tvs(pcols) integer :: ncdate,yr,mon,day,sec real(r8) :: wind_speed(pcols) ! surface wind speed (m/s) logical, parameter :: dyn_soilw = .false. logical :: table_soilw real(r8) :: soilw(pcols) real(r8) :: prect(pcols) real(r8) :: sflx(pcols,gas_pcnst) real(r8) :: wetdepflx_diag(pcols,gas_pcnst) real(r8) :: dust_vmr(ncol,pver,ndust) real(r8) :: dt_diag(pcols,8) ! od diagnostics real(r8) :: fracday(pcols) ! fraction of day real(r8) :: o2mmr(ncol,pver) ! o2 concentration (kg/kg) real(r8) :: ommr(ncol,pver) ! o concentration (kg/kg) real(r8) :: mmr(pcols,pver,gas_pcnst) ! chem working concentrations (kg/kg) real(r8) :: mmr_new(pcols,pver,gas_pcnst) ! chem working concentrations (kg/kg) real(r8) :: hno3_gas(ncol,pver) ! hno3 gas phase concentration (mol/mol) real(r8) :: hno3_cond(ncol,pver,2) ! hno3 condensed phase concentration (mol/mol) real(r8) :: hcl_gas(ncol,pver) ! hcl gas phase concentration (mol/mol) real(r8) :: hcl_cond(ncol,pver) ! hcl condensed phase concentration (mol/mol) real(r8) :: h2o_gas(ncol,pver) ! h2o gas phase concentration (mol/mol) real(r8) :: h2o_cond(ncol,pver) ! h2o condensed phase concentration (mol/mol) real(r8) :: cldice(pcols,pver) ! cloud water "ice" (kg/kg) real(r8) :: radius_strat(ncol,pver,3) ! radius of sulfate, nat, & ice ( cm ) real(r8) :: sad_strat(ncol,pver,3) ! surf area density of sulfate, nat, & ice ( cm^2/cm^3 ) real(r8) :: mmr_tend(pcols,pver,gas_pcnst) ! chemistry species tendencies (kg/kg/s) real(r8) :: qh2o(pcols,pver) ! specific humidity (kg/kg) real(r8) :: delta ! for aerosol formation.... real(r8) :: del_h2so4_gasprod(ncol,pver) real(r8) :: vmr0(ncol,pver,gas_pcnst) ! ! CCMI ! real(r8) :: xlat real(r8) :: pm25(ncol) real(r8) :: dlats(ncol) real(r8), dimension(ncol,pver) :: & ! aerosol reaction diagnostics gprob_n2o5, & gprob_cnt_hcl, & gprob_cnt_h2o, & gprob_bnt_h2o, & gprob_hocl_hcl, & gprob_hobr_hcl, & wtper real(r8), pointer :: ele_temp_fld(:,:) ! electron temperature pointer real(r8), pointer :: ion_temp_fld(:,:) ! ion temperature pointer real(r8) :: prod_out(ncol,pver,max(1,clscnt4)) real(r8) :: loss_out(ncol,pver,max(1,clscnt4)) if ( ele_temp_ndx>0 .and. ion_temp_ndx>0 .and. .not.is_first_step()) then call pbuf_get_field(pbuf, ele_temp_ndx, ele_temp_fld) call pbuf_get_field(pbuf, ion_temp_ndx, ion_temp_fld) else ele_temp_fld => tfld ion_temp_fld => tfld endif ! initialize to NaN to hopefully catch user defined rxts that go unset reaction_rates(:,:,:) = nan delt_inverse = 1._r8 / delt !----------------------------------------------------------------------- ! ... Get chunck latitudes and longitudes !----------------------------------------------------------------------- call get_lat_all_p( lchnk, ncol, latndx ) call get_lon_all_p( lchnk, ncol, lonndx ) call get_rlat_all_p( lchnk, ncol, rlats ) call get_rlon_all_p( lchnk, ncol, rlons ) tim_ndx = pbuf_old_tim_idx() call pbuf_get_field(pbuf, ndx_prain, prain, start=(/1,1/), kount=(/ncol,pver/)) call pbuf_get_field(pbuf, ndx_cldfr, cldfr, start=(/1,1,tim_ndx/), kount=(/ncol,pver,1/) ) call pbuf_get_field(pbuf, ndx_cmfdqr, cmfdqr, start=(/1,1/), kount=(/ncol,pver/)) call pbuf_get_field(pbuf, ndx_nevapr, nevapr, start=(/1,1/), kount=(/ncol,pver/)) call pbuf_get_field(pbuf, ndx_cldtop, cldtop ) reff_strat(:,:) = 0._r8 dlats(:) = rlats(:)*rad2deg ! convert to degrees !----------------------------------------------------------------------- ! ... Calculate cosine of zenith angle ! then cast back to angle (radians) !----------------------------------------------------------------------- call zenith( calday, rlats, rlons, zen_angle, ncol ) zen_angle(:) = acos( zen_angle(:) ) sza(:) = zen_angle(:) * rad2deg call outfld( 'SZA', sza, ncol, lchnk ) !----------------------------------------------------------------------- ! ... Xform geopotential height from m to km ! and pressure from Pa to mb !----------------------------------------------------------------------- zsurf(:ncol) = rga * phis(:ncol) do k = 1,pver zintr(:ncol,k) = m2km * zi(:ncol,k) zmidr(:ncol,k) = m2km * zm(:ncol,k) zmid(:ncol,k) = m2km * (zm(:ncol,k) + zsurf(:ncol)) zint(:ncol,k) = m2km * (zi(:ncol,k) + zsurf(:ncol)) pmb(:ncol,k) = Pa2mb * pmid(:ncol,k) end do zint(:ncol,pver+1) = m2km * (zi(:ncol,pver+1) + zsurf(:ncol)) zintr(:ncol,pver+1)= m2km * zi(:ncol,pver+1) !----------------------------------------------------------------------- ! ... map incoming concentrations to working array !----------------------------------------------------------------------- do m = 1,pcnst n = map2chm(m) if( n > 0 ) then mmr(:ncol,:,n) = q(:ncol,:,m) end if end do call get_short_lived_species( mmr, lchnk, ncol, pbuf ) !----------------------------------------------------------------------- ! ... Set atmosphere mean mass !----------------------------------------------------------------------- call set_mean_mass( ncol, mmr, mbar ) !----------------------------------------------------------------------- ! ... Xform from mmr to vmr !----------------------------------------------------------------------- call mmr2vmr( mmr(:ncol,:,:), vmr(:ncol,:,:), mbar(:ncol,:), ncol ) ! ! CCMI ! ! reset STE tracer to specific vmr of 200 ppbv ! if ( st80_25_ndx > 0 ) then where ( pmid(:ncol,:) < 80.e+2_r8 ) vmr(:ncol,:,st80_25_ndx) = 200.e-9_r8 end where end if ! ! reset AOA_NH, NH_5, NH_50, NH_50W surface mixing ratios between 30N and 50N ! if ( aoa_nh_ndx>0 .and. nh_5_ndx>0 .and. nh_50_ndx>0 .and. nh_50w_ndx>0 ) then do j=1,ncol xlat = dlats(j) if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then vmr(j,pver,nh_5_ndx) = 100.e-9_r8 vmr(j,pver,nh_50_ndx) = 100.e-9_r8 vmr(j,pver,nh_50w_ndx) = 100.e-9_r8 vmr(j,pver,aoa_nh_ndx) = 0._r8 end if end do end if if (h2o_ndx>0) then !----------------------------------------------------------------------- ! ... store water vapor in wrk variable !----------------------------------------------------------------------- qh2o(:ncol,:) = mmr(:ncol,:,h2o_ndx) h2ovmr(:ncol,:) = vmr(:ncol,:,h2o_ndx) else qh2o(:ncol,:) = q(:ncol,:,1) !----------------------------------------------------------------------- ! ... Xform water vapor from mmr to vmr and set upper bndy values !----------------------------------------------------------------------- call h2o_to_vmr( q(:ncol,:,1), h2ovmr(:ncol,:), mbar(:ncol,:), ncol ) call set_fstrat_h2o( h2ovmr, pmid, troplev, calday, ncol, lchnk ) endif !----------------------------------------------------------------------- ! ... force ion/electron balance !----------------------------------------------------------------------- call charge_balance( ncol, vmr ) !----------------------------------------------------------------------- ! ... Set the "invariants" !----------------------------------------------------------------------- call setinv( invariants, tfld, h2ovmr, vmr, pmid, ncol, lchnk, pbuf ) !----------------------------------------------------------------------- ! ... stratosphere aerosol surface area !----------------------------------------------------------------------- if (sad_pbf_ndx>0) then call pbuf_get_field(pbuf, sad_pbf_ndx, strato_sad) else allocate(strato_sad(pcols,pver)) strato_sad(:,:) = 0._r8 ! Prognostic modal stratospheric sulfate: compute dry strato_sad call aero_model_strat_surfarea( ncol, mmr, pmid, tfld, troplevchem, pbuf, strato_sad, reff_strat ) endif stratochem: if ( has_strato_chem ) then !----------------------------------------------------------------------- ! ... initialize condensed and gas phases; all hno3 to gas !----------------------------------------------------------------------- hcl_cond(:,:) = 0.0_r8 hcl_gas (:,:) = 0.0_r8 do k = 1,pver hno3_gas(:,k) = vmr(:,k,hno3_ndx) h2o_gas(:,k) = h2ovmr(:,k) hcl_gas(:,k) = vmr(:,k,hcl_ndx) wrk(:,k) = h2ovmr(:,k) if (snow_ndx>0) then cldice(:ncol,k) = q(:ncol,k,cldice_ndx) + q(:ncol,k,snow_ndx) else cldice(:ncol,k) = q(:ncol,k,cldice_ndx) endif end do do m = 1,2 do k = 1,pver hno3_cond(:,k,m) = 0._r8 end do end do call mmr2vmri( cldice(:ncol,:), h2o_cond(:ncol,:), mbar(:ncol,:), cnst_mw(cldice_ndx), ncol ) !----------------------------------------------------------------------- ! ... call SAD routine !----------------------------------------------------------------------- call sad_strat_calc( lchnk, invariants(:ncol,:,indexm), pmb, tfld, hno3_gas, & hno3_cond, h2o_gas, h2o_cond, hcl_gas, hcl_cond, strato_sad(:ncol,:), radius_strat, & sad_strat, ncol, pbuf ) ! NOTE: output of total HNO3 is before vmr is set to gas-phase. call outfld( 'HNO3_TOTAL', vmr(:ncol,:,hno3_ndx), ncol ,lchnk ) do k = 1,pver vmr(:,k,hno3_ndx) = hno3_gas(:,k) h2ovmr(:,k) = h2o_gas(:,k) vmr(:,k,h2o_ndx) = h2o_gas(:,k) wrk(:,k) = (h2ovmr(:,k) - wrk(:,k))*delt_inverse end do call outfld( 'QDSAD', wrk(:,:), ncol, lchnk ) ! call outfld( 'SAD_STRAT', strato_sad(:ncol,:), ncol, lchnk ) call outfld( 'SAD_SULFC', sad_strat(:,:,1), ncol, lchnk ) call outfld( 'SAD_LNAT', sad_strat(:,:,2), ncol, lchnk ) call outfld( 'SAD_ICE', sad_strat(:,:,3), ncol, lchnk ) ! call outfld( 'RAD_SULFC', radius_strat(:,:,1), ncol, lchnk ) call outfld( 'RAD_LNAT', radius_strat(:,:,2), ncol, lchnk ) call outfld( 'RAD_ICE', radius_strat(:,:,3), ncol, lchnk ) ! call outfld( 'HNO3_GAS', vmr(:ncol,:,hno3_ndx), ncol, lchnk ) call outfld( 'HNO3_STS', hno3_cond(:,:,1), ncol, lchnk ) call outfld( 'HNO3_NAT', hno3_cond(:,:,2), ncol, lchnk ) ! call outfld( 'HCL_TOTAL', vmr(:ncol,:,hcl_ndx), ncol, lchnk ) call outfld( 'HCL_GAS', hcl_gas (:,:), ncol ,lchnk ) call outfld( 'HCL_STS', hcl_cond(:,:), ncol ,lchnk ) !----------------------------------------------------------------------- ! ... call aerosol reaction rates !----------------------------------------------------------------------- call ratecon_sfstrat( ncol, invariants(:,:,indexm), pmid, tfld, & radius_strat(:,:,1), sad_strat(:,:,1), sad_strat(:,:,2), & sad_strat(:,:,3), h2ovmr, vmr, reaction_rates, & gprob_n2o5, gprob_cnt_hcl, gprob_cnt_h2o, gprob_bnt_h2o, & gprob_hocl_hcl, gprob_hobr_hcl, wtper ) call outfld( 'GAMMA_HET1', gprob_n2o5 (:ncol,:), ncol, lchnk ) call outfld( 'GAMMA_HET2', gprob_cnt_h2o (:ncol,:), ncol, lchnk ) call outfld( 'GAMMA_HET3', gprob_bnt_h2o (:ncol,:), ncol, lchnk ) call outfld( 'GAMMA_HET4', gprob_cnt_hcl (:ncol,:), ncol, lchnk ) call outfld( 'GAMMA_HET5', gprob_hocl_hcl(:ncol,:), ncol, lchnk ) call outfld( 'GAMMA_HET6', gprob_hobr_hcl(:ncol,:), ncol, lchnk ) call outfld( 'WTPER', wtper (:ncol,:), ncol, lchnk ) endif stratochem ! NOTE: For gas-phase solver only. ! ratecon_sfstrat needs total hcl. if (hcl_ndx>0) then vmr(:,:,hcl_ndx) = hcl_gas(:,:) endif !----------------------------------------------------------------------- ! ... Set the column densities at the upper boundary !----------------------------------------------------------------------- call set_ub_col( col_delta, vmr, invariants, pint(:,1), pdel, ncol, lchnk) !----------------------------------------------------------------------- ! ... Set rates for "tabular" and user specified reactions !----------------------------------------------------------------------- call setrxt( reaction_rates, tfld, invariants(1,1,indexm), ncol ) sulfate(:,:) = 0._r8 if ( .not. carma_hetchem_feedback ) then if( so4_ndx < 1 ) then ! get offline so4 field if not prognostic call sulf_interp( ncol, lchnk, sulfate ) else sulfate(:,:) = vmr(:,:,so4_ndx) endif endif !----------------------------------------------------------------- ! ... zero out sulfate above tropopause !----------------------------------------------------------------- do k = 1, pver do i = 1, ncol if (k < troplevchem(i)) then sulfate(i,k) = 0.0_r8 end if end do end do call outfld( 'SULF_TROP', sulfate(:ncol,:), ncol, lchnk ) !----------------------------------------------------------------- ! ... compute the relative humidity !----------------------------------------------------------------- call qsat(tfld(:ncol,:), pmid(:ncol,:), satv, satq) do k = 1,pver relhum(:,k) = .622_r8 * h2ovmr(:,k) / satq(:,k) relhum(:,k) = max( 0._r8,min( 1._r8,relhum(:,k) ) ) end do cwat(:ncol,:pver) = cldw(:ncol,:pver) call usrrxt( reaction_rates, tfld, ion_temp_fld, ele_temp_fld, invariants, h2ovmr, & pmid, invariants(:,:,indexm), sulfate, mmr, relhum, strato_sad, & troplevchem, dlats, ncol, sad_trop, reff, cwat, mbar, pbuf ) call outfld( 'SAD_TROP', sad_trop(:ncol,:), ncol, lchnk ) ! Add trop/strat components of SAD for output sad_trop(:ncol,:)=sad_trop(:ncol,:)+strato_sad(:ncol,:) call outfld( 'SAD_AERO', sad_trop(:ncol,:), ncol, lchnk ) ! Add trop/strat components of effective radius for output reff(:ncol,:)=reff(:ncol,:)+reff_strat(:ncol,:) call outfld( 'REFF_AERO', reff(:ncol,:), ncol, lchnk ) if (het1_ndx>0) then call outfld( 'het1_total', reaction_rates(:,:,het1_ndx), ncol, lchnk ) endif if (ghg_chem) then call ghg_chem_set_rates( reaction_rates, latmapback, zen_angle, ncol, lchnk ) endif do i = phtcnt+1,rxntot call outfld( rxn_names(i-phtcnt), reaction_rates(:,:,i), ncol, lchnk ) enddo call adjrxt( reaction_rates, invariants, invariants(1,1,indexm), ncol,pver ) !----------------------------------------------------------------------- ! ... Compute the photolysis rates at time = t(n+1) !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ... Set the column densities !----------------------------------------------------------------------- call setcol( col_delta, col_dens, vmr, pdel, ncol ) !----------------------------------------------------------------------- ! ... Calculate the photodissociation rates !----------------------------------------------------------------------- esfact = 1._r8 call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr , & delta, esfact ) if ( xactive_prates ) then if ( dst_ndx > 0 ) then dust_vmr(:ncol,:,1:ndust) = vmr(:ncol,:,dst_ndx:dst_ndx+ndust-1) else dust_vmr(:ncol,:,:) = 0._r8 endif !----------------------------------------------------------------- ! ... compute the photolysis rates !----------------------------------------------------------------- call xactive_photo( reaction_rates, vmr, tfld, cwat, cldfr, & pmid, zmidr, col_dens, zen_angle, asdir, & invariants(1,1,indexm), ps, ts, & esfact, relhum, dust_vmr, dt_diag, fracday, ncol, lchnk ) call outfld('DTCBS', dt_diag(:ncol,1), ncol, lchnk ) call outfld('DTOCS', dt_diag(:ncol,2), ncol, lchnk ) call outfld('DTSO4', dt_diag(:ncol,3), ncol, lchnk ) call outfld('DTANT', dt_diag(:ncol,4), ncol, lchnk ) call outfld('DTSAL', dt_diag(:ncol,5), ncol, lchnk ) call outfld('DTDUST', dt_diag(:ncol,6), ncol, lchnk ) call outfld('DTSOA', dt_diag(:ncol,7), ncol, lchnk ) call outfld('DTTOTAL', dt_diag(:ncol,8), ncol, lchnk ) call outfld('FRACDAY', fracday(:ncol), ncol, lchnk ) else !----------------------------------------------------------------- ! ... lookup the photolysis rates from table !----------------------------------------------------------------- call table_photo( reaction_rates, pmid, pdel, tfld, zmid, zint, & col_dens, zen_angle, asdir, cwat, cldfr, & esfact, vmr, invariants, ncol, lchnk, pbuf ) endif do i = 1,phtcnt call outfld( pht_names(i), reaction_rates(:ncol,:,i), ncol, lchnk ) call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk ) enddo !----------------------------------------------------------------------- ! ... Adjust the photodissociation rates !----------------------------------------------------------------------- call O1D_to_2OH_adj( reaction_rates, invariants, invariants(:,:,indexm), ncol, tfld ) call phtadj( reaction_rates, invariants, invariants(:,:,indexm), ncol,pver ) !----------------------------------------------------------------------- ! ... Compute the extraneous frcing at time = t(n+1) !----------------------------------------------------------------------- if ( o2_ndx > 0 .and. o_ndx > 0 ) then do k = 1,pver o2mmr(:ncol,k) = mmr(:ncol,k,o2_ndx) ommr(:ncol,k) = mmr(:ncol,k,o_ndx) end do endif call setext( extfrc, zint, zintr, cldtop, & zmid, lchnk, tfld, o2mmr, ommr, & pmid, mbar, rlats, calday, ncol, rlons, pbuf ) ! include forcings from fire emissions ... call fire_emissions_vrt( ncol, lchnk, zint, fire_sflx, fire_ztop, extfrc ) do m = 1,extcnt if( m /= synoz_ndx .and. m /= aoa_nh_ext_ndx ) then do k = 1,pver extfrc(:ncol,k,m) = extfrc(:ncol,k,m) / invariants(:ncol,k,indexm) end do endif call outfld( extfrc_name(m), extfrc(:ncol,:,m), ncol, lchnk ) end do !----------------------------------------------------------------------- ! ... Form the washout rates !----------------------------------------------------------------------- if ( gas_wetdep_method=='MOZ' ) then call sethet( het_rates, pmid, zmid, phis, tfld, & cmfdqr, prain, nevapr, delt, invariants(:,:,indexm), & vmr, ncol, lchnk ) if (.not. convproc_do_aer) then call het_diags( het_rates(:ncol,:,:), mmr(:ncol,:,:), pdel(:ncol,:), lchnk, ncol ) endif else het_rates = 0._r8 end if ! ! CCMI ! ! set loss to below the tropopause only ! if ( st80_25_tau_ndx > 0 ) then do i = 1,ncol reaction_rates(i,1:troplev(i),st80_25_tau_ndx) = 0._r8 enddo end if ! do i = phtcnt+1,rxt_tag_cnt call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk ) enddo if ( has_linoz_data ) then ltrop_sol(:ncol) = troplev(:ncol) else ltrop_sol(:ncol) = 0 ! apply solver to all levels endif ! save h2so4 before gas phase chem (for later new particle nucleation) if (ndx_h2so4 > 0) then del_h2so4_gasprod(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) else del_h2so4_gasprod(:,:) = 0.0_r8 endif vmr0(:ncol,:,:) = vmr(:ncol,:,:) ! mixing ratios before chemistry changes !======================================================================= ! ... Call the class solution algorithms !======================================================================= !----------------------------------------------------------------------- ! ... Solve for "Explicit" species !----------------------------------------------------------------------- call exp_sol( vmr, reaction_rates, het_rates, extfrc, delt, invariants(1,1,indexm), ncol, lchnk, ltrop_sol ) !----------------------------------------------------------------------- ! ... Solve for "Implicit" species !----------------------------------------------------------------------- if ( has_strato_chem ) wrk(:,:) = vmr(:,:,h2o_ndx) call t_startf('imp_sol') ! call imp_sol( vmr, reaction_rates, het_rates, extfrc, delt, & ncol,pver, lchnk, prod_out, loss_out ) call t_stopf('imp_sol') call chem_prod_loss_diags_out( ncol, lchnk, vmr, reaction_rates, prod_out, loss_out, invariants(:ncol,:,indexm) ) if( h2o_ndx>0) call outfld( 'H2O_GAS', vmr(1,1,h2o_ndx), ncol ,lchnk ) ! reset O3S to O3 in the stratosphere ... if ( o3_ndx > 0 .and. o3s_ndx > 0 ) then do i = 1,ncol vmr(i,1:troplev(i),o3s_ndx) = vmr(i,1:troplev(i),o3_ndx) end do end if if (convproc_do_aer) then call vmr2mmr( vmr(:ncol,:,:), mmr_new(:ncol,:,:), mbar(:ncol,:), ncol ) ! mmr_new = average of mmr values before and after imp_sol mmr_new(:ncol,:,:) = 0.5_r8*( mmr(:ncol,:,:) + mmr_new(:ncol,:,:) ) call het_diags( het_rates(:ncol,:,:), mmr_new(:ncol,:,:), pdel(:ncol,:), lchnk, ncol ) endif ! save h2so4 change by gas phase chem (for later new particle nucleation) if (ndx_h2so4 > 0) then del_h2so4_gasprod(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) - del_h2so4_gasprod(1:ncol,:) endif ! ! Aerosol processes ... ! call aero_model_gasaerexch( imozart-1, ncol, lchnk, troplevchem, delt, reaction_rates, & tfld, pmid, pdel, mbar, relhum, & zm, qh2o, cwat, cldfr, ncldwtr, & invariants(:,:,indexm), invariants, del_h2so4_gasprod, & vmr0, vmr, pbuf ) if ( has_strato_chem ) then wrk(:ncol,:) = (vmr(:ncol,:,h2o_ndx) - wrk(:ncol,:))*delt_inverse call outfld( 'QDCHEM', wrk(:ncol,:), ncol, lchnk ) call outfld( 'HNO3_GAS', vmr(:ncol,:,hno3_ndx), ncol ,lchnk ) !----------------------------------------------------------------------- ! ... aerosol settling ! first settle hno3(2) using radius ice ! secnd settle hno3(3) using radius large nat !----------------------------------------------------------------------- wrk(:,:) = vmr(:,:,h2o_ndx) #ifdef ALT_SETTL where( h2o_cond(:,:) > 0._r8 ) settl_rad(:,:) = radius_strat(:,:,3) elsewhere settl_rad(:,:) = 0._r8 endwhere call strat_aer_settling( invariants(1,1,indexm), pmid, delt, zmid, tfld, & hno3_cond(1,1,2), settl_rad, ncol, lchnk, 1 ) where( h2o_cond(:,:) == 0._r8 ) settl_rad(:,:) = radius_strat(:,:,2) elsewhere settl_rad(:,:) = 0._r8 endwhere call strat_aer_settling( invariants(1,1,indexm), pmid, delt, zmid, tfld, & hno3_cond(1,1,2), settl_rad, ncol, lchnk, 2 ) #else call strat_aer_settling( invariants(1,1,indexm), pmid, delt, zmid, tfld, & hno3_cond(1,1,2), radius_strat(1,1,2), ncol, lchnk, 2 ) #endif !----------------------------------------------------------------------- ! ... reform total hno3 and hcl = gas + all condensed !----------------------------------------------------------------------- ! NOTE: vmr for hcl and hno3 is gas-phase at this point. ! hno3_cond(:,k,1) = STS; hno3_cond(:,k,2) = NAT do k = 1,pver vmr(:,k,hno3_ndx) = vmr(:,k,hno3_ndx) + hno3_cond(:,k,1) & + hno3_cond(:,k,2) vmr(:,k,hcl_ndx) = vmr(:,k,hcl_ndx) + hcl_cond(:,k) end do wrk(:,:) = (vmr(:,:,h2o_ndx) - wrk(:,:))*delt_inverse call outfld( 'QDSETT', wrk(:,:), ncol, lchnk ) endif ! ! LINOZ ! if ( do_lin_strat_chem ) then call lin_strat_chem_solve( ncol, lchnk, vmr(:,:,o3_ndx), col_dens(:,:,1), tfld, zen_angle, pmid, delt, rlats, troplev ) end if !----------------------------------------------------------------------- ! ... Check for negative values and reset to zero !----------------------------------------------------------------------- call negtrc( 'After chemistry ', vmr, ncol ) !----------------------------------------------------------------------- ! ... Set upper boundary mmr values !----------------------------------------------------------------------- call set_fstrat_vals( vmr, pmid, pint, troplev, calday, ncol,lchnk ) !----------------------------------------------------------------------- ! ... Set fixed lower boundary mmr values !----------------------------------------------------------------------- call flbc_set( vmr, ncol, lchnk, map2chm ) !----------------------------------------------------------------------- ! set NOy UBC !----------------------------------------------------------------------- call noy_ubc_set( lchnk, ncol, vmr ) if ( ghg_chem ) then call ghg_chem_set_flbc( vmr, ncol ) endif !----------------------------------------------------------------------- ! force ion/electron balance -- ext forcings likely do not conserve charge !----------------------------------------------------------------------- call charge_balance( ncol, vmr ) !----------------------------------------------------------------------- ! ... Xform from vmr to mmr !----------------------------------------------------------------------- call vmr2mmr( vmr(:ncol,:,:), mmr_tend(:ncol,:,:), mbar(:ncol,:), ncol ) call set_short_lived_species( mmr_tend, lchnk, ncol, pbuf ) !----------------------------------------------------------------------- ! ... Form the tendencies !----------------------------------------------------------------------- do m = 1,gas_pcnst mmr_new(:ncol,:,m) = mmr_tend(:ncol,:,m) mmr_tend(:ncol,:,m) = (mmr_tend(:ncol,:,m) - mmr(:ncol,:,m))*delt_inverse enddo do m = 1,pcnst n = map2chm(m) if( n > 0 ) then qtend(:ncol,:,m) = qtend(:ncol,:,m) + mmr_tend(:ncol,:,n) end if end do tvs(:ncol) = tfld(:ncol,pver) * (1._r8 + qh2o(:ncol,pver)) sflx(:,:) = 0._r8 call get_ref_date(yr, mon, day, sec) ncdate = yr*10000 + mon*100 + day wind_speed(:ncol) = sqrt( ufld(:ncol,pver)*ufld(:ncol,pver) + vfld(:ncol,pver)*vfld(:ncol,pver) ) prect(:ncol) = precc(:ncol) + precl(:ncol) if ( drydep_method == DD_XLND ) then soilw = -99 call drydep( ocnfrac, icefrac, ncdate, ts, ps, & wind_speed, qh2o(:,pver), tfld(:,pver), pmid(:,pver), prect, & snowhland, fsds, depvel, sflx, mmr, & tvs, soilw, relhum(:,pver:pver), ncol, lonndx, latndx, lchnk ) else if ( drydep_method == DD_XATM ) then table_soilw = has_drydep( 'H2' ) .or. has_drydep( 'CO' ) if( .not. dyn_soilw .and. table_soilw ) then call set_soilw( soilw, lchnk, calday ) end if call drydep( ncdate, ts, ps, & wind_speed, qh2o(:,pver), tfld(:,pver), pmid(:,pver), prect, & snowhland, fsds, depvel, sflx, mmr, & tvs, soilw, relhum(:,pver:pver), ncol, lonndx, latndx, lchnk ) else if ( drydep_method == DD_TABL ) then call drydep( calday, ts, zen_angle, & depvel, sflx, mmr, pmid(:,pver), & tvs, ncol, icefrac, ocnfrac, lchnk ) endif drydepflx(:,:) = 0._r8 do m = 1,pcnst n = map2chm( m ) if ( n > 0 ) then cflx(:ncol,m) = cflx(:ncol,m) - sflx(:ncol,n) drydepflx(:ncol,m) = sflx(:ncol,n) wetdepflx_diag(:ncol,n) = wetdepflx(:ncol,m) endif end do call chm_diags( lchnk, ncol, vmr(:ncol,:,:), mmr_new(:ncol,:,:), & reaction_rates(:ncol,:,:), invariants(:ncol,:,:), depvel(:ncol,:), sflx(:ncol,:), & mmr_tend(:ncol,:,:), pdel(:ncol,:), pmid(:ncol,:), troplev(:ncol), wetdepflx_diag(:ncol,:), & nhx_nitrogen_flx(:ncol), noy_nitrogen_flx(:ncol) ) call rate_diags_calc( reaction_rates(:,:,:), vmr(:,:,:), invariants(:,:,indexm), ncol, lchnk ) ! ! jfl ! ! surface vmr ! if ( pm25_srf_diag ) then pm25(:ncol) = mmr_new(:ncol,pver,cb1_ndx) & + mmr_new(:ncol,pver,cb2_ndx) & + mmr_new(:ncol,pver,oc1_ndx) & + mmr_new(:ncol,pver,oc2_ndx) & + mmr_new(:ncol,pver,dst1_ndx) & + mmr_new(:ncol,pver,dst2_ndx) & + mmr_new(:ncol,pver,sslt1_ndx) & + mmr_new(:ncol,pver,sslt2_ndx) & + mmr_new(:ncol,pver,soa_ndx) & + mmr_new(:ncol,pver,so4_ndx) call outfld('PM25_SRF',pm25(:ncol) , ncol, lchnk ) endif if ( pm25_srf_diag_soa ) then pm25(:ncol) = mmr_new(:ncol,pver,cb1_ndx) & + mmr_new(:ncol,pver,cb2_ndx) & + mmr_new(:ncol,pver,oc1_ndx) & + mmr_new(:ncol,pver,oc2_ndx) & + mmr_new(:ncol,pver,dst1_ndx) & + mmr_new(:ncol,pver,dst2_ndx) & + mmr_new(:ncol,pver,sslt1_ndx) & + mmr_new(:ncol,pver,sslt2_ndx) & + mmr_new(:ncol,pver,soam_ndx) & + mmr_new(:ncol,pver,soai_ndx) & + mmr_new(:ncol,pver,soat_ndx) & + mmr_new(:ncol,pver,soab_ndx) & + mmr_new(:ncol,pver,soax_ndx) & + mmr_new(:ncol,pver,so4_ndx) call outfld('PM25_SRF',pm25(:ncol) , ncol, lchnk ) endif ! ! call outfld('Q_SRF',qh2o(:ncol,pver) , ncol, lchnk ) call outfld('U_SRF',ufld(:ncol,pver) , ncol, lchnk ) call outfld('V_SRF',vfld(:ncol,pver) , ncol, lchnk ) ! if (.not.sad_pbf_ndx>0) then deallocate(strato_sad) endif end subroutine gas_phase_chemdr end module mo_gas_phase_chemdr