chem/module_chem_cup.F

   1 ! *** changes needed involving soa
   2 !
   3 ! include soa in calculation of oa_a_1to4_ic_cup and oa_cw_1to4_ic_cup
   4 !
   5 ! account for soa hygroscopicity for droplet activation (??)
   6
   7
   8 !**********************************************************************************
   9 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
  10 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
  11 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
  12 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
  13 !
  14 ! MOSAIC module: see module_mosaic_driver.F for references and terms of use
  15
  16 ! *** NOTE ***
  17 ! This module should in theory work with any chem_opt package that uses MOSAIC
  18 !and has cloud-borne aerosols. However, it was only tested with chem_opt=203
  19 !(saprc99_mosaic_8bin_vbs2_aq_kpp) during implementation into WRF-Chem in
  20 !April 2017 at PNNL.
  21
  22
  23 !**********************************************************************************
  24
  25 !  file module_chem_cup.F
  26       module module_chem_cup
  27
  28
  29       implicit none
  30
  31
  32       private
  33
  34       public :: chem_cup_driver
  35
  36       logical, parameter :: do_deep = .true.
  37
  38       integer, parameter :: r4=4
  39       integer, parameter :: r8=8
  40
  41       integer, parameter :: mode_chemcup_timeinteg = 2
  42       ! when mode_chemcup_timeinteg = 1, simulate the entire cumulus cloud fraction
  43       !    with multiple time substeps of dt = dtstepc
  44       ! when mode_chemcup_timeinteg = 2, simulate the entire cumulus cloud fraction
  45       !    with a single timestep of dt = tcloud_cup, then set tcloud_cup negative
  46       !    to signify that chem_cup calcs. are done for this cloud
  47       !
  48       ! *** note 1 - mode_chemcup_timeinteg = 2 is recommended
  49       ! *** note 2 - mode_chemcup_timeinteg = 1 may need lower values for the various ..._smallaa parameters
  50
  51       real(r8), parameter :: cldfra_ls_testvalue = -0.5
  52       ! *** used for testing *** -- when cldfra_ls_testvalue >= 0.0, it overrides cldfra
  53       real(r8), parameter :: wact_cup_testvalue = -0.1
  54       ! *** used for testing *** -- when wact_cup_testvalue > 0.0, it overrides wact_cup
  55
  56       real(r8), parameter :: air_outflow_limit = 0.90_r8
  57       ! maximum faction of layer mass that can be transported out (to other layers) in a transport sub-time-step
  58
  59       real(r8), parameter :: af_cucld_smallaa = 0.003_r8  ! "cutoff" value for cumulus cloud fractional area
  60       real(r8), parameter :: af_cucld_maxaa = 0.8_r8      !  maximum value for cumulus cloud fractional area
  61
  62 ! subarea-average vertical mass fluxes (kg/m2/s) smaller than
  63 !    aw_up_smallaa*rhoair are treated as zero
  64 ! note that with a*w = 3e-5 m/s, dz over 1 hour = 0.11 m which is small
  65       real(r8), parameter :: aw_up_smallaa = 3.0e-5_r8   ! m/s
  66 ! maximum expected updraft
  67       real(r8), parameter :: w_up_maxaa = 50.0_r8   ! m/s
  68 ! updraft fractional areas below af_up_smallaa are ignored
  69       real(r8), parameter :: af_up_smallaa = aw_up_smallaa/w_up_maxaa
  70
  71       real(r8), parameter :: af_up_maxaa = 0.2_r8        !  maximum value for cumulus updraft fractional area
  72
  73       real(r8), parameter :: qci_incu_smallaa = 1.0e-6_r8  ! "cutoff" value for in-cloud   cloud-ice   (kg/kg)
  74       real(r8), parameter :: qci_inup_smallaa = 1.0e-6_r8  ! "cutoff" value for in-updraft cloud-ice   (kg/kg)
  75       real(r8), parameter :: qcw_incu_smallaa = 1.0e-6_r8  ! "cutoff" value for in-cloud   cloud-water (kg/kg)
  76       real(r8), parameter :: qcw_inup_smallaa = 1.0e-5_r8  ! "cutoff" value for in-updraft cloud-water (kg/kg)
  77
  78       real(r8), parameter :: tau_active_smallaa   = 30.0_r8  ! "cutoff" value for active cloud calculations (s)
  79       real(r8), parameter :: tau_inactive_smallaa = 30.0_r8  ! "cutoff" value for inactive cloud calculations (s)
  80
  81
  82       contains
  83
  84
  85 !----------------------------------------------------------------------
  86       subroutine chem_cup_driver(                                     &
  87             grid_id, ktau, ktauc, dtstep, dtstepc, config_flags,      &
  88             t_phy, p_phy, rho_phy, alt, dz8w, zmid, z_at_w,           &
  89             moist, cldfra, ph_no2,                                    &
  90             chem,                                                     &
  91             chem_cupflag, cupflag, shall, tcloud_cup, nca, wact_cup,  &
  92             cldfra_cup, updfra_cup, qc_ic_cup, qc_iu_cup,             &
  93             mfup_cup, mfup_ent_cup, mfdn_cup, mfdn_ent_cup,           &
  94             fcvt_qc_to_pr_cup, fcvt_qc_to_qi_cup, fcvt_qi_to_pr_cup,  &
  95             co_a_ic_cup, hno3_a_ic_cup,                               &
  96             so4_a_1to4_ic_cup, so4_cw_1to4_ic_cup,                      &
  97             nh4_a_1to4_ic_cup, nh4_cw_1to4_ic_cup,                      &
  98             no3_a_1to4_ic_cup, no3_cw_1to4_ic_cup,                      &
  99             oa_a_1to4_ic_cup,  oa_cw_1to4_ic_cup,                       &
 100             oin_a_1to4_ic_cup,  oin_cw_1to4_ic_cup,                     &
 101             bc_a_1to4_ic_cup,  bc_cw_1to4_ic_cup,                       &
 102             na_a_1to4_ic_cup,  na_cw_1to4_ic_cup,                       &
 103             cl_a_1to4_ic_cup,  cl_cw_1to4_ic_cup,                       &
 104             water_1to4_ic_cup,                                          &
 105             so4_a_5to6_ic_cup, so4_cw_5to6_ic_cup,                      &
 106             nh4_a_5to6_ic_cup, nh4_cw_5to6_ic_cup,                      &
 107             no3_a_5to6_ic_cup, no3_cw_5to6_ic_cup,                      &
 108             oa_a_5to6_ic_cup,  oa_cw_5to6_ic_cup,                       &
 109             oin_a_5to6_ic_cup,  oin_cw_5to6_ic_cup,                     &
 110             bc_a_5to6_ic_cup,  bc_cw_5to6_ic_cup,                       &
 111             na_a_5to6_ic_cup,  na_cw_5to6_ic_cup,                       &
 112             cl_a_5to6_ic_cup,  cl_cw_5to6_ic_cup,                       &
 113             water_5to6_ic_cup,                                          &
 114             ids,ide, jds,jde, kds,kde,                                &
 115             ims,ime, jms,jme, kms,kme,                                &
 116             its,ite, jts,jte, kts,kte                                 )
 117 !
 118 !  convective cloud processing (vertical transport,
 119 !     activation/resuspension, cloud chemistry, wet removal eventually)
 120 !     by **cup** convective clouds
 121 !
 122 !  currently only works with following chem_opt packages
 123 !     CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ )
 124 !  currently does not do wet removal
 125 !
 126
 127       use module_configure
 128       use module_state_description
 129       use module_model_constants
 130       use module_scalar_tables, only:  chem_dname_table
 131
 132 !----------------------------------------------------------------------
 133 ! arguments
 134
 135       type(grid_config_rec_type), intent(in) :: config_flags
 136
 137       integer, intent(in) ::              &
 138             ids,ide, jds,jde, kds,kde,    &
 139             ims,ime, jms,jme, kms,kme,    &
 140             its,ite, jts,jte, kts,kte
 141
 142       integer, intent(in) ::              &
 143             grid_id, ktau, ktauc
 144
 145       integer, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: &
 146             chem_cupflag               !  1 if cupflag=.true. and chem_cup calcs were successful
 147                                        ! -1 if cupflag=.true. and chem_cup calcs failed
 148                                        !  0 if cupflag=.false.
 149
 150       real, intent(in) :: dtstep, dtstepc  ! model and chemistry time-steps (s)
 151
 152       real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: &
 153             t_phy,                   & ! temp (K)
 154             p_phy,                   & ! pressure (Pa)
 155             rho_phy,                 & ! moist air density (kg/m3)
 156             alt,                     & ! dry air 1/density (m3/kg)
 157             dz8w,                    & ! layer thickness (m)
 158             zmid,                    & ! height at layer center (m)
 159             z_at_w,                  & ! height at layer boundary (m)
 160             cldfra,                  & ! grid-resolved cloud fraction
 161             ph_no2                     ! no2 photolysis rate (1/s)
 162
 163 ! advected moisture variables
 164       real, dimension( ims:ime, kms:kme, jms:jme, num_moist ), intent(in) :: &
 165             moist
 166
 167 ! all advected chemical species
 168       real, dimension( ims:ime, kms:kme, jms:jme, num_chem ), intent(inout) :: &
 169             chem
 170
 171       logical, dimension( ims:ime, jms:jme ), intent(in) :: &
 172             cupflag                    ! .true. if cup convection is present in column
 173
 174       real, dimension( ims:ime, jms:jme ), intent(in) :: &
 175             nca,                     & ! time remaining for this cloud (s)
 176             shall,                   & ! cumulus type:  0=deep,  1=shallow,  2=none
 177             wact_cup                   ! vertical velocity for cloud-base aerosol activation (m/s)
 178
 179       real, dimension( ims:ime, jms:jme ), intent(inout) :: &
 180             tcloud_cup                 ! cumulus cloud duration (s)
 181
 182       real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: &
 183             cldfra_cup,              & ! cumulus cloud fractional area (-)
 184             updfra_cup,              & ! cumulus updraft fractional area (-)
 185             qc_ic_cup,               & ! cloud-water in cumulus cloud (kg/kg)
 186             qc_iu_cup,               & ! cloud-water in updraft (kg/kg)
 187             mfup_cup,                & ! updraft mass-flux (kg/m2/s)
 188             mfup_ent_cup,            & ! updraft mass-flux change from entrainment (kg/m2/s)
 189             mfdn_cup,                & ! downdraft mass-flux (kg/m2/s)
 190             mfdn_ent_cup,            & ! downdraft mass-flux change from entrainment (kg/m2/s)
 191             fcvt_qc_to_pr_cup,       & ! fraction of cloud-water converted to precip    as air move thru updraft layer (-)
 192             fcvt_qc_to_qi_cup,       & ! fraction of cloud-water converted to cloud-ice as air move thru updraft layer (-)
 193             fcvt_qi_to_pr_cup          ! fraction of cloud-ice   converted to precip    as air move thru updraft layer (-)
 194
 195 ! interstitial and cloudborne mixing ratios within the convective cloud,
 196 ! summed over size bins that the ams can see (Dp < 625 nm)
 197       real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: &
 198             co_a_ic_cup, hno3_a_ic_cup,               &
 199             so4_a_1to4_ic_cup, so4_cw_1to4_ic_cup,  &
 200             nh4_a_1to4_ic_cup, nh4_cw_1to4_ic_cup,  &
 201             no3_a_1to4_ic_cup, no3_cw_1to4_ic_cup,  &
 202             oa_a_1to4_ic_cup,  oa_cw_1to4_ic_cup,   &
 203             oin_a_1to4_ic_cup,  oin_cw_1to4_ic_cup,   &
 204             bc_a_1to4_ic_cup,  bc_cw_1to4_ic_cup,   &
 205             na_a_1to4_ic_cup,  na_cw_1to4_ic_cup,   &
 206             cl_a_1to4_ic_cup,  cl_cw_1to4_ic_cup,   &
 207             water_1to4_ic_cup,                       &
 208             so4_a_5to6_ic_cup, so4_cw_5to6_ic_cup,  &
 209             nh4_a_5to6_ic_cup, nh4_cw_5to6_ic_cup,  &
 210             no3_a_5to6_ic_cup, no3_cw_5to6_ic_cup,  &
 211             oa_a_5to6_ic_cup,  oa_cw_5to6_ic_cup,   &
 212             oin_a_5to6_ic_cup,  oin_cw_5to6_ic_cup,   &
 213             bc_a_5to6_ic_cup,  bc_cw_5to6_ic_cup,   &
 214             na_a_5to6_ic_cup,  na_cw_5to6_ic_cup,   &
 215             cl_a_5to6_ic_cup,  cl_cw_5to6_ic_cup,   &
 216             water_5to6_ic_cup
 217
 218
 219
 220 !-----------------------------------------------------------------
 221 ! local variables
 222       integer :: iok
 223       integer :: ii, jj, kk
 224
 225       character(len=12) :: chem_name(num_chem)
 226
 227 !-----------------------------------------------------------------
 228
 229 ! check for correct options
 230       if ( config_flags%chem_conv_tr <= 0 .or. &
 231            config_flags%cu_physics /= kfcupscheme ) then
 232          call wrf_debug( 15, 'chem_cup_driver skipped because - ' // &
 233                              'chem_conv_tr or cu_physics' )
 234          return
 235       end if
 236
 237       chem_name(1:num_chem) = chem_dname_table(grid_id,1:num_chem)
 238
 239
 240       chem_opt_select: select case(config_flags%chem_opt)
 241
 242          case ( CBMZ_MOSAIC_4BIN, CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ, &
 243                 SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP ) !BSINGH(12/11/13): Got rid of SAPRC99_MOSAIC_4BIN_VBS2_AQ_KPP pkg statement
 244 !-------------------------------------------------------------------------------
 245 ! Below lines commented by Manish Shrivastava to skip CUP chemistry
 246 !------------------------------------------------------------------------------
 247             call wrf_debug(15, &
 248                'chem_cup_driver calling mosaic_chem_cup_driver')
 249             call mosaic_chem_cup_driver(                                    &
 250                   grid_id, ktau, ktauc, dtstep, dtstepc, config_flags,      &
 251                   t_phy, p_phy, rho_phy, alt, dz8w, zmid, z_at_w,           &
 252                   moist, cldfra, ph_no2,                                    &
 253                   chem, chem_name,                                          &
 254                   chem_cupflag, cupflag, shall, tcloud_cup, nca, wact_cup,  &
 255                   cldfra_cup, updfra_cup, qc_ic_cup, qc_iu_cup,             &
 256                   mfup_cup, mfup_ent_cup, mfdn_cup, mfdn_ent_cup,           &
 257                   fcvt_qc_to_pr_cup, fcvt_qc_to_qi_cup, fcvt_qi_to_pr_cup,  &
 258                   co_a_ic_cup, hno3_a_ic_cup,                               &
 259                   so4_a_1to4_ic_cup, so4_cw_1to4_ic_cup,                      &
 260                   nh4_a_1to4_ic_cup, nh4_cw_1to4_ic_cup,                      &
 261                   no3_a_1to4_ic_cup, no3_cw_1to4_ic_cup,                      &
 262                   oa_a_1to4_ic_cup,  oa_cw_1to4_ic_cup,                       &
 263                   oin_a_1to4_ic_cup,  oin_cw_1to4_ic_cup,                     &
 264                   bc_a_1to4_ic_cup,  bc_cw_1to4_ic_cup,                       &
 265                   na_a_1to4_ic_cup,  na_cw_1to4_ic_cup,                       &
 266                   cl_a_1to4_ic_cup,  cl_cw_1to4_ic_cup,                       &
 267                   water_1to4_ic_cup,                                          &
 268                   so4_a_5to6_ic_cup, so4_cw_5to6_ic_cup,                      &
 269                   nh4_a_5to6_ic_cup, nh4_cw_5to6_ic_cup,                      &
 270                   no3_a_5to6_ic_cup, no3_cw_5to6_ic_cup,                      &
 271                   oa_a_5to6_ic_cup,  oa_cw_5to6_ic_cup,                       &
 272                   oin_a_5to6_ic_cup,  oin_cw_5to6_ic_cup,                     &
 273                   bc_a_5to6_ic_cup,  bc_cw_5to6_ic_cup,                       &
 274                   na_a_5to6_ic_cup,  na_cw_5to6_ic_cup,                       &
 275                   cl_a_5to6_ic_cup,  cl_cw_5to6_ic_cup,                       &
 276                   water_5to6_ic_cup,                                          &
 277                   ids,ide, jds,jde, kds,kde,                                &
 278                   ims,ime, jms,jme, kms,kme,                                &
 279                   its,ite, jts,jte, kts,kte, kte+1                          )
 280 !-----------------------------------------------------------------------------------------
 281 ! Above lines commented by Manish Shrivastava to skip CUP chemistry
 282 !-------------------------------------------------------------------------------------------
 283          case default
 284             chem_cupflag = 0
 285             call wrf_debug( 15, 'chem_cup_driver skipped because - ' // &
 286                                 'chem_opt' )
 287
 288       end select chem_opt_select
 289
 290       return
 291       end subroutine chem_cup_driver
 292
 293
 294 !----------------------------------------------------------------------
 295       subroutine mosaic_chem_cup_driver(                              &
 296             grid_id, ktau, ktauc, dtstep, dtstepc, config_flags,      &
 297             t_phy, p_phy, rho_phy, alt, dz8w, zmid, z_at_w,           &
 298             moist, cldfra, ph_no2,                                    &
 299             chem, chem_name,                                          &
 300             chem_cupflag, cupflag, shall, tcloud_cup, nca, wact_cup,  &
 301             cldfra_cup, updfra_cup, qc_ic_cup, qc_iu_cup,             &
 302             mfup_cup, mfup_ent_cup, mfdn_cup, mfdn_ent_cup,           &
 303             fcvt_qc_to_pr_cup, fcvt_qc_to_qi_cup, fcvt_qi_to_pr_cup,  &
 304             co_a_ic_cup, hno3_a_ic_cup,                               &
 305             so4_a_1to4_ic_cup, so4_cw_1to4_ic_cup,                      &
 306             nh4_a_1to4_ic_cup, nh4_cw_1to4_ic_cup,                      &
 307             no3_a_1to4_ic_cup, no3_cw_1to4_ic_cup,                      &
 308             oa_a_1to4_ic_cup,  oa_cw_1to4_ic_cup,                       &
 309             oin_a_1to4_ic_cup,  oin_cw_1to4_ic_cup,                     &
 310             bc_a_1to4_ic_cup,  bc_cw_1to4_ic_cup,                       &
 311             na_a_1to4_ic_cup,  na_cw_1to4_ic_cup,                       &
 312             cl_a_1to4_ic_cup,  cl_cw_1to4_ic_cup,                       &
 313             water_1to4_ic_cup,                                          &
 314             so4_a_5to6_ic_cup, so4_cw_5to6_ic_cup,                      &
 315             nh4_a_5to6_ic_cup, nh4_cw_5to6_ic_cup,                      &
 316             no3_a_5to6_ic_cup, no3_cw_5to6_ic_cup,                      &
 317             oa_a_5to6_ic_cup,  oa_cw_5to6_ic_cup,                       &
 318             oin_a_5to6_ic_cup,  oin_cw_5to6_ic_cup,                     &
 319             bc_a_5to6_ic_cup,  bc_cw_5to6_ic_cup,                       &
 320             na_a_5to6_ic_cup,  na_cw_5to6_ic_cup,                       &
 321             cl_a_5to6_ic_cup,  cl_cw_5to6_ic_cup,                       &
 322             water_5to6_ic_cup,                                          &
 323             ids,ide, jds,jde, kds,kde,                                &
 324             ims,ime, jms,jme, kms,kme,                                &
 325             its,ite, jts,jte, kts,kte, ktep1                          )
 326
 327 !  wet removal by grid-resolved precipitation
 328 !  scavenging of cloud-phase aerosols and gases by collection, freezing, ...
 329 !  scavenging of interstitial-phase aerosols by impaction
 330 !  scavenging of gas-phase gases by mass transfer and reaction
 331
 332       use module_configure
 333       use module_state_description
 334       use module_model_constants
 335
 336       use module_data_mosaic_asect, only:  &
 337          dcen_sect, &
 338          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
 339          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
 340          ai_phase, cw_phase, msectional, &
 341          massptr_aer, numptr_aer, &
 342          lptr_cl_aer, lptr_na_aer, lptr_nh4_aer, lptr_no3_aer, &
 343          lptr_oc_aer, lptr_oin_aer, lptr_so4_aer, lptr_bc_aer, &
 344          lptr_pcg1_b_c_aer,lptr_pcg1_b_o_aer,lptr_opcg1_b_c_aer, &
 345          lptr_opcg1_b_o_aer,lptr_pcg1_f_c_aer,lptr_pcg1_f_o_aer, &
 346          lptr_opcg1_f_c_aer, lptr_opcg1_f_o_aer, lptr_ant1_c_aer, &
 347          lptr_biog1_c_aer, waterptr_aer, &
 348          dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer
 349
 350 !----------------------------------------------------------------------
 351 ! arguments
 352
 353       type(grid_config_rec_type), intent(in) :: config_flags
 354
 355       integer, intent(in) ::              &
 356             ids,ide, jds,jde, kds,kde,    &
 357             ims,ime, jms,jme, kms,kme,    &
 358             its,ite, jts,jte, kts,kte, ktep1
 359
 360       integer, intent(in) ::              &
 361             grid_id, ktau, ktauc
 362
 363       integer, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: &
 364             chem_cupflag
 365
 366       real, intent(in) :: dtstep, dtstepc  ! model and chemistry time-steps (s)
 367
 368       real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: &
 369             t_phy,                   & ! temp (K)
 370             p_phy,                   & ! pressure (Pa)
 371             rho_phy,                 & ! moist air density (kg/m3)
 372             alt,                     & ! dry air 1/density (m3/kg)
 373             dz8w,                    & ! layer thickness (m)
 374             zmid,                    & ! height at layer center (m)
 375             z_at_w,                  & ! height at layer boundary (m)
 376             cldfra,                  & ! grid-resolved cloud fraction
 377             ph_no2                     ! no2 photolysis rate (1/s)
 378
 379 ! advected moisture variables
 380       real, dimension( ims:ime, kms:kme, jms:jme, num_moist ), intent(in) :: &
 381             moist
 382
 383 ! all advected chemical species
 384       real, dimension( ims:ime, kms:kme, jms:jme, num_chem ), intent(inout) :: &
 385             chem
 386
 387 ! interstitial and cloudborne mixing ratios within the convective cloud,
 388 ! summed over size bins that the ams can see (Dp < 625 nm)
 389       real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: &
 390             co_a_ic_cup, hno3_a_ic_cup,                                 &
 391             so4_a_1to4_ic_cup, so4_cw_1to4_ic_cup,                      &
 392             nh4_a_1to4_ic_cup, nh4_cw_1to4_ic_cup,                      &
 393             no3_a_1to4_ic_cup, no3_cw_1to4_ic_cup,                      &
 394             oa_a_1to4_ic_cup,  oa_cw_1to4_ic_cup,                       &
 395             oin_a_1to4_ic_cup,  oin_cw_1to4_ic_cup,                     &
 396             bc_a_1to4_ic_cup,  bc_cw_1to4_ic_cup,                       &
 397             na_a_1to4_ic_cup,  na_cw_1to4_ic_cup,                       &
 398             cl_a_1to4_ic_cup,  cl_cw_1to4_ic_cup,                       &
 399             water_1to4_ic_cup,                                          &
 400             so4_a_5to6_ic_cup, so4_cw_5to6_ic_cup,                      &
 401             nh4_a_5to6_ic_cup, nh4_cw_5to6_ic_cup,                      &
 402             no3_a_5to6_ic_cup, no3_cw_5to6_ic_cup,                      &
 403             oa_a_5to6_ic_cup,  oa_cw_5to6_ic_cup,                       &
 404             oin_a_5to6_ic_cup,  oin_cw_5to6_ic_cup,                     &
 405             bc_a_5to6_ic_cup,  bc_cw_5to6_ic_cup,                       &
 406             na_a_5to6_ic_cup,  na_cw_5to6_ic_cup,                       &
 407             cl_a_5to6_ic_cup,  cl_cw_5to6_ic_cup,                       &
 408             water_5to6_ic_cup
 409
 410
 411       character(len=12), intent(in) :: chem_name(num_chem)
 412
 413
 414       logical, dimension( ims:ime, jms:jme ), intent(in) :: &
 415             cupflag                    ! .true. if cup convection is present in column
 416
 417       real, dimension( ims:ime, jms:jme ), intent(in) :: &
 418             nca,                     & ! time remaining for this cloud (s)
 419             shall,                   & ! cumulus type:  0=deep,  1=shallow,  2=none
 420             wact_cup                   ! vertical velocity for cloud-base aerosol activation (m/s)
 421
 422       real, dimension( ims:ime, jms:jme ), intent(inout) :: &
 423             tcloud_cup                 ! cumulus cloud duration (s)
 424
 425       real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: &
 426             cldfra_cup,              & ! cumulus cloud fractional area (-)
 427             updfra_cup,              & ! cumulus updraft fractional area (-)
 428             qc_ic_cup,               & ! cloud-water in cumulus cloud (kg/kg)
 429             qc_iu_cup,               & ! cloud-water in updraft (kg/kg)
 430             mfup_cup,                & ! updraft mass-flux (kg/m2/s)
 431             mfup_ent_cup,            & ! updraft mass-flux change from entrainment (kg/m2/s)
 432             mfdn_cup,                & ! downdraft mass-flux (kg/m2/s)
 433             mfdn_ent_cup,            & ! downdraft mass-flux change from entrainment (kg/m2/s)
 434             fcvt_qc_to_pr_cup,       & ! fraction of cloud-water converted to precip    as air move thru updraft layer (-)
 435             fcvt_qc_to_qi_cup,       & ! fraction of cloud-water converted to cloud-ice as air move thru updraft layer (-)
 436             fcvt_qi_to_pr_cup          ! fraction of cloud-ice   converted to precip    as air move thru updraft layer (-)
 437
 438 !-----------------------------------------------------------------
 439 ! local variables
 440       integer :: aer_mech_id
 441       integer :: chem_cupflag_1d(kts:kte)
 442       integer :: i, icalcflagaa, idiagee, idiagff, ishall, isize, itype
 443       integer :: j
 444       integer :: k, kcldbot_1d, kcldtop_1d
 445       integer :: lundiag, lunerr
 446
 447       real(r8) :: chem_1d(kts:kte,num_chem)
 448       real(r8) :: chem_incu(kts:kte,num_chem)
 449       real(r8) :: cldfra_ls_1d(kts:kte)
 450       real(r8) :: cldfra_cup_1d(kts:kte)
 451       real(r8) :: dz_1d(kts:kte)
 452       real(r8) :: fcvt_qc_to_pr_cup_1d(kts:kte)
 453       real(r8) :: fcvt_qc_to_qi_cup_1d(kts:kte)
 454       real(r8) :: fcvt_qi_to_pr_cup_1d(kts:kte)
 455       real(r8) :: mfup_cup_1d(kts:ktep1)
 456       real(r8) :: mfup_ent_cup_1d(kts:kte)
 457       real(r8) :: mfdn_cup_1d(kts:ktep1)
 458       real(r8) :: mfdn_ent_cup_1d(kts:kte)
 459       real(r8) :: pcen_1d(kts:kte)
 460       real(r8) :: ph_no2_1d(kts:kte)
 461       real(r8) :: qc_ic_cup_1d(kts:kte)
 462       real(r8) :: qc_iu_cup_1d(kts:kte)
 463       real(r8) :: qi_ic_cup_1d(kts:kte)
 464       real(r8) :: qi_iu_cup_1d(kts:kte)
 465       real(r8) :: rhocen_1d(kts:kte)
 466       real(r8) :: tcloud_cup_1d
 467       real(r8) :: tcen_1d(kts:kte)
 468       real(r8) :: tmpa
 469       real(r8) :: updfra_cup_1d(kts:kte)
 470       real(r8) :: wact_cup_1d
 471       real(r8) :: zbnd_1d(kts:ktep1)
 472       real(r8) :: zcen_1d(kts:kte)
 473
 474
 475 !-----------------------------------------------------------------
 476
 477       lunerr  = 6
 478       lundiag = 6
 479       lundiag = 121
 480
 481       idiagff = 0 ; idiagee = 0
 482       if ((ide-ids <= 3) .and. (jde-jds <= 3)) then
 483          idiagff = 1  ! turn on diagnostics for single column runs
 484 !        idiagff = 0  ! (do this to turn off extra diagnostics)
 485       end if
 486
 487       aer_mech_id = 3
 488
 489       if (ktau <= 1) then
 490          write(*,'(a)')
 491          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 492             'cldfra_ls_testvalue, wact_cup_testvalue   ', &
 493              cldfra_ls_testvalue, wact_cup_testvalue
 494          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 495             'air_outflow_limit                         ', &
 496              air_outflow_limit
 497          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 498             'af_cucld_smallaa, af_cucld_maxaa          ', &
 499              af_cucld_smallaa, af_cucld_maxaa
 500          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 501             'aw_up_smallaa, w_up_maxaa                 ', &
 502              aw_up_smallaa, w_up_maxaa
 503          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 504             'af_up_smallaa, af_up_maxaa                ', &
 505              af_up_smallaa, af_up_maxaa
 506          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 507             'qci_incu_smallaa, qci_inup_smallaa        ', &
 508              qci_incu_smallaa, qci_inup_smallaa
 509          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 510             'qcw_incu_smallaa, qcw_inup_smallaa        ', &
 511              qcw_incu_smallaa, qcw_inup_smallaa
 512          write(*,'(2a,1p,4e12.4)') 'chemcup_control -- ', &
 513             'tau_active_smallaa, tau_inactive_smallaa  ', &
 514              tau_active_smallaa, tau_inactive_smallaa
 515          write(*,'(a,2i5/(a,3(i9,i5)))') &
 516             'chemcup_control -- grid_id, ktau', grid_id, ktau, &
 517             'chemcup_control -- d indices', ids,ide, jds,jde, kds,kde, &
 518             'chemcup_control -- m indices', ims,ime, jms,jme, kms,kme, &
 519             'chemcup_control -- e indices', its,ite, jts,jte, kts,kte
 520          write(*,'(a)')
 521       end if
 522
 523
 524 ! loop over grid columns and do chem_cup where shall = 1 (shallow cu) or 2 (deep)
 525 main_j_loop: &
 526       do j = jts, jte
 527 main_i_loop: &
 528       do i = its, ite
 529
 530       idiagee = 0
 531       if (idiagff > 0) then
 532          ! turn on diagnostics at i=j=1 for single column runs
 533          if (i==its .and. j==jts) idiagee = 1
 534       end if
 535
 536       if (idiagee > 0) write(*,'(a,i7,l5,i5,20x,4f10.1)') &
 537          'chcup_a20 ktau, cupflag, ishall,      tcloud, nca dtc', &
 538          ktau, cupflag(i,j), nint(shall(i,j)), tcloud_cup(i,j), nca(i,j), dtstepc
 539
 540
 541       ! icalcflagaa > 0 if chem_cup calculations will be done for current i,j and timestep
 542       icalcflagaa = 0
 543       if (abs(shall(i,j)-1.0) < 0.1) then
 544          ishall = 1
 545          icalcflagaa = 1
 546       else if (abs(shall(i,j)) < 0.1) then
 547          ishall = 0
 548          if ( do_deep ) icalcflagaa = 1
 549       end if
 550       ! no calculations if cloud has "expired"
 551       if (nca(i,j) < 0.01) icalcflagaa = 0
 552       ! no calculations if cloud lifetime is too small
 553       if (abs(tcloud_cup(i,j)) < tau_active_smallaa) icalcflagaa = 0
 554
 555       ! also no calculations when mode_chemcup_timeinteg=2 and the "one time only"
 556       !    calculations have already been done (indicated by tcloud_cup < 0)
 557       if (icalcflagaa > 0) then
 558          if (mode_chemcup_timeinteg == 2 .and. tcloud_cup(i,j) <= 0.0) then
 559             icalcflagaa = -1
 560          end if
 561       end if
 562
 563       if (icalcflagaa >= 0) then
 564          ! zero out these history variables at current i,j except when icalcflagaa = -1
 565          chem_cupflag(i,:,j) = 0
 566          co_a_ic_cup(i,:,j) = 0.0  ;  hno3_a_ic_cup(i,:,j) = 0.0
 567          so4_a_1to4_ic_cup(i,:,j) = 0.0 ; so4_cw_1to4_ic_cup(i,:,j) = 0.0
 568          nh4_a_1to4_ic_cup(i,:,j) = 0.0 ; nh4_cw_1to4_ic_cup(i,:,j) = 0.0
 569          no3_a_1to4_ic_cup(i,:,j) = 0.0 ; no3_cw_1to4_ic_cup(i,:,j) = 0.0
 570          oa_a_1to4_ic_cup (i,:,j) = 0.0 ;  oa_cw_1to4_ic_cup(i,:,j) = 0.0
 571          oin_a_1to4_ic_cup (i,:,j) = 0.0 ;  oin_cw_1to4_ic_cup(i,:,j) = 0.0
 572          bc_a_1to4_ic_cup (i,:,j) = 0.0 ;  bc_cw_1to4_ic_cup(i,:,j) = 0.0
 573          na_a_1to4_ic_cup (i,:,j) = 0.0 ;  na_cw_1to4_ic_cup(i,:,j) = 0.0
 574          cl_a_1to4_ic_cup (i,:,j) = 0.0 ;  cl_cw_1to4_ic_cup(i,:,j) = 0.0
 575          water_1to4_ic_cup (i,:,j) = 0.0
 576          so4_a_5to6_ic_cup(i,:,j) = 0.0 ; so4_cw_5to6_ic_cup(i,:,j) = 0.0
 577          nh4_a_5to6_ic_cup(i,:,j) = 0.0 ; nh4_cw_5to6_ic_cup(i,:,j) = 0.0
 578          no3_a_5to6_ic_cup(i,:,j) = 0.0 ; no3_cw_5to6_ic_cup(i,:,j) = 0.0
 579          oa_a_5to6_ic_cup (i,:,j) = 0.0 ;  oa_cw_5to6_ic_cup(i,:,j) = 0.0
 580          oin_a_5to6_ic_cup (i,:,j) = 0.0 ;  oin_cw_5to6_ic_cup(i,:,j) = 0.0
 581          bc_a_5to6_ic_cup (i,:,j) = 0.0 ;  bc_cw_5to6_ic_cup(i,:,j) = 0.0
 582          na_a_5to6_ic_cup (i,:,j) = 0.0 ;  na_cw_5to6_ic_cup(i,:,j) = 0.0
 583          cl_a_5to6_ic_cup (i,:,j) = 0.0 ;  cl_cw_5to6_ic_cup(i,:,j) = 0.0
 584          water_5to6_ic_cup (i,:,j) = 0.0
 585
 586       end if
 587
 588       if (icalcflagaa <= 0) cycle main_i_loop
 589
 590
 591       write(*,'(/a,i10,4i5)') &
 592          'mosaic_chem_cup_driver doing ktau, id, i, j, ishall =', ktau, grid_id, i, j, ishall
 593
 594       chem_cupflag_1d(kts:kte) = chem_cupflag(i,kts:kte,j)
 595
 596       chem_1d(kts:kte,1:num_chem) = chem(i,kts:kte,j,1:num_chem)
 597
 598       cldfra_ls_1d(kts:kte)    = cldfra(i,kts:kte,j)
 599       dz_1d(kts:kte)           = dz8w(i,kts:kte,j)
 600       pcen_1d(kts:kte)         = p_phy(i,kts:kte,j)
 601       ph_no2_1d(kts:kte)       = ph_no2(i,kts:kte,j)
 602       rhocen_1d(kts:kte)       = 1.0_r8/alt(i,kts:kte,j)
 603       tcen_1d(kts:kte)         = t_phy(i,kts:kte,j)
 604       zbnd_1d(kts:ktep1)          = z_at_w(i,kts:ktep1,j)
 605       zcen_1d(kts:kte)         = zmid(i,kts:kte,j)
 606
 607       qc_ic_cup_1d(kts:kte)    = qc_ic_cup(i,kts:kte,j)
 608       qc_iu_cup_1d(kts:kte)    = qc_iu_cup(i,kts:kte,j)
 609       qi_ic_cup_1d(kts:kte)    = 0.0_r8
 610       qi_iu_cup_1d(kts:kte)    = 0.0_r8
 611       wact_cup_1d                = wact_cup(i,j)
 612
 613       fcvt_qc_to_pr_cup_1d(kts:kte) = fcvt_qc_to_pr_cup(i,kts:kte,j)
 614       fcvt_qc_to_qi_cup_1d(kts:kte) = fcvt_qc_to_qi_cup(i,kts:kte,j)
 615       fcvt_qi_to_pr_cup_1d(kts:kte) = fcvt_qi_to_pr_cup(i,kts:kte,j)
 616
 617       if ( mode_chemcup_timeinteg == 1 ) then
 618          ! simulate the entire cumulus cloud fraction using substeps of dt = dtstepc
 619          !    also force dtsub <= nca+dtstep because tcloud_cup may not be an
 620          !       integer multiple of dtstepc
 621          !    use nca+dtstep because nca is decremented in phys before chem sees it
 622          tcloud_cup_1d = min( nca(i,j)+dtstep, dtstepc )
 623       else if ( mode_chemcup_timeinteg == 2 ) then
 624          ! simulate the entire cumulus cloud fraction in a "one time only" calculation
 625          !    with dt = tcloud_cup
 626          tcloud_cup_1d = tcloud_cup(i,j)
 627       else
 628          call wrf_error_fatal( &
 629             '*** mosaic_chem_cup_driver -- bad value for mode_chemcup_timeinteg' )
 630       end if
 631       cldfra_cup_1d(kts:kte)   = cldfra_cup(i,kts:kte,j)
 632       updfra_cup_1d(kts:kte)   = updfra_cup(i,kts:kte,j)
 633       mfup_cup_1d(kts:ktep1)      = mfup_cup(i,kts:ktep1,j)
 634       mfup_ent_cup_1d(kts:kte) = mfup_ent_cup(i,kts:kte,j)
 635       mfdn_cup_1d(kts:ktep1)     = mfdn_cup(i,kts:ktep1,j)
 636       mfdn_ent_cup_1d(kts:kte) = mfdn_ent_cup(i,kts:kte,j)
 637
 638       if (idiagee > 0) write(*,'(a,i6,l5,i5,20x,2f10.1)') &
 639          'chcup_a20 ktau, cupflag, ishall,      tcloud, tcloud1d', &
 640          ktau, cupflag(i,j), nint(shall(i,j)), tcloud_cup(i,j), tcloud_cup_1d
 641
 642 ! *** following are for testing
 643       if (cldfra_ls_testvalue >= 0.0) &
 644          cldfra_ls_1d(kts:kte) = max( 0.0, min( 1.0, cldfra_ls_testvalue ) )
 645       if (wact_cup_testvalue >= 0.0) &
 646          wact_cup_1d = wact_cup_testvalue
 647
 648
 649       kcldbot_1d = kts-1
 650       kcldtop_1d = kts-1
 651       do k = kts, kte
 652          if (cldfra_cup_1d(k) > 0.0_r8) then
 653             kcldtop_1d = k
 654             if (kcldbot_1d < kts) kcldbot_1d = k
 655          end if
 656       end do
 657
 658
 659 !     subr. chem_cup_1d( &
 660 !        config_flags, aer_mech_id, &
 661 !        lundiag, lunerr, &
 662 !        kts, kte, ktep1, p1st, num_chem, num_moist, &
 663 !        ktau, grid_id, i, j, &
 664 !        ishall, kcldbot_inp, kcldtop_inp, &
 665 !        tau_active, tau_inactive, &
 666 !        dz, zcen, zbnd, pcen, tcen, rhocen, ph_no2, &
 667 !        af_lscld, af_cucld_inp, af_up_inp, &
 668 !        qcw_incu_inp, qci_incu_inp, &
 669 !        qcw_inup_inp, qci_inup_inp, &
 670 !        mf_up_inp, mf_up_ent_inp, &
 671 !        mf_dn_inp, mf_dn_ent_inp, &
 672 !        fcvt_qc_to_pr_inp, fcvt_qc_to_qi_inp, fcvt_qi_to_pr_inp, &
 673 !        wact_inp, &
 674 !        chem, chem_incu, chem_name, chem_cupflag, &
 675 !        maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
 676 !        ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
 677 !        ai_phase, cw_phase, msectional, &
 678 !        massptr_aer, numptr_aer, &
 679 !        lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
 680 !        dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer
 681
 682       call chem_cup_1d( &
 683          config_flags, aer_mech_id, &
 684          lunerr, lundiag, idiagee, &
 685          kts, kte, ktep1, param_first_scalar, num_chem, num_moist, &
 686          ktau, grid_id, i, j, &
 687          ishall, kcldbot_1d, kcldtop_1d, &
 688          tcloud_cup_1d, tcloud_cup_1d, &
 689          dz_1d, zcen_1d, zbnd_1d, pcen_1d, tcen_1d, rhocen_1d, ph_no2_1d, &
 690          cldfra_ls_1d, cldfra_cup_1d, updfra_cup_1d, &
 691          qc_ic_cup_1d, qi_ic_cup_1d, &
 692          qc_iu_cup_1d, qi_iu_cup_1d, &
 693          mfup_cup_1d, mfup_ent_cup_1d, &
 694          mfdn_cup_1d, mfdn_ent_cup_1d, &
 695          fcvt_qc_to_pr_cup_1d, fcvt_qc_to_qi_cup_1d, fcvt_qi_to_pr_cup_1d, &
 696          wact_cup_1d, &
 697          chem_1d, chem_incu, chem_name, chem_cupflag_1d, &
 698          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
 699          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
 700          ai_phase, cw_phase, msectional, &
 701          massptr_aer, numptr_aer, &
 702          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
 703          dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer )
 704
 705       write(*,'(/a,4i10)') &
 706          'mosaic_chem_cup_driver back  ktau, id, i, j =', ktau, grid_id, i, j
 707
 708       chem_cupflag(i,kts:kte,j) = chem_cupflag_1d(kts:kte)
 709
 710       chem(i,kts:kte,j,1:num_chem) = chem_1d(kts:kte,1:num_chem)
 711
 712 !     call wrf_error_fatal( 'aborting after first call to chem_cup_1d' )
 713 !     if (ktau >= 465) call wrf_error_fatal( &
 714 !        'aborting after call to chem_cup_1d and ktau>=465' )
 715
 716       if (mode_chemcup_timeinteg == 2) then
 717          ! calculations for this cloud are totally finished
 718          ! set tcloud_cup to a negative value to indicate this on next timestep
 719          tcloud_cup(i,j) = -tcloud_cup(i,j)
 720       end if
 721
 722 !  Do the gas variables co and hno3 here
 723          do k = kcldbot_1d, kcldtop_1d
 724             if (chem_cupflag(i,k,j) <= 0) cycle
 725             co_a_ic_cup(i,k,j)  = co_a_ic_cup(i,k,j) &
 726                                      + chem_incu(k,p_co)
 727             hno3_a_ic_cup(i,k,j)  = hno3_a_ic_cup(i,k,j) &
 728                                      + chem_incu(k,p_hno3)
 729          enddo
 730
 731       do itype = 1, ntype_aer
 732       do isize = 1, nsize_aer(itype)
 733          if (dcen_sect(isize,itype) .le. 0.625e-4) then
 734          do k = kcldbot_1d, kcldtop_1d
 735             if (chem_cupflag(i,k,j) <= 0) cycle
 736             so4_a_1to4_ic_cup(i,k,j)  = so4_a_1to4_ic_cup(i,k,j) &
 737                                      + chem_incu(k,lptr_so4_aer(isize,itype,ai_phase))
 738             so4_cw_1to4_ic_cup(i,k,j) = so4_cw_1to4_ic_cup(i,k,j) &
 739                                      + chem_incu(k,lptr_so4_aer(isize,itype,cw_phase))
 740             nh4_a_1to4_ic_cup(i,k,j)  = nh4_a_1to4_ic_cup(i,k,j) &
 741                                      + chem_incu(k,lptr_nh4_aer(isize,itype,ai_phase))
 742             nh4_cw_1to4_ic_cup(i,k,j) = nh4_cw_1to4_ic_cup(i,k,j) &
 743                                      + chem_incu(k,lptr_nh4_aer(isize,itype,cw_phase))
 744             no3_a_1to4_ic_cup(i,k,j)  = no3_a_1to4_ic_cup(i,k,j) &
 745                                      + chem_incu(k,lptr_no3_aer(isize,itype,ai_phase))
 746             no3_cw_1to4_ic_cup(i,k,j) = no3_cw_1to4_ic_cup(i,k,j) &
 747                                      + chem_incu(k,lptr_no3_aer(isize,itype,cw_phase))
 748             oa_a_1to4_ic_cup(i,k,j)   = oa_a_1to4_ic_cup(i,k,j) &
 749                                      + chem_incu(k,lptr_oc_aer(isize,itype,ai_phase)) &
 750                                      + chem_incu(k,lptr_pcg1_b_c_aer(isize,itype,ai_phase)) &
 751                                      + chem_incu(k,lptr_pcg1_b_o_aer(isize,itype,ai_phase)) &
 752                                      + chem_incu(k,lptr_opcg1_b_c_aer(isize,itype,ai_phase)) &
 753                                      + chem_incu(k,lptr_opcg1_b_o_aer(isize,itype,ai_phase)) &
 754                                      + chem_incu(k,lptr_pcg1_f_c_aer(isize,itype,ai_phase)) &
 755                                      + chem_incu(k,lptr_pcg1_f_o_aer(isize,itype,ai_phase)) &
 756                                      + chem_incu(k,lptr_opcg1_f_c_aer(isize,itype,ai_phase)) &
 757                                      + chem_incu(k,lptr_opcg1_f_o_aer(isize,itype,ai_phase)) &
 758                                      + chem_incu(k,lptr_ant1_c_aer(isize,itype,ai_phase)) &
 759                                      + chem_incu(k,lptr_biog1_c_aer(isize,itype,ai_phase))
 760
 761
 762             oa_cw_1to4_ic_cup(i,k,j)  = oa_cw_1to4_ic_cup(i,k,j) &
 763                                      + chem_incu(k,lptr_oc_aer(isize,itype,cw_phase)) &
 764                                      + chem_incu(k,lptr_pcg1_b_c_aer(isize,itype,cw_phase)) &
 765                                      + chem_incu(k,lptr_pcg1_b_o_aer(isize,itype,cw_phase)) &
 766                                      + chem_incu(k,lptr_opcg1_b_c_aer(isize,itype,cw_phase)) &
 767                                      + chem_incu(k,lptr_opcg1_b_o_aer(isize,itype,cw_phase)) &
 768                                      + chem_incu(k,lptr_pcg1_f_c_aer(isize,itype,cw_phase)) &
 769                                      + chem_incu(k,lptr_pcg1_f_o_aer(isize,itype,cw_phase)) &
 770                                      + chem_incu(k,lptr_opcg1_f_c_aer(isize,itype,cw_phase)) &
 771                                      + chem_incu(k,lptr_opcg1_f_o_aer(isize,itype,cw_phase)) &
 772                                      + chem_incu(k,lptr_ant1_c_aer(isize,itype,cw_phase)) &
 773                                      + chem_incu(k,lptr_biog1_c_aer(isize,itype,cw_phase))
 774
 775             oin_a_1to4_ic_cup(i,k,j)  = oin_a_1to4_ic_cup(i,k,j) &
 776                                      + chem_incu(k,lptr_oin_aer(isize,itype,ai_phase))
 777             oin_cw_1to4_ic_cup(i,k,j) = oin_cw_1to4_ic_cup(i,k,j) &
 778                                      + chem_incu(k,lptr_oin_aer(isize,itype,cw_phase))
 779
 780             bc_a_1to4_ic_cup(i,k,j)  = bc_a_1to4_ic_cup(i,k,j) &
 781                                      + chem_incu(k,lptr_bc_aer(isize,itype,ai_phase))
 782             bc_cw_1to4_ic_cup(i,k,j) = bc_cw_1to4_ic_cup(i,k,j) &
 783                                      + chem_incu(k,lptr_bc_aer(isize,itype,cw_phase))
 784
 785             na_a_1to4_ic_cup(i,k,j)  = na_a_1to4_ic_cup(i,k,j) &
 786                                      + chem_incu(k,lptr_na_aer(isize,itype,ai_phase))
 787             na_cw_1to4_ic_cup(i,k,j) = na_cw_1to4_ic_cup(i,k,j) &
 788                                      + chem_incu(k,lptr_na_aer(isize,itype,cw_phase))
 789
 790             cl_a_1to4_ic_cup(i,k,j)  = cl_a_1to4_ic_cup(i,k,j) &
 791                                      + chem_incu(k,lptr_cl_aer(isize,itype,ai_phase))
 792             cl_cw_1to4_ic_cup(i,k,j) = cl_cw_1to4_ic_cup(i,k,j) &
 793                                      + chem_incu(k,lptr_cl_aer(isize,itype,cw_phase))
 794
 795             water_1to4_ic_cup(i,k,j)  = water_1to4_ic_cup(i,k,j) &
 796                                      + chem_incu(k,waterptr_aer(isize,itype))
 797
 798          end do ! k
 799
 800            elseif (dcen_sect(isize,itype) .gt. 0.625e-4 .and. &
 801                      dcen_sect(isize,itype) .le. 2.5e-4) then
 802          do k = kcldbot_1d, kcldtop_1d
 803             if (chem_cupflag(i,k,j) <= 0) cycle
 804             so4_a_5to6_ic_cup(i,k,j)  = so4_a_5to6_ic_cup(i,k,j) &
 805                                      + chem_incu(k,lptr_so4_aer(isize,itype,ai_phase))
 806             so4_cw_5to6_ic_cup(i,k,j) = so4_cw_5to6_ic_cup(i,k,j) &
 807                                      + chem_incu(k,lptr_so4_aer(isize,itype,cw_phase))
 808             nh4_a_5to6_ic_cup(i,k,j)  = nh4_a_5to6_ic_cup(i,k,j) &
 809                                      + chem_incu(k,lptr_nh4_aer(isize,itype,ai_phase))
 810             nh4_cw_5to6_ic_cup(i,k,j) = nh4_cw_5to6_ic_cup(i,k,j) &
 811                                      + chem_incu(k,lptr_nh4_aer(isize,itype,cw_phase))
 812             no3_a_5to6_ic_cup(i,k,j)  = no3_a_5to6_ic_cup(i,k,j) &
 813                                      + chem_incu(k,lptr_no3_aer(isize,itype,ai_phase))
 814             no3_cw_5to6_ic_cup(i,k,j) = no3_cw_5to6_ic_cup(i,k,j) &
 815                                      + chem_incu(k,lptr_no3_aer(isize,itype,cw_phase))
 816             oa_a_5to6_ic_cup(i,k,j)   = oa_a_5to6_ic_cup(i,k,j) &
 817                                      + chem_incu(k,lptr_oc_aer(isize,itype,ai_phase)) &
 818                                      + chem_incu(k,lptr_pcg1_b_c_aer(isize,itype,ai_phase)) &
 819                                      + chem_incu(k,lptr_pcg1_b_o_aer(isize,itype,ai_phase)) &
 820                                      + chem_incu(k,lptr_opcg1_b_c_aer(isize,itype,ai_phase)) &
 821                                      + chem_incu(k,lptr_opcg1_b_o_aer(isize,itype,ai_phase)) &
 822                                      + chem_incu(k,lptr_pcg1_f_c_aer(isize,itype,ai_phase)) &
 823                                      + chem_incu(k,lptr_pcg1_f_o_aer(isize,itype,ai_phase)) &
 824                                      + chem_incu(k,lptr_opcg1_f_c_aer(isize,itype,ai_phase)) &
 825                                      + chem_incu(k,lptr_opcg1_f_o_aer(isize,itype,ai_phase)) &
 826                                      + chem_incu(k,lptr_ant1_c_aer(isize,itype,ai_phase)) &
 827                                      + chem_incu(k,lptr_biog1_c_aer(isize,itype,ai_phase))
 828
 829
 830             oa_cw_5to6_ic_cup(i,k,j)  = oa_cw_5to6_ic_cup(i,k,j) &
 831                                      + chem_incu(k,lptr_oc_aer(isize,itype,cw_phase)) &
 832                                      + chem_incu(k,lptr_pcg1_b_c_aer(isize,itype,cw_phase)) &
 833                                      + chem_incu(k,lptr_pcg1_b_o_aer(isize,itype,cw_phase)) &
 834                                      + chem_incu(k,lptr_opcg1_b_c_aer(isize,itype,cw_phase)) &
 835                                      + chem_incu(k,lptr_opcg1_b_o_aer(isize,itype,cw_phase)) &
 836                                      + chem_incu(k,lptr_pcg1_f_c_aer(isize,itype,cw_phase)) &
 837                                      + chem_incu(k,lptr_pcg1_f_o_aer(isize,itype,cw_phase)) &
 838                                      + chem_incu(k,lptr_opcg1_f_c_aer(isize,itype,cw_phase)) &
 839                                      + chem_incu(k,lptr_opcg1_f_o_aer(isize,itype,cw_phase)) &
 840                                      + chem_incu(k,lptr_ant1_c_aer(isize,itype,cw_phase)) &
 841                                      + chem_incu(k,lptr_biog1_c_aer(isize,itype,cw_phase))
 842
 843             oin_a_5to6_ic_cup(i,k,j)  = oin_a_5to6_ic_cup(i,k,j) &
 844                                      + chem_incu(k,lptr_oin_aer(isize,itype,ai_phase))
 845             oin_cw_5to6_ic_cup(i,k,j) = oin_cw_5to6_ic_cup(i,k,j) &
 846                                      + chem_incu(k,lptr_oin_aer(isize,itype,cw_phase))
 847
 848             bc_a_5to6_ic_cup(i,k,j)  = bc_a_5to6_ic_cup(i,k,j) &
 849                                      + chem_incu(k,lptr_bc_aer(isize,itype,ai_phase))
 850             bc_cw_5to6_ic_cup(i,k,j) = bc_cw_5to6_ic_cup(i,k,j) &
 851                                      + chem_incu(k,lptr_bc_aer(isize,itype,cw_phase))
 852
 853             na_a_5to6_ic_cup(i,k,j)  = na_a_5to6_ic_cup(i,k,j) &
 854                                      + chem_incu(k,lptr_na_aer(isize,itype,ai_phase))
 855             na_cw_5to6_ic_cup(i,k,j) = na_cw_5to6_ic_cup(i,k,j) &
 856                                      + chem_incu(k,lptr_na_aer(isize,itype,cw_phase))
 857
 858             cl_a_5to6_ic_cup(i,k,j)  = cl_a_5to6_ic_cup(i,k,j) &
 859                                      + chem_incu(k,lptr_cl_aer(isize,itype,ai_phase))
 860             cl_cw_5to6_ic_cup(i,k,j) = cl_cw_5to6_ic_cup(i,k,j) &
 861                                      + chem_incu(k,lptr_cl_aer(isize,itype,cw_phase))
 862
 863             water_5to6_ic_cup(i,k,j)  = water_5to6_ic_cup(i,k,j) &
 864                                      + chem_incu(k,waterptr_aer(isize,itype))
 865
 866
 867          end do ! k
 868            end if ! dcen size
 869
 870       end do ! isize
 871       end do ! itype
 872
 873
 874       end do main_i_loop
 875       end do main_j_loop
 876
 877
 878       return
 879       end subroutine mosaic_chem_cup_driver
 880
 881
 882 !-------------------------------------------------------------------------------
 883       subroutine chem_cup_1d( &
 884          config_flags, aer_mech_id, &
 885          lunerr, lundiag, idiagaa_inp, &
 886          kts, kte, ktep1, p1st, num_chem, num_moist, &
 887          ktau, grid_id, i, j, &
 888          ishall, kcldbot_inp, kcldtop_inp, &
 889          tau_active, tau_inactive, &
 890          dz, zcen, zbnd, pcen, tcen, rhocen, ph_no2, &
 891          af_lscld, af_cucld_inp, af_up_inp, &
 892          qcw_incu_inp, qci_incu_inp, &
 893          qcw_inup_inp, qci_inup_inp, &
 894          mf_up_inp, mf_up_ent_inp, &
 895          mf_dn_inp, mf_dn_ent_inp, &
 896          fcvt_qc_to_pr, fcvt_qc_to_qi, fcvt_qi_to_pr, &
 897          wact_inp, &
 898          chem, chem_incu, chem_name, chem_cupflag, &
 899          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
 900          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
 901          ai_phase, cw_phase, msectional, &
 902          massptr_aer, numptr_aer, &
 903          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
 904          dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer )
 905
 906       use module_configure, only:  grid_config_rec_type
 907       use module_configure, only:  &
 908          p_qc, &
 909          p_h2o2, p_hcl, p_hno3, p_nh3, p_so2, p_sulf
 910
 911 ! arguments
 912 !
 913 ! note: arguments ending in "_inp" may be adjusted for consistency
 914       type(grid_config_rec_type), intent(in) :: config_flags
 915
 916       integer, intent(in) :: aer_mech_id
 917       integer, intent(in) :: grid_id
 918       integer, intent(in) :: i, idiagaa_inp, ishall
 919       integer, intent(in) :: j
 920       integer, intent(in) :: ktau, kts, kte, ktep1
 921       integer, intent(in) :: kcldbot_inp, kcldtop_inp
 922       integer, intent(in) :: lunerr, lundiag
 923       integer, intent(in) :: num_chem, num_moist
 924       integer, intent(in) :: p1st
 925
 926       integer, intent(inout) :: chem_cupflag(kts:kte)  ! the "in" value is 0
 927
 928       integer, intent(in) :: maxd_acomp, maxd_aphase, maxd_atype, maxd_asize
 929       integer, intent(in) :: &
 930          nphase_aer, ntype_aer, &
 931          ai_phase, cw_phase, msectional, &
 932          nsize_aer( maxd_atype ), &
 933          ncomp_aer( maxd_atype ), &
 934          massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
 935          numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
 936          lptr_so4_aer(maxd_asize, maxd_atype, maxd_aphase), &
 937          lptr_oin_aer(maxd_asize, maxd_atype, maxd_aphase), &
 938          lptr_no3_aer(maxd_asize, maxd_atype, maxd_aphase), &
 939          lptr_nh4_aer(maxd_asize, maxd_atype, maxd_aphase), &
 940          lptr_cl_aer(maxd_asize, maxd_atype, maxd_aphase)
 941       real(r4), intent(in) :: &
 942          dlo_sect( maxd_asize, maxd_atype ), &
 943          dhi_sect( maxd_asize, maxd_atype ), &
 944          dens_aer( maxd_acomp, maxd_atype ), &
 945          hygro_aer( maxd_acomp, maxd_atype ), &
 946          sigmag_aer(maxd_asize, maxd_atype)
 947
 948
 949       real(r8), intent(in) :: tau_active, tau_inactive
 950       real(r8), intent(in) :: wact_inp  ! vertical velocity for activation at cloud-base (m/s)
 951
 952       real(r8), intent(in), dimension(kts:kte) :: &
 953          dz, &            ! layer thickness (m)
 954          zcen, &          ! height at layer center (m)
 955          pcen, &          ! pressure at layer center (pa)
 956          tcen, &          ! temperature at layer center (K)
 957          rhocen, &        ! dry air density at layer center (kg/m3)
 958          ph_no2, &        ! no2 photolysis rate (1/s)
 959          af_lscld, &      ! grid-resolved cloud fractional area
 960          af_cucld_inp, &  ! cumulus cloud fractional area
 961          af_up_inp, &     ! updraft fractional area
 962          qcw_incu_inp, &  ! cloud-water mixing ratio in cumulus cloud (kg/kg)
 963          qci_incu_inp, &  ! cloud-ice   mixing ratio in updraft       (kg/kg)
 964          qcw_inup_inp, &  ! cloud-water mixing ratio in cumulus cloud (kg/kg)
 965          qci_inup_inp, &  ! cloud-ice   mixing ratio in updraft       (kg/kg)
 966          mf_up_ent_inp, & ! change to updraft mass flux in layer from entrainment (kg/m2/s)
 967          mf_dn_ent_inp, & ! change to updraft mass flux in layer from entrainment (kg/m2/s)
 968          fcvt_qc_to_pr, & ! fraction of cloud-water converted to precip    as air move thru updraft layer (-)
 969          fcvt_qc_to_qi, & ! fraction of cloud-water converted to cloud-ice as air move thru updraft layer (-)
 970          fcvt_qi_to_pr    ! fraction of cloud-ice   converted to precip    as air move thru updraft layer (-)
 971
 972       real(r8), intent(in), dimension(kts:ktep1) :: &
 973          zbnd, &          ! height at layer boundary (m)
 974          mf_up_inp, &     ! updraft mass flux at layer boundary (kg/m2/s)
 975          mf_dn_inp        ! updraft mass flux at layer boundary (kg/m2/s)
 976
 977       real(r8), intent(inout), dimension(kts:kte,num_chem) :: &
 978          chem, &          ! grid-average aerosol and trace-gas mixing ratios (ug/kg, #/kg, or ppmv)
 979          chem_incu        ! in-convective-cloud aerosol and trace-gas mixing ratios (ug/kg, #/kg, or ppmv)
 980
 981       character(len=12), intent(in) :: chem_name(num_chem)
 982
 983 ! local variables
 984       integer :: aip, cwp, typ
 985       integer :: idiagaa, idiagbb, iflagaa, iok, icomp, isize, itype, itmpa, itmpb
 986       integer :: ido_inact(kts:kte)
 987       integer :: jtsub
 988       integer :: k, kaa, kzz
 989       integer :: kcldbot, kcldtop, kcldbotliq
 990       integer :: kdiagbot, kdiagtop
 991       integer :: kupdrbot, kupdrtop, kdndrbot, kdndrtop
 992       integer :: l, la, lc, l2, l3, lundiagbb
 993       integer :: m
 994       integer :: n
 995       integer :: ntsub  ! number of time sub-steps for active cloud integration
 996
 997       logical, parameter :: do_activa  = .true.
 998
 999 !------To turn off cloud chemistry set do_aqchem-false : Manish Shrivastava -------------------
1000       logical, parameter :: do_aqchem  = .true.
1001 !---------------------------------------------------------------------------
1002       logical, parameter :: do_inact   = .true.
1003       logical, parameter :: do_resusp  = .true.
1004       logical, parameter :: do_updraft = .true.
1005
1006       logical, parameter :: do_2ndact_deep  = .true.
1007       logical, parameter :: do_2ndact_shal  = .false.
1008       logical, parameter :: do_dndraft_deep = .true.
1009       logical, parameter :: do_dndraft_shal = .false.
1010       logical, parameter :: do_wetrem_deep  = .true.
1011       logical, parameter :: do_wetrem_shal  = .false.
1012
1013       logical :: do_2ndact, do_dndraft, do_wetrem
1014
1015       real(r8), parameter :: rerrtol1_mbal = 3.0e-6
1016
1017       real(r8) :: af_cucld(kts:kte)  ! cumulus cloud fractional area
1018       real(r8) :: af_dn(kts:kte)     ! downdraft fractional area
1019       real(r8) :: af_ev(kts:kte)     ! environment fractional area
1020       real(r8) :: af_up(kts:kte)     ! updraft fractional area
1021       real(r8) :: af_inact(kts:kte)  ! inactive cumulus cloud fractional area
1022       real(r8) :: chem_av_new(kts:kte,num_chem)      ! grid-average chem values at end   of sub time-step
1023       real(r8) :: chem_av_old(kts:kte,num_chem)      ! grid-average chem values at start of sub time-step
1024       real(r8) :: chem_inact_dsp(kts:kte,num_chem)   ! inactive cloud chem values at end   of time-step
1025                                                      ! but after  resuspension/dissipation calculations
1026       real(r8) :: chem_inact_new(kts:kte,num_chem)   ! inactive cloud chem values at end   of time-step
1027                                                      ! but before resuspension/dissipation calculations
1028       real(r8) :: chem_inact_old(kts:kte,num_chem)   ! inactive cloud chem values at start of time-step
1029       real(r8) :: chem_dn(kts:ktep1,num_chem)        ! steady-state dndraft chem values (at layer boundaries)
1030       real(r8) :: chem_up(kts:ktep1,num_chem)        ! steady-state updraft chem values (at layer boundaries)
1031       ! the following dchemdt_xx_yyyyyy arrays are d(chem)/dt in environment or updraft
1032       ! due to resuspension, activation, aqueous-chemistry, or wet-removal [(ug/kg or ...)/s]
1033       real(r8) :: dchemdt_ev_resusp(kts:kte,num_chem)
1034       real(r8) :: dchemdt_up_activa(kts:kte,num_chem)
1035       real(r8) :: dchemdt_up_aqchem(kts:kte,num_chem)
1036       real(r8) :: dchemdt_up_wetrem(kts:kte,num_chem)
1037       ! the following del_chem_yyyyyy arrays are change (delta)
1038       ! of grid-average chem due to various processes (ug/kg or ...)
1039       real(r8) :: del_chem_activa(kts:kte,num_chem)  ! activation
1040       real(r8) :: del_chem_aqchem(kts:kte,num_chem)  ! aqueous-chemistry
1041       real(r8) :: del_chem_residu(kts:kte,num_chem)  ! residual
1042       real(r8) :: del_chem_resusp(kts:kte,num_chem)  ! resuspension
1043       real(r8) :: del_chem_totall(kts:kte,num_chem)  ! total change
1044       real(r8) :: del_chem_wetrem(kts:kte,num_chem)  ! wet-removal
1045       real(r8) :: del_chem_ztrans(kts:kte,num_chem)  ! vertical-transport
1046       real(r8) :: del_chem_actvbb(kts:kte,num_chem)  ! activation for inactive cloud
1047       real(r8) :: del_chem_aqchbb(kts:kte,num_chem)  ! aqueous-chemistry for inactive cloud
1048       real(r8) :: del_chem_resdbb(kts:kte,num_chem)  ! residual for inactive cloud
1049       real(r8) :: del_chem_respbb(kts:kte,num_chem)  ! resuspension for inactive cloud
1050       real(r8) :: del_chem_totlbb(kts:kte,num_chem)  ! total change for inactive cloud
1051       real(r8) :: dtsub                 ! time sub-step for active cloud integration (s)
1052       real(r8) :: mf_dn_det(kts:kte)    ! change to downdraft mass flux in layer from detrainment (kg/m2/s)
1053       real(r8) :: mf_dn_ent(kts:kte)    ! change to downdraft mass flux in layer from entrainment (kg/m2/s)
1054       real(r8) :: mf_dn(kts:ktep1)      ! downdraft mass flux at layer boundary (kg/m2/s)
1055       real(r8) :: mf_ev(kts:ktep1)      ! environment mass flux at layer boundary (kg/m2/s)
1056       real(r8) :: mf_up_det(kts:kte)    ! change to updraft mass flux in layer from detrainment (kg/m2/s)
1057       real(r8) :: mf_up_ent(kts:kte)    ! adjusted change to updraft mass flux in layer from entrainment (kg/m2/s)
1058       real(r8) :: mf_up(kts:ktep1)      ! adjusted updraft mass flux at layer boundary (kg/m2/s)
1059       real(r8) :: qci_incu(kts:kte), qci_inup(kts:kte)
1060       real(r8) :: qcw_incu(kts:kte), qcw_inup(kts:kte)
1061       real(r8), parameter :: qcw_cldchem_cutoff = 1.0e-6_r8
1062       real(r8) :: rhodz(kts:kte)        ! rho*dz (kg/m2)
1063       real(r8) :: rhodzsum              ! sum( rho*dz )
1064       real(r8) :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg
1065       real(r8) :: tmp_chem_dn(num_chem)     ! temporary/working value of chem_dn
1066       real(r8) :: tmp_chem_up(num_chem)     ! temporary/working value of chem_up
1067       real(r8) :: tmp_dt                    ! temporary/working value of dt for cloud chemistry
1068       real(r8) :: tmp_fmact(maxd_asize,maxd_atype), tmp_fnact(maxd_asize,maxd_atype)
1069       real(r8) :: tmp_gas_aqfrac_up(num_chem)  ! temporary/working value of gas_aqfrac in updraft
1070       real(r8) :: tmp_mfxchem_dn(num_chem)  ! temporary/working value of mf_up*chem_dn
1071       real(r8) :: tmp_mfxchem_up(num_chem)  ! temporary/working value of mf_up*chem_up
1072       real(r8) :: tmp_mf_dn ! temporary/working value of mf_dn
1073       real(r8) :: tmp_mf_up ! temporary/working value of mf_up
1074       real(r8) :: tmp_w_up  ! temporary/working value of w in updraft
1075       real(r8) :: tmpveca(201), tmpvecb(201)
1076       ! tmp_zflux_... - temporary/working values of trace-species vertical fluxes [(ug/kg or ...)/m2/s]
1077       real(r8) :: tmp_zflux_bot, tmp_zflux_botev, tmp_zflux_botup, tmp_zflux_botdn
1078       real(r8) :: tmp_zflux_top, tmp_zflux_topev, tmp_zflux_topup, tmp_zflux_topdn
1079       real(r8) :: wact  ! vertical velocity for activation at cloud-base (m/s)
1080       real(r8), parameter :: wact_min = 0.2  ! (m/s)
1081       ! the following zav_yyyyyy arrays are vertical averages of
1082       ! the corresponding yyyyyy arrays, weighted by rhodz
1083       real(r8) :: zav_chem_av_new(num_chem)
1084       real(r8) :: zav_chem_av_old(num_chem)
1085       real(r8) :: zav_chem_dn(num_chem)
1086       real(r8) :: zav_chem_up(num_chem)
1087       real(r8) :: zav_del_chem_activa(num_chem)
1088       real(r8) :: zav_del_chem_aqchem(num_chem)
1089       real(r8) :: zav_del_chem_residu(num_chem)
1090       real(r8) :: zav_del_chem_resusp(num_chem)
1091       real(r8) :: zav_del_chem_totall(num_chem)
1092       real(r8) :: zav_del_chem_wetrem(num_chem)
1093       real(r8) :: zav_del_chem_ztrans(num_chem)
1094       real(r8) :: zav_del_chem_actvbb(num_chem)
1095       real(r8) :: zav_del_chem_aqchbb(num_chem)
1096       real(r8) :: zav_del_chem_resdbb(num_chem)
1097       real(r8) :: zav_del_chem_respbb(num_chem)
1098       real(r8) :: zav_del_chem_totlbb(num_chem)
1099
1100       character(len=160) :: msg
1101
1102
1103 ! sanity checks
1104       if ( (ai_phase < 1) .or. (ai_phase > nphase_aer) .or. &
1105            (cw_phase < 1) .or. (cw_phase > nphase_aer) ) then
1106          write(msg,'(a,3(1x,i10))') &
1107             'chem_cup_1d - bad ai_phase, cw_phase, nphase_aer =', &
1108             ai_phase, cw_phase, nphase_aer
1109          call wrf_message( msg )
1110          call wrf_error_fatal( msg )
1111       end if
1112       if (aer_mech_id /= 3) then
1113          write(msg,'(a,3(1x,i10))') &
1114             'chem_cup_1d - bad aer_mech_id = ', aer_mech_id
1115          call wrf_message( msg )
1116          call wrf_error_fatal( msg )
1117       end if
1118
1119 ! checks on inputs
1120       idiagaa = 0
1121       if (lundiag > 0) idiagaa = idiagaa_inp
1122       if (idiagaa > 0) write(lundiag,'(//a,i10,4i5)') &
1123          'chem_cup_1d doing ktau, id, i, j, ishall =', ktau, grid_id, i, j, ishall
1124
1125       idiagbb = 0
1126       if (idiagaa_inp >      0) idiagbb = idiagaa_inp
1127       if (idiagaa_inp <= -1000) idiagbb = -idiagaa_inp/1000
1128       lundiagbb = 6
1129       if (lundiag > 0) lundiagbb = lundiag
1130
1131       if (ishall == 1) then
1132          do_2ndact  = do_2ndact_shal
1133          do_dndraft = do_dndraft_shal
1134          do_wetrem  = do_wetrem_shal
1135       else
1136          do_2ndact  = do_2ndact_deep
1137          do_dndraft = do_dndraft_deep
1138          do_wetrem  = do_wetrem_deep
1139       end if
1140
1141       kcldbot = kcldbot_inp
1142       kcldtop = kcldtop_inp
1143
1144       af_cucld(:)  = max( af_cucld_inp(:), 0.0_r8 )
1145       af_up(:)     = max( af_up_inp(:), 0.0_r8 )
1146       qcw_incu(:)  = max( qcw_incu_inp(:), 0.0_r8 )
1147       qci_incu(:)  = max( qci_incu_inp(:), 0.0_r8 )
1148       qcw_inup(:)  = max( qcw_inup_inp(:), 0.0_r8 )
1149       qci_inup(:)  = max( qci_inup_inp(:), 0.0_r8 )
1150       mf_up(:)     = max( mf_up_inp(:), 0.0_r8 )
1151       mf_up_ent(:) = max( mf_up_ent_inp(:), 0.0_r8 )
1152       mf_up_det(:) = 0.0_r8
1153       wact = max( wact_inp, wact_min )
1154
1155       if ( do_dndraft ) then
1156 ! the cloud parameterization does not provide (?) downdraft area
1157 ! the chem_cup only uses af_dn for af_ev = 1 - af_up - af_dn,
1158 !    and it only uses af_ev for determining resuspension fraction of
1159 !    detrained cloud-borne aerosol
1160 ! for now, just leave af_dn = 0
1161          af_dn(:)     = 0.0_r8
1162          mf_dn(:)     = min( mf_dn_inp(:), 0.0_r8 )
1163          mf_dn_ent(:) = max( mf_dn_ent_inp(:), 0.0_r8 )
1164          mf_dn_det(:) = 0.0_r8
1165       else
1166          af_dn(:)     = 0.0_r8
1167          mf_dn(:)     = 0.0_r8
1168          mf_dn_ent(:) = 0.0_r8
1169          mf_dn_det(:) = 0.0_r8
1170       end if
1171
1172       if ( 1 .eq. 1 ) then
1173       if (idiagaa > 0) then
1174       write(lundiag,'(a,2i10)')      'kcldbot/top', kcldbot, kcldtop
1175       write(lundiag,'(a,1p,2e10.2)') 'tau_...    ', tau_active, tau_inactive
1176       write(lundiag,'(a,1p,2e10.2)') 'wact_inp   ', wact_inp
1177       write(lundiag,'(a)') 'zbnd'
1178       write(lundiag,'(1p,15e10.2)') zbnd
1179       write(lundiag,'(a)') 'zcen'
1180       write(lundiag,'(1p,15e10.2)') zcen
1181       write(lundiag,'(a)') 'dz'
1182       write(lundiag,'(1p,15e10.2)') dz
1183       write(lundiag,'(a)') 'pcen'
1184       write(lundiag,'(1p,15e10.2)') pcen
1185       write(lundiag,'(a)') 'tcen'
1186       write(lundiag,'(1p,15e10.2)') tcen
1187       write(lundiag,'(a)') 'rhocen'
1188       write(lundiag,'(1p,15e10.2)') rhocen
1189       write(lundiag,'(a)') 'af_lscld'
1190       write(lundiag,'(1p,15e10.2)') af_lscld
1191       write(lundiag,'(a)') 'af_cucld'
1192       write(lundiag,'(1p,15e10.2)') af_cucld
1193       write(lundiag,'(a)') 'af_up'
1194       write(lundiag,'(1p,15e10.2)') af_up
1195       write(lundiag,'(a)') 'qcw_incu'
1196       write(lundiag,'(1p,15e10.2)') qcw_incu
1197       write(lundiag,'(a)') 'qcw_inup'
1198       write(lundiag,'(1p,15e10.2)') qcw_inup
1199       write(lundiag,'(a)') 'mf_up'
1200       write(lundiag,'(1p,15e10.2)') mf_up
1201       write(lundiag,'(a)') 'mf_up_ent'
1202       write(lundiag,'(1p,15e10.2)') mf_up_ent
1203       if ( do_dndraft ) then
1204       write(lundiag,'(a)') 'mf_dn'
1205       write(lundiag,'(1p,15e10.2)') mf_dn
1206       write(lundiag,'(a)') 'mf_dn_ent'
1207       write(lundiag,'(1p,15e10.2)') mf_dn_ent
1208       end if
1209       end if
1210       end if
1211
1212 !
1213 ! apply "reality checks" to inputs
1214 !
1215       call chem_cup_check_adjust_inputs( &
1216          lunerr, lundiag, idiagaa, &
1217          kts, kte, ktep1, &
1218          ktau, grid_id, i, j, &
1219          ishall, do_dndraft, &
1220          kcldbot, kcldtop, kcldbotliq, &
1221          kupdrbot, kupdrtop, kdndrbot, kdndrtop, &
1222          iok, &
1223          tau_active, tau_inactive, &
1224          dz, zcen, zbnd, pcen, tcen, rhocen, &
1225          af_lscld, af_cucld, af_up, af_dn, &
1226          qcw_incu, qci_incu, &
1227          qcw_inup, qci_inup, &
1228          mf_up, mf_up_ent, mf_up_det, &
1229          mf_dn, mf_dn_ent, mf_dn_det )
1230
1231       if (idiagaa > 0) write(lundiag,'(//a,i10)') &
1232          'chem_cup_check_adjust_inputs iok =', iok
1233
1234       if ( do_dndraft ) then
1235          kdiagbot = max( min(kupdrbot,kdndrbot)-2, kts )
1236          kdiagtop = min( max(kupdrtop,kdndrtop)+2, kte )
1237       else
1238          kdiagbot = max( kupdrbot-2, kts )
1239          kdiagtop = min( kupdrtop+2, kte )
1240       end if
1241 !     kdiagbot = kts
1242 !     kdiagtop = kte
1243
1244       if ( .not. do_updraft ) then
1245          mf_up(:)     = 0.0_r8
1246          mf_up_ent(:) = 0.0_r8
1247          mf_up_det(:) = 0.0_r8
1248       end if
1249
1250 ! diagnostic output
1251       if (idiagaa > 0) then
1252
1253       write(lundiag,'(/4a)') 'k,   ', &
1254             'qcw_incu*1e3 a/b,   qci_incu*1e3 a/b,   ', &
1255             'qcw_inup*1e3 a/b,   qci_inup*1e3 a/b'
1256       do k = kdiagtop, kdiagbot, -1
1257          if (mod(kdiagtop-k,1) == 0) write(lundiag,'(a)')
1258          write(lundiag,'(i2,6(3x,1p,2e10.2))') k, &
1259             qcw_incu_inp(k)*1.0e3, qcw_incu(k)*1.0e3, &
1260             qci_incu_inp(k)*1.0e3, qci_incu(k)*1.0e3, &
1261             qcw_inup_inp(k)*1.0e3, qcw_inup(k)*1.0e3, &
1262             qci_inup_inp(k)*1.0e3, qci_inup(k)*1.0e3
1263       end do
1264
1265       if ( do_dndraft ) then
1266          write(lundiag,'(/a2,2a23,2(2a23,a13))') &
1267             'k', 'af_cucld a/b', 'af_up a/b', &
1268             'mf_up a/b', 'mf_up_ent a/b', 'mf_up_det b', &
1269             'mf_dn a/b', 'mf_dn_ent a/b', 'mf_dn_det b'
1270       else
1271          write(lundiag,'(/a2,2a23,2(2a23,a13))') &
1272             'k', 'af_cucld a/b', 'af_up a/b', &
1273             'mf_up a/b', 'mf_up_ent a/b', 'mf_up_det b'
1274       end if
1275       do k = kdiagtop, kdiagbot, -1
1276          if (mod(kdiagtop-k,1) == 0) write(lundiag,'(a)')
1277          if ( do_dndraft ) then
1278             write(lundiag, &
1279                '(i2,1p, 2(3x,2e10.2), 2(2(3x,2e10.2), 3x,e10.2))') k, &
1280                af_cucld_inp(k), af_cucld(k), &
1281                af_up_inp(k), af_up(k), &
1282                mf_up_inp(k), mf_up(k), &
1283                mf_up_ent_inp(k), mf_up_ent(k), &
1284                mf_up_det(k), &
1285                mf_dn_inp(k), mf_dn(k), &
1286                mf_dn_ent_inp(k), mf_dn_ent(k), &
1287                mf_dn_det(k)
1288          else
1289             write(lundiag, &
1290                '(i2,1p, 2(3x,2e10.2), 2(2(3x,2e10.2), 3x,e10.2))') k, &
1291                af_cucld_inp(k), af_cucld(k), &
1292                af_up_inp(k), af_up(k), &
1293                mf_up_inp(k), mf_up(k), &
1294                mf_up_ent_inp(k), mf_up_ent(k), &
1295                mf_up_det(k)
1296          end if
1297       end do
1298
1299       write(lundiag,'(/a,2i5)') 'kcldbot,  top inp', kcldbot_inp, kcldtop_inp
1300       write(lundiag,'( a,2i5)') 'kcldbot,  top    ', kcldbot, kcldtop
1301       write(lundiag,'( a,2i5)') 'kupdrbot, top    ', kupdrbot, kupdrtop
1302       if ( do_dndraft ) &
1303       write(lundiag,'( a,2i5)') 'kdndrbot, top    ', kdndrbot, kdndrtop
1304       write(lundiag,'(a)')
1305
1306       end if ! (idiagaa > 0)
1307
1308       chem_incu = 0.0
1309       if (iok < 0) then
1310          chem_cupflag = -1
1311          goto 89000
1312       end if
1313
1314
1315 !
1316 ! determine time step for active cloud calculations
1317 ! the mass of air moving out of a layer
1318 !    into updraft and downdraft (through their entrainment)
1319 !    and into adjacent layers (through environment subsidence)
1320 !    cannot exceed the mass of air in the layer
1321 !
1322       rhodz(kts:kte) = rhocen(kts:kte)*dz(kts:kte)
1323       rhodzsum = sum( rhodz(kts:kte) )
1324
1325       mf_ev(kts:ktep1) = -( mf_up(kts:ktep1) + mf_dn(kts:ktep1) )
1326       af_ev(kts:kte) = 1.0_r8 - ( af_up(kts:kte) + af_dn(kts:kte) )
1327
1328       tmpa = 1.0e10   ! tmpa is the smallest of the
1329                       ! individual layer maximum time steps
1330       do k = kupdrbot, kupdrtop
1331          tmpb = mf_up_ent(k) + mf_dn_ent(k) &
1332               + max( mf_ev(k+1), 0.0_r8 ) &
1333               + max( -mf_ev(k),  0.0_r8 )           ! tmpb is air-mass flux out of this layer
1334          tmpc = rhodz(k) / max( tmpb, 1.0e-10_r8 )  ! tmpc is max. time step for this layer
1335          tmpa = min( tmpa, tmpc )                   ! tmpa is the smallest of all the tmpc
1336          if (idiagaa > 0) write(lundiag,'(a,1x,i10,1p,e11.3)') &
1337             'k, dtmax', k, tmpc
1338       end do
1339       tmpd = tmpa
1340       tmpa = tmpa * air_outflow_limit
1341       ntsub = floor( tau_active/tmpa ) + 1
1342       dtsub = tau_active/ntsub
1343       if (idiagaa > 0) then
1344          write(lundiag,'(a,1x,i10,1p,2e11.3)') 'k, dtmax', -1, tmpd, tmpa
1345          write(lundiag,'(a,1x,i10,1p,2e11.3)') &
1346             'ntsub, dtsub, tau_active', ntsub, dtsub, tau_active
1347          write(lundiag,'(2a,1x,10l5)') &
1348             'do_activa, _2ndact, _resusp, _aqchem, ', &
1349             '_wetrem, _updraft, _dndraft', &
1350             do_activa, do_2ndact, do_resusp, do_aqchem, &
1351             do_wetrem, do_updraft, do_dndraft
1352       end if
1353
1354 !
1355 !
1356 ! do multiple substeps of active cloud
1357 !
1358 !
1359       chem_av_new(:,:) = chem(:,:)
1360       zav_chem_av_new(:) = 0.0_r8
1361       do l = p1st, num_chem
1362          zav_chem_av_new(l) = sum( rhodz(kts:kte)*chem_av_new(kts:kte,l) ) / rhodzsum
1363       end do
1364       del_chem_activa(:,:) = 0.0_r8
1365       del_chem_aqchem(:,:) = 0.0_r8
1366       del_chem_residu(:,:) = 0.0_r8
1367       del_chem_resusp(:,:) = 0.0_r8
1368       del_chem_totall(:,:) = 0.0_r8
1369       del_chem_wetrem(:,:) = 0.0_r8
1370       del_chem_ztrans(:,:) = 0.0_r8
1371       del_chem_actvbb(:,:) = 0.0_r8
1372       del_chem_aqchbb(:,:) = 0.0_r8
1373       del_chem_resdbb(:,:) = 0.0_r8
1374       del_chem_respbb(:,:) = 0.0_r8
1375       del_chem_totlbb(:,:) = 0.0_r8
1376
1377       ido_inact = 0
1378
1379 active_cloud_jtsub_loop: &
1380       do jtsub = 1, ntsub
1381
1382
1383 !
1384 ! calculate aerosol and gas profiles in the updraft
1385 !
1386       chem_av_old(:,:) = chem_av_new(:,:)
1387       zav_chem_av_old(:) = zav_chem_av_new(:)
1388
1389       chem_up(:,:) = 0.0_r8
1390       dchemdt_ev_resusp(:,:) = 0.0_r8
1391       dchemdt_up_activa(:,:) = 0.0_r8
1392       dchemdt_up_aqchem(:,:) = 0.0_r8
1393       dchemdt_up_wetrem(:,:) = 0.0_r8
1394
1395       zav_chem_up(:) = 0.0_r8
1396       tmp_mfxchem_up = 0.0_r8
1397
1398 do_updraft_mixratio_calc: &
1399       if ( do_updraft ) then
1400
1401 updraft_mixratio_k_loop: &
1402       do k = kupdrbot, kupdrtop
1403
1404          tmp_mf_up = mf_up(k)
1405
1406 ! do activation at first cloudy layer
1407 ! when first cloudy layer > first updraft layer
1408          if ( do_activa ) then
1409          if ((k == kcldbotliq) .and. (k > kupdrbot)) then
1410             iflagaa = 1
1411             call chem_cup_activate_up( &
1412                lunerr, lundiag, idiagaa, &
1413                kts, kte, p1st, num_chem, &
1414                ktau, grid_id, i, j, k, iflagaa, &
1415                pcen(k), tcen(k), rhocen(k), qcw_inup(k), &
1416                rhodz(k), af_up(k), wact, &
1417                tmp_mf_up, tmp_mfxchem_up, dchemdt_up_activa, &
1418                maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
1419                ncomp_aer, nsize_aer, ntype_aer, &
1420                ai_phase, cw_phase, msectional, &
1421                massptr_aer, numptr_aer, &
1422                dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer )
1423          end if
1424          end if ! ( do_activa ) then
1425
1426 ! do entrainment
1427          if (mf_up_ent(k) > 0.0_r8) then
1428             do l = p1st, num_chem
1429                tmp_mfxchem_up(l) = tmp_mfxchem_up(l) + chem_av_old(k,l)*mf_up_ent(k)
1430             end do
1431             tmp_mf_up = tmp_mf_up + mf_up_ent(k)
1432          end if
1433
1434 ! do activation at first cloudy layer
1435 ! when first cloudy layer = first updraft layer
1436          if ( do_activa ) then
1437          if ((k == kcldbotliq) .and. (k == kupdrbot)) then
1438             iflagaa = 2
1439          else if (( do_2ndact ) .and. (k > kcldbotliq)) then
1440             iflagaa = 10
1441          else
1442             iflagaa = 0
1443          end if
1444          if (iflagaa > 0) then
1445             call chem_cup_activate_up( &
1446                lunerr, lundiag, idiagaa, &
1447                kts, kte, p1st, num_chem, &
1448                ktau, grid_id, i, j, k, iflagaa, &
1449                pcen(k), tcen(k), rhocen(k), qcw_inup(k), &
1450                rhodz(k), af_up(k), wact, &
1451                tmp_mf_up, tmp_mfxchem_up, dchemdt_up_activa, &
1452                maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
1453                ncomp_aer, nsize_aer, ntype_aer, &
1454                ai_phase, cw_phase, msectional, &
1455                massptr_aer, numptr_aer, &
1456                dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer )
1457          end if
1458          end if ! ( do_activa ) then
1459
1460 ! do cloud chemistry
1461          if ( do_aqchem ) then
1462          if (qcw_inup(k) > qcw_cldchem_cutoff) then
1463             tmp_w_up = tmp_mf_up / (af_up(k) * rhocen(k))
1464             tmp_w_up = max( tmp_w_up, 0.1_r8 )
1465             tmp_dt = dz(k)/tmp_w_up
1466
1467 !           subr. chem_cup_aqchem( &
1468 !              config_flags, aer_mech_id, &
1469 !              lunerr, lundiag, idiagaa, &
1470 !              kts, kte, p1st, num_chem, &
1471 !              p_qc, num_moist, &
1472 !              ktau, grid_id, i, j, &
1473 !              iflagaa, ido_aqchem, &
1474 !              dt_aqchem, &
1475 !              pcen, tcen, rhocen, rhodz, qcw, ph_no2, &
1476 !              af_up, tmp_gas_aqfrac_up, tmp_mf_up, tmp_mfxchem_up, &
1477 !              dchemdt_up_aqchem, chem_inact )
1478
1479             call chem_cup_aqchem( &
1480                config_flags, aer_mech_id, &
1481                lunerr, lundiag, idiagaa, &
1482                kts, kte, p1st, num_chem, &
1483                p_qc, num_moist, &
1484                ktau, grid_id, i, j, &
1485                k, ido_inact, &
1486                tmp_dt, &
1487                pcen, tcen, rhocen, rhodz, qcw_inup, ph_no2, &
1488                af_up(k), tmp_gas_aqfrac_up, tmp_mf_up, tmp_mfxchem_up, &
1489                dchemdt_up_aqchem, chem_inact_new )
1490          end if
1491          end if ! ( do_aqchem ) then
1492
1493 ! do wet removal
1494          if ( do_wetrem ) then
1495 ! for now just do
1496 !    fcvt_qc_to_pr transfers cloud-borne to precip, which is assumed to fall out
1497 ! eventually should track ice-borne aerosol and do
1498 !    fcvt_qc_to_qi transfers cloud-borne to ice-borne
1499 !    fcvt_qi_to_pr transfers   ice-borne to precip, which is assumed to fall out
1500 !    detrained ice-borne goes to interstitial (?)
1501             tmpf = min( 1.0_r8, max( 0.0_r8, fcvt_qc_to_pr(k) ) )
1502             if (tmpf > 1.0e-10_r8) then
1503                if ( do_aqchem .and. (qcw_inup(k) > qcw_cldchem_cutoff) ) then
1504                   do l = p1st, num_chem
1505                      if (tmp_gas_aqfrac_up(l) <= 1.0e-10_r8) cycle
1506                      tmpg = min( 1.0_r8, max( 0.0_r8, tmp_gas_aqfrac_up(l)*tmpf ) )
1507                      ! tmpd = wet-removal change to (mf_up * chem_up) as air moves thru the layer
1508                      tmpd = -tmp_mfxchem_up(l)*tmpg
1509                      tmp_mfxchem_up(l) = tmp_mfxchem_up(l) + tmpd
1510                      dchemdt_up_wetrem(k,l) = tmpd/(rhodz(k)*af_up(k))
1511 !                    if (idiagaa > 0) &
1512 !                       write(lundiag,'(a,2i5,1p,2e10.2,2x,a)') &
1513 !                       'gas_aqfrac', k, l, tmp_gas_aqfrac_up(l), qcw_inup(k), chem_name(l)
1514                   end do
1515                end if
1516
1517                do itype = 1, ntype_aer
1518                do isize = 1, nsize_aer(itype)
1519                do icomp = 0, ncomp_aer(itype)
1520                   if (icomp == 0) then
1521                      l = numptr_aer(isize,itype,cw_phase)
1522                   else
1523                      l = massptr_aer(icomp,isize,itype,cw_phase)
1524                   end if
1525                   if ((l < p1st) .or. (l > num_chem)) cycle
1526                   ! tmpd = wet-removal change to (mf_up * chem_up) as air moves thru the layer
1527                   tmpd = -tmp_mfxchem_up(l)*tmpf
1528                   tmp_mfxchem_up(l) = tmp_mfxchem_up(l) + tmpd
1529                   dchemdt_up_wetrem(k,l) = tmpd/(rhodz(k)*af_up(k))
1530                end do
1531                end do
1532                end do
1533             end if ! (tmpf > 1.0e-10_r8)
1534          end if ! ( do_wetrem )
1535
1536
1537 ! calculate updraft mixing ratios after above processes
1538          tmp_chem_up(p1st:num_chem) = tmp_mfxchem_up(p1st:num_chem)/tmp_mf_up
1539
1540 ! do detrainment (which does not change the updraft mixing ratios
1541          tmp_mf_up = max( 0.0_r8, tmp_mf_up - mf_up_det(k) )
1542 !        tmp_mf_up = mf_up(k+1)
1543          do l = p1st, num_chem
1544             tmp_mfxchem_up(l) = tmp_mf_up*tmp_chem_up(l)
1545          end do
1546
1547 ! tmp_chem_up at this point is the chem_up at top of layer k
1548          chem_up(k+1,p1st:num_chem) = tmp_chem_up(p1st:num_chem)
1549
1550 ! calculate resuspension of detrained aerosol
1551 ! this occurs in occurs in the environment, so does not affect chem_up
1552          if ( do_resusp ) then
1553             tmpa = 1.0_r8 - af_lscld(k)   ! fraction resuspended = 1 - ls_cloud_fraction
1554             tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
1555
1556             tmpb = tmpa * mf_up_det(k)/(rhodz(k)*af_ev(k))
1557
1558 !           do l1 = kts, kte
1559 !              l2 = l1 + lgrp*3
1560 !              l3 = l1 + lgrp*4
1561
1562 !              tmpc = tmp_chem_up(l2)*tmpb
1563 !              dchemdt_ev_resusp(k,l2) = dchemdt_ev_resusp(k,l2) - tmpc
1564 !              dchemdt_ev_resusp(k,l3) = dchemdt_ev_resusp(k,l3) + tmpc
1565 !           end do ! l1
1566
1567             do itype = 1, ntype_aer
1568             do isize = 1, nsize_aer(itype)
1569                do l2 = 0, ncomp_aer(itype)
1570                   if (l2 == 0) then
1571                      la = numptr_aer(isize,itype,ai_phase)
1572                      lc = numptr_aer(isize,itype,cw_phase)
1573                   else
1574                      la = massptr_aer(l2,isize,itype,ai_phase)
1575                      lc = massptr_aer(l2,isize,itype,cw_phase)
1576                   end if
1577                   if ((la < p1st) .or. (la > num_chem)) cycle
1578                   if ((lc < p1st) .or. (lc > num_chem)) cycle
1579
1580                   tmpc = tmp_chem_up(lc)*tmpb
1581                   dchemdt_ev_resusp(k,lc) = dchemdt_ev_resusp(k,lc) - tmpc
1582                   dchemdt_ev_resusp(k,la) = dchemdt_ev_resusp(k,la) + tmpc
1583                end do ! l2
1584             end do ! isize
1585             end do ! itype
1586
1587          end if ! ( do_resusp ) then
1588
1589       end do updraft_mixratio_k_loop
1590
1591       do l = 1, num_chem
1592          zav_chem_up(l) = sum( rhodz(kts:kte)*chem_up(kts:kte,l) ) / rhodzsum
1593       end do
1594
1595       end if do_updraft_mixratio_calc
1596
1597
1598       chem_dn(:,:) = 0.0_r8   !BSINGH(03/11/2015): Moved out of the if(do_dndraft) condition to avoid uninitilized chem_dn
1599       zav_chem_dn(:) = 0.0_r8 !BSINGH(03/11/2015): Moved out of the if(do_dndraft) condition to avoid uninitilized zav_chem_dn
1600 !
1601 ! calculate aerosol and gas profiles in the dndraft
1602 !
1603 do_dndraft_mixratio_calc: &
1604       if ( do_dndraft ) then
1605
1606       tmp_mfxchem_dn = 0.0_r8
1607
1608 dndraft_mixratio_k_loop: &
1609       do k = kdndrtop, kdndrbot, -1
1610
1611          tmp_mf_dn = mf_dn(k+1)
1612 ! at this point, tmp_mf_dn and tmp_mvxchem_dn are at the layer top
1613
1614 ! do entrainment
1615          if (mf_dn_ent(k) > 0.0_r8) then
1616             do l = p1st, num_chem
1617                tmp_mfxchem_dn(l) = tmp_mfxchem_dn(l) - chem_av_old(k,l)*mf_dn_ent(k)
1618             end do
1619             tmp_mf_dn = tmp_mf_dn - mf_dn_ent(k)
1620          end if
1621 ! at this point, tmp_mf_dn and tmp_mvxchem_dn are at the layer bottom
1622 ! (but before detrainment)
1623
1624 ! calculate dndraft mixing ratios after above processes
1625          tmp_chem_dn(p1st:num_chem) = tmp_mfxchem_dn(p1st:num_chem)/tmp_mf_dn
1626
1627 ! do detrainment (which does not change the dndraft mixing ratios
1628          tmp_mf_dn = min( 0.0_r8, tmp_mf_dn + mf_dn_det(k) )
1629          do l = p1st, num_chem
1630             tmp_mfxchem_dn(l) = tmp_mf_dn*tmp_chem_dn(l)
1631          end do
1632
1633 ! tmp_chem_dn at this point is the chem_dn at bottom of layer k
1634          chem_dn(k,p1st:num_chem) = tmp_chem_dn(p1st:num_chem)
1635
1636 ! calculate resuspension of detrained aerosol
1637 ! this occurs in occurs in the environment, so does not affect chem_dn
1638          if ( do_resusp ) then
1639 !           tmpa = 1.0_r8 - af_lscld(k)   ! fraction resuspended = 1 - ls_cloud_fraction
1640 !           tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
1641             tmpa = 0.0_r8  ! assume that downdraft is subsaturated at detrainment level
1642                            ! so any cloud-borne aerosol is resuspended
1643
1644             tmpb = tmpa * mf_dn_det(k)/(rhodz(k)*af_ev(k))
1645
1646             do itype = 1, ntype_aer
1647             do isize = 1, nsize_aer(itype)
1648                do l2 = 0, ncomp_aer(itype)
1649                   if (l2 == 0) then
1650                      la = numptr_aer(isize,itype,ai_phase)
1651                      lc = numptr_aer(isize,itype,cw_phase)
1652                   else
1653                      la = massptr_aer(l2,isize,itype,ai_phase)
1654                      lc = massptr_aer(l2,isize,itype,cw_phase)
1655                   end if
1656                   if ((la < p1st) .or. (la > num_chem)) cycle
1657                   if ((lc < p1st) .or. (lc > num_chem)) cycle
1658
1659                   tmpc = tmp_chem_dn(lc)*tmpb
1660                   dchemdt_ev_resusp(k,lc) = dchemdt_ev_resusp(k,lc) - tmpc
1661                   dchemdt_ev_resusp(k,la) = dchemdt_ev_resusp(k,la) + tmpc
1662                end do ! l2
1663             end do ! isize
1664             end do ! itype
1665
1666          end if ! ( do_resusp ) then
1667
1668       end do dndraft_mixratio_k_loop
1669
1670       do l = 1, num_chem
1671          zav_chem_dn(l) = sum( rhodz(kts:kte)*chem_dn(kts:kte,l) ) / rhodzsum
1672       end do
1673
1674       end if do_dndraft_mixratio_calc
1675
1676
1677 !
1678 !
1679 ! now solve mass conservation for timestep dtsub
1680 !
1681 !
1682       if ( do_dndraft ) then
1683          kaa = min( kupdrbot, kdndrbot ) ; kzz = max( kupdrtop, kdndrtop )
1684       else
1685          kaa = kupdrbot ; kzz = kupdrtop
1686       end if
1687
1688       do l = p1st, num_chem
1689
1690       tmp_zflux_top = 0.0_r8
1691
1692       do k = kaa, kzz
1693
1694       tmp_zflux_bot   = tmp_zflux_top
1695
1696       !BSINGH(12/19/2013) - Commented out the following (tmp_zflux_botev, tmp_zflux_botup and tmp_zflux_botdn)
1697       !variables as "tmp_zflux_topev" is undefined at this point and other variables
1698       !(tmp_zflux_botev, tmp_zflux_botup and tmp_zflux_botdn) are acting as dummy variables in the code
1699
1700       !tmp_zflux_botev = tmp_zflux_topev
1701       !tmp_zflux_botup = tmp_zflux_topup
1702       !tmp_zflux_botdn = tmp_zflux_topdn
1703
1704       if (mf_ev(k+1) >= 0.0_r8) then
1705          tmp_zflux_topev = mf_ev(k+1)*chem_av_old(k,l)
1706       else
1707          tmp_zflux_topev = mf_ev(k+1)*chem_av_old(k+1,l)
1708       end if
1709
1710       tmp_zflux_topup = mf_up(k+1)*chem_up(k+1,l)
1711       tmp_zflux_topdn = mf_dn(k+1)*chem_dn(k+1,l)  ! this will be zero when do_dndraft=.false.
1712
1713       tmp_zflux_top = tmp_zflux_topev + tmp_zflux_topup + tmp_zflux_topdn
1714
1715       tmpa = (tmp_zflux_bot - tmp_zflux_top)*dtsub/rhodz(k)
1716       tmpb = dchemdt_up_activa(k,l)*af_up(k)*dtsub
1717       tmpc = dchemdt_up_aqchem(k,l)*af_up(k)*dtsub
1718       tmpd = dchemdt_up_wetrem(k,l)*af_up(k)*dtsub
1719       tmpe = dchemdt_ev_resusp(k,l)*af_ev(k)*dtsub
1720
1721       del_chem_ztrans(k,l) = del_chem_ztrans(k,l) + tmpa
1722       del_chem_activa(k,l) = del_chem_activa(k,l) + tmpb
1723       del_chem_aqchem(k,l) = del_chem_aqchem(k,l) + tmpc
1724       del_chem_wetrem(k,l) = del_chem_wetrem(k,l) + tmpd
1725       del_chem_resusp(k,l) = del_chem_resusp(k,l) + tmpe
1726
1727       chem_av_new(k,l) = chem_av_old(k,l) + (tmpa + tmpb + tmpc + tmpd + tmpe)
1728
1729       if (chem_av_new(k,l) < 0.0_r8) then
1730          del_chem_residu(k,l) = del_chem_residu(k,l) - chem_av_new(k,l)
1731          chem_av_new(k,l) = 0.0_r8
1732       end if
1733       del_chem_totall(k,l) = del_chem_totall(k,l) + (chem_av_new(k,l) - chem_av_old(k,l))
1734
1735 !     if (idiagaa > 0) then
1736 !     if ((l == numptr_aer(2,1,ai_phase)) .or. (l == numptr_aer(2,1,cw_phase))) then
1737 !        tmpveca(1) = (tmp_zflux_botev - tmp_zflux_topev)*dtsub/rhodz(k)
1738 !        tmpveca(2) = (tmp_zflux_botup - tmp_zflux_topup)*dtsub/rhodz(k)
1739 !        tmpveca(3) = (tmp_zflux_botdn - tmp_zflux_topdn)*dtsub/rhodz(k)
1740 !        tmpveca(4) = tmpb
1741 !        tmpveca(5) = tmpc
1742 !        tmpveca(6) = tmpd
1743 !        tmpveca(7) = tmpe
1744 !        write(lundiag,'(i4,1p,7e12.4,3x,2a)') k, tmpveca(1:7), 'tends ', chem_name(l)
1745 !     end if
1746 !     end if
1747
1748       end do ! k
1749
1750       zav_chem_av_new(l)     = sum( rhodz(kts:kte)*chem_av_new(    kts:kte,l) ) / rhodzsum
1751       zav_del_chem_activa(l) = sum( rhodz(kts:kte)*del_chem_activa(kts:kte,l) ) / rhodzsum
1752       zav_del_chem_aqchem(l) = sum( rhodz(kts:kte)*del_chem_aqchem(kts:kte,l) ) / rhodzsum
1753       zav_del_chem_residu(l) = sum( rhodz(kts:kte)*del_chem_residu(kts:kte,l) ) / rhodzsum
1754       zav_del_chem_resusp(l) = sum( rhodz(kts:kte)*del_chem_resusp(kts:kte,l) ) / rhodzsum
1755       zav_del_chem_totall(l) = sum( rhodz(kts:kte)*del_chem_totall(kts:kte,l) ) / rhodzsum
1756       zav_del_chem_wetrem(l) = sum( rhodz(kts:kte)*del_chem_wetrem(kts:kte,l) ) / rhodzsum
1757       zav_del_chem_ztrans(l) = sum( rhodz(kts:kte)*del_chem_ztrans(kts:kte,l) ) / rhodzsum
1758
1759       end do ! l
1760
1761
1762 ! diagnostic outputs
1763       if (idiagaa > 0) then
1764          call chem_cup_1d_diags_pt21( &
1765          lundiag, kdiagbot, kdiagtop, &
1766          kts, kte, ktep1, p1st, num_chem, &
1767          ktau, grid_id, i, j, jtsub, ntsub, &
1768          ishall, do_dndraft, &
1769          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
1770          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
1771          ai_phase, cw_phase, &
1772          massptr_aer, numptr_aer, &
1773          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
1774          rhodz, rhodzsum, &
1775          chem_av_new, chem_av_old, chem_up, chem_dn, &
1776          zav_chem_av_new, zav_chem_av_old, zav_chem_up, zav_chem_dn, &
1777          zav_del_chem_activa, zav_del_chem_aqchem, &
1778          zav_del_chem_residu, zav_del_chem_resusp, &
1779          zav_del_chem_wetrem, zav_del_chem_ztrans, &
1780          zav_del_chem_totall, &
1781          chem_name )
1782       end if ! (idiagaa > 0)
1783
1784
1785       end do active_cloud_jtsub_loop
1786
1787
1788 ! mass balance check after active cloud calcs
1789       itmpa = 0
1790       tmpveca(1:20) = 0.0_r8
1791       do l = p1st, num_chem
1792          tmpa = 0.0_r8 ; tmpb = 0.0_r8 ; tmpd = 0.0_r8
1793          do k = kts, kte
1794             tmpa = tmpa +  rhodz(k)*chem(k,l)
1795             tmpb = tmpb +  rhodz(k)*chem_av_new(k,l)
1796             tmpd = tmpd +  rhodz(k)*(chem_av_new(k,l) - chem(k,l))
1797          end do ! k
1798          tmpa = tmpa/rhodzsum      ! initial column-avg mix-ratio
1799          tmpb = tmpb/rhodzsum      ! ending column-avg mix-ratio
1800          tmpd = tmpd/rhodzsum      ! column-avg mix-ratio change
1801
1802          tmpvecb(11) = zav_del_chem_activa(l)
1803          tmpvecb(12) = zav_del_chem_resusp(l)
1804          tmpvecb(13) = zav_del_chem_aqchem(l)
1805          tmpvecb(14) = zav_del_chem_wetrem(l)
1806          tmpvecb(15) = zav_del_chem_ztrans(l)
1807          tmpvecb(16) = zav_del_chem_residu(l)
1808
1809          tmpe = sum( tmpvecb(11:16) )  ! column-avg m.r. changes by processes
1810          tmpf = sum( max(  tmpvecb(11:16), 0.0 ) )
1811          tmpg = sum( max( -tmpvecb(11:16), 0.0 ) )
1812          tmpf = max( tmpf, tmpg )
1813          tmpg = (tmpd-tmpe)/max( tmpa, tmpb, tmpf, 1.0e-30_r8 )  ! normalized residual
1814
1815          if (abs(tmpg) > abs(tmpveca(1))) then
1816             itmpa = l
1817             tmpveca(1) = tmpg
1818             tmpveca(2) = tmpd-tmpe
1819             tmpveca(3) = tmpd
1820             tmpveca(4) = tmpe
1821             tmpveca(5) = tmpa
1822             tmpveca(6) = tmpb
1823             tmpveca(11:16) = tmpvecb(11:16)
1824          end if
1825          if ((idiagbb > 0) .and. (abs(tmpg) > rerrtol1_mbal)) then
1826             write(lundiagbb,'(/a,i10,3i5,1p,6e11.3,2x,a)') &
1827                'chem_cup_1d massbal       active -', ktau, grid_id, i, j, &
1828                tmpg, (tmpd-tmpe), tmpd, tmpe, tmpa, tmpb, chem_name(l)
1829             write(lundiagbb,'(2a,1p,8e11.3)') &
1830                 'zav_del_chem_activa, resusp, ', &
1831                 'aqchem, wetrem, ztrans, residu', tmpvecb(11:16)
1832          end if
1833       end do ! l
1834
1835       if (idiagaa > 0) then
1836          msg = 'perfect' ; if (itmpa > 0) msg = chem_name(itmpa)
1837          write(lundiag,'(/a,i10,3i5,1p,6e11.3,2x,a)') &
1838             'chem_cup_1d massbal worst active -', &
1839             ktau, grid_id, i, j, tmpveca(1:6), msg(1:12)
1840          write(lundiag,'(2a,1p,8e11.3)') &
1841              'zav_del_chem_activa, resusp, ', &
1842              'aqchem, wetrem, ztrans, residu', tmpveca(11:16)
1843       end if
1844
1845
1846 !
1847 !
1848 ! now do calculations for inactive cloud
1849 !
1850 ! > start with chem_av_new in the entire grid-area
1851 ! > divide the grid area into
1852 !    > inactive cloud portion with area fraction = af_cucld(k) - af_up(k)
1853 !    > "other" portion = the rest
1854 ! > in inactive cloud portion
1855 !    > start with chem_inact_old = chem_av_new
1856 !      (cannot use chem_up here because in the steady-state updraft approach,
1857 !       the updrafts process aerosol/gas mass,
1858 !       but they do not really contain any aerosol/gas mass)
1859 !    > set interstitial/activated fractions based on
1860 !         the fractions in the updraft
1861 !    > do cloud chemistry
1862 !    > do partial aerosol resuspension before "dissipating"
1863 !         the inactive cloud back to the grid-average
1864 ! > "dissipate" the inactive cloud back to the grid-average
1865 !
1866       af_inact = 0.0
1867       ido_inact = 0
1868       chem_inact_new = -1.0e10
1869       chem_inact_dsp = -1.0e10
1870       if ( .not. do_inact ) goto 79000
1871
1872       itmpa = 0
1873       do k = kcldbot, kcldtop
1874          ! count levels that have qcw>0 both in cumulus and in updraft
1875          ! and also have cumulus area fraction > updraft area fraction
1876          tmpa = af_cucld(k)-af_up(k)
1877          if ( (qcw_incu(k) > 0.0) .and. &
1878               (qcw_inup(k) > 0.0) .and. &
1879               (tmpa >= af_cucld_smallaa*0.5) ) then
1880             af_inact(k) = max( 0.0_r8, min( 1.0_r8, tmpa ) )
1881             ido_inact(k) = 1
1882             itmpa = itmpa + 1
1883          end if
1884       end do
1885
1886       if (itmpa <= 0)then
1887          if (idiagaa > 0) write(lundiag,'(/a,4i10)') &
1888             'chem_cup_1d - no inactive cloud calcs - ktau, id, i, j =', &
1889             ktau, grid_id, i, j
1890
1891          chem(:,:) = chem_av_new(:,:)  ! put updated mix-ratios into chem
1892          goto 79000
1893       end if
1894
1895
1896 ! set initial mixing ratios for inactive cloud
1897       chem_av_old(:,:) = chem_av_new(:,:)
1898       chem_inact_old(:,:) = chem_av_new(:,:)
1899
1900       if (idiagaa > 0) write(lundiag,'(//a)') 'inactive k, fmact, fnact'
1901       do k = kcldtop, kcldbot, -1
1902          if (ido_inact(k) <= 0) cycle
1903
1904          tmp_fmact = 0.0_r8 ; tmp_fnact = 0.0_r8
1905          do itype = 1, ntype_aer
1906          do isize = 1, nsize_aer(itype)
1907             tmpg = 0.0_r8
1908             do l2 = 0, ncomp_aer(itype)
1909                if (l2 == 0) then
1910                   la = numptr_aer(isize,itype,ai_phase)
1911                   lc = numptr_aer(isize,itype,cw_phase)
1912                else
1913                   la = massptr_aer(l2,isize,itype,ai_phase)
1914                   lc = massptr_aer(l2,isize,itype,cw_phase)
1915                end if
1916                if ((la < p1st) .or. (la > num_chem)) cycle
1917                if ((lc < p1st) .or. (lc > num_chem)) cycle
1918
1919                tmpa = max( chem_up(k+1,la), 1.0e-35_r8 )
1920                tmpc = max( chem_up(k+1,lc), 1.0e-35_r8 )
1921                tmpd = tmpa + tmpc
1922                tmpe = max( chem_inact_old(k,la) + chem_inact_old(k,lc), 0.0_r8 )
1923                ! apply interstitial/activated fractions from the updraft
1924                chem_inact_old(k,la) = tmpe*(tmpa/tmpd)
1925                chem_inact_old(k,lc) = tmpe*(tmpc/tmpd)
1926
1927                del_chem_actvbb(k,la) = af_inact(k)*(chem_inact_old(k,la) - chem_av_new(k,la))
1928                del_chem_actvbb(k,lc) = af_inact(k)*(chem_inact_old(k,lc) - chem_av_new(k,lc))
1929
1930                if (l2 == 0) then
1931                   tmp_fnact(isize,itype) = tmpc/tmpd
1932                else
1933                   tmpe = max( tmpe, 1.0e-35_r8 )
1934                   tmp_fmact(isize,itype) = tmp_fmact(isize,itype) + tmpe*(tmpc/tmpd)
1935                   tmpg = tmpg + tmpe
1936                end if
1937             end do ! l2
1938             tmp_fmact(isize,itype) = tmp_fmact(isize,itype)/tmpg
1939          end do ! isize
1940          if (idiagaa > 0) &
1941             write(lundiag,'( i3,2(2x,8f8.4) / (3x,2(2x,8f8.4)) )') k, &
1942                tmp_fmact(1:nsize_aer(itype),itype), tmp_fnact(1:nsize_aer(itype),itype)
1943 ! Added ',itype' in tmp_fmact and tmp_fnact arrays by Manish Shrivastava
1944          end do ! itype
1945
1946          chem_inact_new(k,:) = chem_inact_old(k,:)
1947       end do ! k
1948
1949
1950       if ( do_aqchem ) then
1951
1952 ! do cloud chemistry
1953 !     subr. chem_cup_aqchem( &
1954 !        config_flags, aer_mech_id, &
1955 !        lunerr, lundiag, idiagaa, &
1956 !        kts, kte, p1st, num_chem, &
1957 !        p_qc, num_moist, &
1958 !        ktau, grid_id, i, j, &
1959 !        iflagaa, ido_aqchem, &
1960 !        dt_aqchem, &
1961 !        pcen, tcen, rhocen, rhodz, qcw, ph_no2, &
1962 !        af_up, tmp_gas_aqfrac_up, tmp_mf_up, tmp_mfxchem_up, &
1963 !        dchemdt_up_aqchem, chem_inact )
1964
1965       itmpb = kts-1
1966       call chem_cup_aqchem( &
1967          config_flags, aer_mech_id, &
1968          lunerr, lundiag, idiagaa, &
1969          kts, kte, p1st, num_chem, &
1970          p_qc, num_moist, &
1971          ktau, grid_id, i, j, &
1972          itmpb, ido_inact, &
1973          tau_inactive, &
1974          pcen, tcen, rhocen, rhodz, qcw_incu, ph_no2, &
1975          af_up(kts), tmp_gas_aqfrac_up, tmp_mf_up, tmp_mfxchem_up, &
1976          dchemdt_up_aqchem, chem_inact_new )
1977
1978       do k = kcldbot, kcldtop
1979          if (ido_inact(k) <= 0) cycle
1980          do l = p1st, num_chem
1981             del_chem_aqchbb(k,l) = af_inact(k)*(chem_inact_new(k,l) - chem_inact_old(k,l))
1982          end do
1983       end do
1984
1985       end if ! ( do_aqchem )
1986
1987
1988 ! do "dissipation" of the cumulus cloud,
1989 !    because there is no "memory" for the aerosols/gases within the cumulus cloud
1990 !
1991 ! as part of this, do partial resuspension of activated aerosols in the inactive cloud
1992 !    part of the cumulus cloud goes to grid-resolved cloud, and part goes to clear air
1993 !    if the grid-resolved cloud fraction is 0.0,
1994 !       then all  of the activated aerosols in the cumulus cloud are resuspended
1995 !    if the grid-resolved cloud fraction is 1.0,
1996 !       then none of the activated aerosols in the cumulus cloud are resuspended
1997
1998       chem_inact_dsp = chem_inact_new
1999       do k = kcldbot, kcldtop
2000          if (ido_inact(k) <= 0) cycle
2001
2002 ! first do the resuspension
2003          if ( do_resusp ) then
2004             tmpa = 1.0_r8 - af_lscld(k)   ! fraction resuspended = 1 - ls_cloud_fraction
2005             tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
2006
2007             if (tmpa > 0.0) then
2008             do itype = 1, ntype_aer
2009             do isize = 1, nsize_aer(itype)
2010                do l2 = 0, ncomp_aer(itype)
2011                   if (l2 == 0) then
2012                      la = numptr_aer(isize,itype,ai_phase)
2013                      lc = numptr_aer(isize,itype,cw_phase)
2014                   else
2015                      la = massptr_aer(l2,isize,itype,ai_phase)
2016                      lc = massptr_aer(l2,isize,itype,cw_phase)
2017                   end if
2018                   if ((la < p1st) .or. (la > num_chem)) cycle
2019                   if ((lc < p1st) .or. (lc > num_chem)) cycle
2020
2021                   tmpc = chem_inact_new(k,lc)*tmpa
2022                   chem_inact_dsp(k,lc) = chem_inact_new(k,lc) - tmpc
2023                   chem_inact_dsp(k,la) = chem_inact_new(k,la) + tmpc
2024                   del_chem_respbb(k,lc) = -tmpc*af_inact(k)
2025                   del_chem_respbb(k,la) = tmpc*af_inact(k)
2026                end do ! l2
2027             end do ! isize
2028             end do ! itype
2029             end if ! (tmpa > 0.0)
2030          end if ! ( do_resusp )
2031
2032 ! now the dissipation
2033          tmpb = 1.0 - af_inact(k)
2034          do l = p1st, num_chem
2035             chem_av_new(k,l) = tmpb*chem_av_old(k,l) + af_inact(k)*chem_inact_dsp(k,l)
2036             del_chem_totlbb(k,l) = chem_av_new(k,l) - chem_av_old(k,l)
2037             del_chem_resdbb(k,l) = del_chem_totlbb(k,l) &
2038                - ( del_chem_actvbb(k,l) + del_chem_aqchbb(k,l) + del_chem_respbb(k,l) )
2039          end do ! l
2040
2041       end do ! k
2042
2043       do l = p1st, num_chem
2044       zav_del_chem_actvbb(l) = sum( rhodz(kts:kte)*del_chem_actvbb(kts:kte,l) ) / rhodzsum
2045       zav_del_chem_aqchbb(l) = sum( rhodz(kts:kte)*del_chem_aqchbb(kts:kte,l) ) / rhodzsum
2046       zav_del_chem_resdbb(l) = sum( rhodz(kts:kte)*del_chem_resdbb(kts:kte,l) ) / rhodzsum
2047       zav_del_chem_respbb(l) = sum( rhodz(kts:kte)*del_chem_respbb(kts:kte,l) ) / rhodzsum
2048       zav_del_chem_totlbb(l) = sum( rhodz(kts:kte)*del_chem_totlbb(kts:kte,l) ) / rhodzsum
2049       end do
2050
2051
2052 ! diagnostics
2053       if (idiagaa > 0) then
2054          call chem_cup_1d_diags_pt41( &
2055          lundiag, 1, &
2056          kts, kte, ktep1, p1st, num_chem, &
2057          ktau, grid_id, i, j, &
2058          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2059          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
2060          ai_phase, cw_phase, &
2061          massptr_aer, numptr_aer, &
2062          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
2063          rhodz, rhodzsum, &
2064          chem_av_old, chem_av_new, &
2065          zav_del_chem_actvbb, zav_del_chem_aqchbb, &
2066          zav_del_chem_resdbb, zav_del_chem_respbb, &
2067          zav_del_chem_totlbb, zav_del_chem_wetrem, &
2068          chem_name )
2069       end if ! (idiagaa > 0)
2070
2071
2072 ! mass balance check after inactive cloud calcs
2073       itmpa = 0
2074       tmpveca(1:20) = 0.0
2075       do l = p1st, num_chem
2076          tmpa = 0.0_r8 ; tmpb = 0.0_r8 ; tmpd = 0.0_r8
2077          do k = kts, kte
2078             tmpa = tmpa +  rhodz(k)*chem(k,l)
2079             tmpb = tmpb +  rhodz(k)*chem_av_new(k,l)
2080             tmpd = tmpd +  rhodz(k)*(chem_av_new(k,l) - chem(k,l))
2081          end do ! k
2082          tmpa = tmpa/rhodzsum      ! initial column-avg mix-ratio
2083          tmpb = tmpb/rhodzsum      ! ending column-avg mix-ratio
2084          tmpd = tmpd/rhodzsum      ! column-avg mix-ratio change
2085
2086          tmpvecb(11) = zav_del_chem_activa(l)
2087          tmpvecb(12) = zav_del_chem_resusp(l)
2088          tmpvecb(13) = zav_del_chem_aqchem(l)
2089          tmpvecb(14) = zav_del_chem_wetrem(l)
2090          tmpvecb(15) = zav_del_chem_ztrans(l)
2091          tmpvecb(16) = zav_del_chem_residu(l)
2092          tmpvecb(17) = zav_del_chem_actvbb(l)
2093          tmpvecb(18) = zav_del_chem_respbb(l)
2094          tmpvecb(19) = zav_del_chem_aqchbb(l)
2095
2096          tmpe = sum( tmpvecb(11:14) )  + sum( tmpvecb(17:19) )  ! column-avg m.r. changes by processes
2097          tmpf = sum( max(  tmpvecb(11:14), 0.0 ) ) + sum( max(  tmpvecb(17:19), 0.0 ) )
2098          tmpg = sum( max( -tmpvecb(11:14), 0.0 ) ) + sum( max( -tmpvecb(17:19), 0.0 ) )
2099          tmpf = max( tmpf, tmpg )
2100          tmpg = (tmpd-tmpe)/max( tmpa, tmpb, tmpf, 1.0e-30_r8 )  ! normalized residual
2101
2102          if (abs(tmpg) > abs(tmpveca(1))) then
2103             itmpa = l
2104             tmpveca(1) = tmpg
2105             tmpveca(2) = tmpd-tmpe
2106             tmpveca(3) = tmpd
2107             tmpveca(4) = tmpe
2108             tmpveca(5) = tmpa
2109             tmpveca(6) = tmpb
2110             tmpveca(11:19) = tmpvecb(11:19)
2111          end if
2112          if ((idiagbb > 0) .and. (abs(tmpg) > rerrtol1_mbal)) then
2113             write(lundiagbb,'(/a,i10,3i5,1p,6e11.3,2x,a)') &
2114                'chem_cup_1d massbal       final  -', ktau, grid_id, i, j, &
2115                tmpg, (tmpd-tmpe), tmpd, tmpe, tmpa, tmpb, chem_name(l)
2116             write(lundiagbb,'(2a,1p,8e11.3)') &
2117                 'zav_del_chem_activa, resusp, ', &
2118                 'aqchem, wetrem, ztrans, residu', tmpvecb(11:16)
2119             write(lundiagbb,'(2a,1p,8e11.3)') &
2120                 'zav_del_chem_actvbb, respbb, ', &
2121                 'aqchbb                        ', tmpvecb(17:19)
2122          end if
2123       end do ! l
2124
2125       if (idiagaa > 0) then
2126          msg = 'perfect' ; if (itmpa > 0) msg = chem_name(itmpa)
2127          write(lundiag,'(/a,i10,3i5,1p,6e11.3,2x,a)') &
2128             'chem_cup_1d massbal worst final  -', &
2129             ktau, grid_id, i, j, tmpveca(1:6), msg(1:12)
2130          write(lundiag,'(2a,1p,8e11.3)') &
2131              'zav_del_chem_activa, resusp, ', &
2132              'aqchem, wetrem, ztrans, residu', tmpveca(11:16)
2133          write(lundiag,'(2a,1p,8e11.3)') &
2134              'zav_del_chem_actvbb, respbb, ', &
2135              'aqchbb                        ', tmpveca(17:19)
2136       end if
2137
2138
2139 ! diagnostics
2140       if (idiagaa > 0) then
2141          call chem_cup_1d_diags_pt41( &
2142          lundiag, 2, &
2143          kts, kte, ktep1, p1st, num_chem, &
2144          ktau, grid_id, i, j, &
2145          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2146          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
2147          ai_phase, cw_phase, &
2148          massptr_aer, numptr_aer, &
2149          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
2150          rhodz, rhodzsum, &
2151          chem, chem_av_new, &
2152          zav_del_chem_actvbb, zav_del_chem_aqchbb, &
2153          zav_del_chem_resdbb, zav_del_chem_respbb, &
2154          zav_del_chem_totlbb, zav_del_chem_wetrem, &
2155          chem_name )
2156       end if ! (idiagaa > 0)
2157
2158
2159 ! put final mixing ratios into chem
2160       chem(:,:) = chem_av_new(:,:)
2161
2162
2163 79000 continue ! "goto 79000" means no inactive cloud calcs
2164
2165
2166 ! calculate average in-convective-cloud mixing ratios
2167 ! The code below calculates average in-convective-cloud mixing ratios
2168 ! Commented out by Manish Shrivastava on 05/14/2013 because we wanted to look at just the active part values of CO_IC_CUP and other  variables
2169 !--------------------------------------------------------------------------------
2170 !      do k = kcldbot, kcldtop
2171 !         if (af_up(k) < 0.5*af_up_smallaa) then
2172
2173 !            if (ido_inact(k) <= 0) then
2174 !               cycle  ! note that chem_incu and chem_cupflag are initialized
2175 !                      ! to zero, so those will be the values at this level
2176 !            else
2177 !               tmpa = 0.0_r8 ; tmpb = 1.0_r8
2178 !            end if
2179 !         else
2180 !            if (ido_inact(k) <= 0) then
2181 !               tmpa = 1.0_r8 ; tmpb = 0.0_r8
2182 !            else
2183 !               tmpa = af_up(k)/(af_up(k) + af_inact(k))
2184 !               tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
2185 !               tmpb = 1.0_r8 - tmpa
2186 !            end if
2187 !         end if
2188 !---------------------------------------------------------------------------------
2189
2190 !The code below was changed by Easter Richard to calculate just active part cloud values
2191 !It has changes so that the chem_incu will be equal to the updraft mixing ratios,
2192 !Rather than area-weighted average of updraft and inactive cloud.
2193
2194 !------------------------------------------------------------------------------------
2195       do k = kcldbot, kcldtop
2196          if (af_up(k) < 0.5*af_up_smallaa) then
2197
2198             cycle   ! change to make chem_incu = chem_up (but no updraft at this k)
2199          else
2200             tmpa = 1.0_r8 ; tmpb = 0.0_r8   ! change to make chem_incu = chem_up
2201          end if
2202 !-----------------------------------------------------------------------------------------------
2203          ! chem_up(k,:) is at top of layer k
2204          ! when k>kcldbot, use an average of the k and k+1 values
2205          ! when k=kcldbot, use the k+1 values because the k values are essentially
2206          !                 clear air (at the top of the last sub-cloud-base layer)
2207          if (k == kcldbot) then
2208             chem_incu(k,:) = tmpa*chem_up(k+1,:) + tmpb*chem_inact_new(k,:)
2209          else
2210             chem_incu(k,:) = tmpa*0.5*(chem_up(k,:)+chem_up(k+1,:)) + tmpb*chem_inact_new(k,:)
2211          end if
2212          chem_cupflag(k) = 1
2213       end do
2214
2215       if (idiagaa > 0) then
2216          call chem_cup_1d_diags_pt71( &
2217          lundiag, kcldbot, kcldtop, &
2218          kts, kte, ktep1, p1st, num_chem, &
2219          ktau, grid_id, i, j, &
2220          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2221          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
2222          ai_phase, cw_phase, &
2223          massptr_aer, numptr_aer, &
2224          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
2225          af_up, af_inact, &
2226          chem_av_new, chem_up, chem_inact_new, chem_incu, &
2227          chem_name )
2228       end if
2229
2230
2231 89000 continue ! "goto 89000" means check_adjust_inputs failed
2232       if (idiagaa > 0) write(lundiag,'(/a,4i10)') &
2233          'chem_cup_1d done  ktau, id, i, j =', ktau, grid_id, i, j
2234
2235       return
2236       end subroutine chem_cup_1d
2237
2238
2239 !-------------------------------------------------------------------------------
2240       subroutine chem_cup_1d_diags_pt21( &
2241          lundiag_inp, kdiagbot, kdiagtop, &
2242          kts, kte, ktep1, p1st, num_chem, &
2243          ktau, grid_id, i, j, jtsub, ntsub, &
2244          ishall, do_dndraft, &
2245          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2246          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
2247          ai_phase, cw_phase, &
2248          massptr_aer, numptr_aer, &
2249          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
2250          rhodz, rhodzsum, &
2251          chem_av_new, chem_av_old, chem_up, chem_dn, &
2252          zav_chem_av_new, zav_chem_av_old, zav_chem_up, zav_chem_dn, &
2253          zav_del_chem_activa, zav_del_chem_aqchem, &
2254          zav_del_chem_residu, zav_del_chem_resusp, &
2255          zav_del_chem_wetrem, zav_del_chem_ztrans, &
2256          zav_del_chem_totall, &
2257          chem_name )
2258
2259       use module_configure, only:  &
2260          p_qc, &
2261          p_h2o2, p_hcl, p_hno3, p_nh3, p_so2, p_sulf
2262
2263 ! arguments
2264       integer, intent(in) :: lundiag_inp
2265       integer, intent(in) :: kts, kte, ktep1, p1st, num_chem
2266       integer, intent(in) :: ktau, grid_id, i, j, jtsub, ntsub
2267       integer, intent(in) :: kdiagbot, kdiagtop
2268       integer, intent(in) :: ishall
2269
2270       integer, intent(in) :: maxd_acomp, maxd_aphase, maxd_atype, maxd_asize
2271       integer, intent(in) :: &
2272          nphase_aer, ntype_aer, &
2273          ai_phase, cw_phase, &
2274          nsize_aer( maxd_atype ), &
2275          ncomp_aer( maxd_atype ), &
2276          massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
2277          numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
2278          lptr_so4_aer(maxd_asize, maxd_atype, maxd_aphase), &
2279          lptr_oin_aer(maxd_asize, maxd_atype, maxd_aphase), &
2280          lptr_no3_aer(maxd_asize, maxd_atype, maxd_aphase), &
2281          lptr_nh4_aer(maxd_asize, maxd_atype, maxd_aphase), &
2282          lptr_cl_aer(maxd_asize, maxd_atype, maxd_aphase)
2283
2284       logical, intent(in) :: do_dndraft
2285
2286       real(r8), intent(in) :: rhodzsum
2287
2288       real(r8), intent(in), dimension(kts:kte) :: rhodz
2289
2290       real(r8), intent(in), dimension(kts:kte,1:num_chem) :: &
2291          chem_av_new, chem_av_old
2292
2293       real(r8), intent(in), dimension(kts:ktep1,1:num_chem) :: chem_up, chem_dn
2294
2295       real(r8), intent(in), dimension(1:num_chem) :: &
2296          zav_chem_av_new, zav_chem_av_old, zav_chem_up, zav_chem_dn, &
2297          zav_del_chem_activa, zav_del_chem_aqchem, &
2298          zav_del_chem_residu, zav_del_chem_resusp, &
2299          zav_del_chem_wetrem, zav_del_chem_ztrans, &
2300          zav_del_chem_totall
2301
2302       character(len=12), intent(in) :: chem_name(num_chem)
2303
2304 ! local variables
2305       integer, parameter :: mxg_max=100, nxg_max=8
2306       integer :: aip, cwp, typ
2307       integer :: isize, itype, itmpa
2308       integer :: l, l2, l3, lundiag
2309       integer :: lxg(nxg_max,mxg_max)
2310       integer :: k
2311       integer :: m, mxg
2312       integer :: n, nxg(mxg_max)
2313
2314       real(r8) :: fxg(nxg_max,mxg_max)
2315       real(r8) :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg
2316       real(r8) :: tmpveca(num_chem)
2317
2318       character(len=80)  :: fmtaa, fmtbb, fmtcc, fmtdd, fmtee
2319
2320
2321       lundiag = lundiag_inp
2322 !     lundiag = 122
2323
2324       m = 1
2325       nxg(m) = 3
2326       lxg(1:nxg(m),m) = (/ p_h2o2, p_so2, p_sulf /)
2327       fxg(1:nxg(m),m) = 1.0e3
2328       m = m + 1
2329       nxg(m) = 2
2330       lxg(1:nxg(m),m) = (/ p_nh3, p_hno3 /)
2331       fxg(1:nxg(m),m) = 1.0e3
2332
2333       aip = ai_phase ; cwp = cw_phase
2334       typ = 1 ! aerosol type
2335       if (cw_phase > 0) then
2336          n = 2 ! size bin
2337          m = m + 1
2338          nxg(m) = 2
2339          lxg(1:nxg(m),m) = (/ lptr_so4_aer(n,typ,aip), lptr_so4_aer(n,typ,cwp) /)
2340          fxg(1:nxg(m),m) = 28.966/96.0
2341          m = m + 1
2342          nxg(m) = 2
2343          lxg(1:nxg(m),m) = (/ lptr_nh4_aer(n,typ,aip), lptr_nh4_aer(n,typ,cwp) /)
2344          fxg(1:nxg(m),m) = 28.966/18.0
2345          m = m + 1
2346          nxg(m) = 2
2347          lxg(1:nxg(m),m) = (/ lptr_no3_aer(n,typ,aip), lptr_no3_aer(n,typ,cwp) /)
2348          fxg(1:nxg(m),m) = 28.966/62.0
2349          m = m + 1
2350          nxg(m) = 2
2351          lxg(1:nxg(m),m) = (/ numptr_aer(n,typ,aip), numptr_aer(n,typ,cwp) /)
2352          fxg(1:nxg(m),m) = 1.0e-6
2353       else
2354          n = 2 ! size bin
2355          m = m + 1
2356          nxg(m) = 2
2357          lxg(1:nxg(m),m) = (/ lptr_so4_aer(n,typ,aip), lptr_nh4_aer(n,typ,aip) /)
2358          fxg(1:nxg(m),m) = (/ 28.966/96.0, 28.966/18.0 /)
2359          m = m + 1
2360          nxg(m) = 2
2361          lxg(1:nxg(m),m) = (/ lptr_no3_aer(n,typ,aip), numptr_aer(n,typ,aip) /)
2362          fxg(1:nxg(m),m) = (/ 28.966/62.0, 1.0e-6 /)
2363       end if
2364       mxg = m
2365
2366
2367       write(lundiag,'(/a,i10,5i5)') &
2368          'chem_cup_1d_diags_pt21 - ktau, id, i, j =', ktau, grid_id, i, j
2369       do m = 1, mxg
2370          n = nxg(m)
2371          write(lundiag,'(a9,2i3,3x,5(i5, 8x  ))') 'm, n, lxg', m, n, (lxg(l2,m), l2=1,n)
2372          write(lundiag,'(15x,   3x,5(a12,1x  ))') (chem_name(lxg(l2,m)), l2=1,n)
2373          write(lundiag,'(15x,   5(1p,e13.3))') (fxg(l2,m), l2=1,n)
2374       end do
2375       write(lundiag,'(/2a)') '*** units for following:  ', &
2376          'trace gas and aerosol mass = ppb,  aerosol number = #/mg'
2377
2378
2379       write(lundiag,'(/a,i10,5i5)') &
2380          'ktau, jtsub, ntsub, id, i, j =', ktau, jtsub, ntsub, grid_id, i, j
2381
2382       do m = 1, mxg
2383 !        cycle
2384          n = nxg(m)
2385          if ( do_dndraft ) then
2386             if (n == 3) then
2387                fmtaa = '(/4x,    4(3x,a33  ))'
2388                fmtbb = '( 4x,    4(3x,3a11 ))'
2389                fmtcc = '( i2, 1p,4(3x,3e11.3))'
2390                fmtdd = '(90x,a20,1p,3x,3e11.3)'
2391             else
2392                fmtaa = '(/4x,    4(3x,a22  ))'
2393                fmtbb = '( 4x,    4(3x,2a11 ))'
2394                fmtcc = '( i2, 1p,4(3x,2e11.3))'
2395                fmtdd = '(57x,a20,1p,3x,2e11.3)'
2396             end if
2397          else
2398             if (n == 3) then
2399                fmtaa = '(/4x,    3(3x,a33  ))'
2400                fmtbb = '( 4x,    3(3x,3a11 ))'
2401                fmtcc = '( i2, 1p,3(3x,3e11.3))'
2402                fmtdd = '(54x,a20,1p,3x,3e11.3)'
2403             else
2404                fmtaa = '(/4x,    3(3x,a22  ))'
2405                fmtbb = '( 4x,    3(3x,2a11 ))'
2406                fmtcc = '( i2, 1p,3(3x,2e11.3))'
2407                fmtdd = '(32x,a20,1p,3x,2e11.3)'
2408             end if
2409          end if
2410
2411          if ( do_dndraft ) then
2412             write(lundiag,fmtaa) &
2413                'chem_up                                 ', &
2414                'chem_dn                                 ', &
2415                'chem_av_old                             ', &
2416                'chem_av_new                             '
2417             write(lundiag,fmtbb) &
2418                ( ( chem_name(lxg(l2,m)), l2=1,n ), l3=1,4 )
2419          else
2420             write(lundiag,fmtaa) &
2421                'chem_up                                 ', &
2422                'chem_av_old                             ', &
2423                'chem_av_new                             '
2424             write(lundiag,fmtbb) &
2425                ( ( chem_name(lxg(l2,m)), l2=1,n ), l3=1,3 )
2426          end if
2427
2428          itmpa = 0
2429          do l2 = 1, n
2430             l = lxg(l2,n)
2431             if (zav_chem_up(l) /= 0.0_r8 ) then
2432                itmpa = 1
2433                cycle
2434             end if
2435             tmpa = 0.0
2436             do k = kts, kte
2437                tmpa = max( tmpa, abs( chem_av_new(k,l)-chem_av_old(k,l) ) )
2438             end do
2439             if (tmpa /= 0.0_r8 ) itmpa = 1
2440          end do
2441
2442          do k = kdiagtop, kdiagbot, -1
2443             if (itmpa == 0) cycle
2444             if ( do_dndraft ) then
2445                write(lundiag,fmtcc) k, &
2446                  ( chem_up(    k,lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2447                  ( chem_dn(    k,lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2448                  ( chem_av_old(k,lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2449                  ( chem_av_new(k,lxg(l2,m))*fxg(l2,m) , l2=1,n )
2450             else
2451                write(lundiag,fmtcc) k, &
2452                  ( chem_up(    k,lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2453                  ( chem_av_old(k,lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2454                  ( chem_av_new(k,lxg(l2,m))*fxg(l2,m) , l2=1,n )
2455             end if
2456          end do ! k
2457          if (itmpa > 0) write(lundiag,'(a)')
2458          if ( do_dndraft ) then
2459             write(lundiag,fmtcc) -1, &
2460                ( zav_chem_up(    lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2461                ( zav_chem_dn(    lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2462                ( zav_chem_av_old(lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2463                ( zav_chem_av_new(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2464          else
2465             write(lundiag,fmtcc) -1, &
2466                ( zav_chem_up(    lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2467                ( zav_chem_av_old(lxg(l2,m))*fxg(l2,m) , l2=1,n ), &
2468                ( zav_chem_av_new(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2469          end if
2470          write(lundiag,fmtdd) &
2471              'zav_del_chem_totall ', &
2472             ( zav_del_chem_totall(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2473          if (itmpa == 0) cycle
2474          write(lundiag,fmtdd) &
2475              'zav_del_chem_aqchem ', &
2476             ( zav_del_chem_aqchem(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2477          write(lundiag,fmtdd) &
2478              'zav_del_chem_wetrem ', &
2479             ( zav_del_chem_wetrem(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2480          write(lundiag,fmtdd) &
2481              'zav_del_chem_activa ', &
2482             ( zav_del_chem_activa(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2483          write(lundiag,fmtdd) &
2484              'zav_del_chem_resusp ', &
2485             ( zav_del_chem_resusp(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2486          write(lundiag,fmtdd) &
2487              'zav_del_chem_ztrans ', &
2488             ( zav_del_chem_ztrans(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2489          write(lundiag,fmtdd) &
2490              'zav_del_chem_residu ', &
2491             ( zav_del_chem_residu(lxg(l2,m))*fxg(l2,m) , l2=1,n )
2492       end do ! m = 1, mxg
2493
2494
2495       tmpveca(p1st:num_chem) = zav_del_chem_aqchem(p1st:num_chem)
2496       write(lundiag,'(/a)') 'zav_del_chem_aqchem summary'
2497
2498       write(lundiag,'(a,1p,10e12.4)') 'h2o2                 ', &
2499          tmpveca(p_h2o2)*1.0e3
2500
2501       tmpb = tmpveca(p_so2) + tmpveca(p_sulf)
2502       write(lundiag,'(a,1p,10e12.4)') 'so2+h2so4, individ   ', &
2503          tmpb*1.0e3, &
2504          tmpveca(p_so2)*1.0e3, &
2505          tmpveca(p_sulf)*1.0e3
2506       itype = 1
2507       tmpa = 28.966/96.0 ; tmpb = 0.0
2508       do isize = 1, nsize_aer(itype)
2509          tmpb = tmpb + tmpveca(lptr_so4_aer(isize,itype,cw_phase))
2510       end do
2511       write(lundiag,'(a,1p,10e12.4)') 'so4_cw total, individ', &
2512          tmpb*tmpa, &
2513          ( tmpveca(lptr_so4_aer(isize,itype,cw_phase))*tmpa, isize=1,nsize_aer(itype) )
2514
2515       write(lundiag,'(a,1p,10e12.4)') 'nh3                  ', &
2516          tmpveca(p_nh3)*1.0e3
2517       tmpa = 28.966/18.0 ; tmpb = 0.0
2518       do isize = 1, nsize_aer(itype)
2519          tmpb = tmpb + tmpveca(lptr_nh4_aer(isize,itype,cw_phase))
2520       end do
2521       write(lundiag,'(a,1p,10e12.4)') 'nh4_cw total, individ', &
2522          tmpb*tmpa, &
2523          ( tmpveca(lptr_nh4_aer(isize,itype,cw_phase))*tmpa, isize=1,nsize_aer(itype) )
2524
2525       write(lundiag,'(a,1p,10e12.4)') 'hno3                 ', &
2526          tmpveca(p_hno3)*1.0e3
2527       tmpa = 28.966/62.0 ; tmpb = 0.0
2528       do isize = 1, nsize_aer(itype)
2529          tmpb = tmpb + tmpveca(lptr_no3_aer(isize,itype,cw_phase))
2530       end do
2531       write(lundiag,'(a,1p,10e12.4)') 'no3_cw total, individ', &
2532          tmpb*tmpa, &
2533          ( tmpveca(lptr_no3_aer(isize,itype,cw_phase))*tmpa, isize=1,nsize_aer(itype) )
2534
2535       write(lundiag,'(a,1p,10e12.4)') 'hcl                  ', &
2536          tmpveca(p_hcl)*1.0e3
2537       tmpa = 28.966/35.5 ; tmpb = 0.0
2538       do isize = 1, nsize_aer(itype)
2539          tmpb = tmpb + tmpveca(lptr_cl_aer(isize,itype,cw_phase))
2540       end do
2541       write(lundiag,'(a,1p,10e12.4)') 'cl_cw  total, individ', &
2542          tmpb*tmpa, &
2543          ( tmpveca(lptr_cl_aer(isize,itype,cw_phase))*tmpa, isize=1,nsize_aer(itype) )
2544
2545       tmpa = 1.0e-6 ; tmpb = 0.0
2546       do isize = 1, nsize_aer(itype)
2547          tmpb = tmpb + tmpveca(numptr_aer(isize,itype,cw_phase))
2548       end do
2549       write(lundiag,'(a,1p,10e12.4)') 'num_cw total, individ', &
2550          tmpb*tmpa, &
2551          ( tmpveca(numptr_aer(isize,itype,cw_phase))*tmpa, isize=1,nsize_aer(itype) )
2552
2553
2554
2555
2556       return
2557       end subroutine chem_cup_1d_diags_pt21
2558
2559
2560 !-------------------------------------------------------------------------------
2561       subroutine chem_cup_1d_diags_pt41( &
2562          lundiag_inp, iflagaa, &
2563          kts, kte, ktep1, p1st, num_chem, &
2564          ktau, grid_id, i, j, &
2565          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2566          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
2567          ai_phase, cw_phase, &
2568          massptr_aer, numptr_aer, &
2569          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
2570          rhodz, rhodzsum, &
2571          chem_av_old, chem_av_new, &
2572          zav_del_chem_actvbb, zav_del_chem_aqchbb, &
2573          zav_del_chem_resdbb, zav_del_chem_respbb, &
2574          zav_del_chem_totlbb, zav_del_chem_wetrem, &
2575          chem_name )
2576
2577       use module_configure, only:  &
2578          p_qc, &
2579          p_h2o2, p_hcl, p_hno3, p_nh3, p_so2, p_sulf
2580
2581 ! arguments
2582       integer, intent(in) :: lundiag_inp, iflagaa
2583       integer, intent(in) :: kts, kte, ktep1, p1st, num_chem
2584       integer, intent(in) :: ktau, grid_id, i, j
2585
2586       integer, intent(in) :: maxd_acomp, maxd_aphase, maxd_atype, maxd_asize
2587       integer, intent(in) :: &
2588          nphase_aer, ntype_aer, &
2589          ai_phase, cw_phase, &
2590          nsize_aer( maxd_atype ), &
2591          ncomp_aer( maxd_atype ), &
2592          massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
2593          numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
2594          lptr_so4_aer(maxd_asize, maxd_atype, maxd_aphase), &
2595          lptr_oin_aer(maxd_asize, maxd_atype, maxd_aphase), &
2596          lptr_no3_aer(maxd_asize, maxd_atype, maxd_aphase), &
2597          lptr_nh4_aer(maxd_asize, maxd_atype, maxd_aphase), &
2598          lptr_cl_aer(maxd_asize, maxd_atype, maxd_aphase)
2599
2600       real(r8), intent(in) :: rhodzsum
2601
2602       real(r8), intent(in), dimension(kts:kte) :: rhodz
2603
2604       real(r8), intent(in), dimension(kts:kte,1:num_chem) :: chem_av_old, chem_av_new
2605
2606       real(r8), intent(in), dimension(1:num_chem) :: &
2607          zav_del_chem_actvbb, zav_del_chem_aqchbb, &
2608          zav_del_chem_resdbb, zav_del_chem_respbb, &
2609          zav_del_chem_totlbb, zav_del_chem_wetrem
2610
2611       character(len=12), intent(in) :: chem_name(num_chem)
2612
2613 ! local variables
2614       integer :: isize, itype
2615       integer :: l, lundiag
2616       integer :: m
2617       real(r8) :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg, tmph, tmpi, tmpj, tmpk
2618       real(r8) :: tmpveca(num_chem)
2619
2620
2621       lundiag = lundiag_inp
2622 !     lundiag = 122
2623
2624       if (iflagaa <= 1) then
2625          write(lundiag,'(/a,i10,5i5)') &
2626             'chem_cup_1d_diags_pt41 - ktau, id, i, j =', ktau, grid_id, i, j
2627       else
2628          write(lundiag,'(/a,i10,5i5)') &
2629             'chem_cup_1d_diags_pt42 - ktau, id, i, j =', ktau, grid_id, i, j
2630       end if
2631
2632       do m = 1, 9
2633
2634       if (iflagaa <= 1) then
2635          if (m > 6) cycle
2636       else
2637          if (m < 7) cycle
2638       end if
2639
2640       if      (m == 1) then
2641          tmpveca(p1st:num_chem) = zav_del_chem_actvbb(p1st:num_chem)
2642          write(lundiag,'(/a)') 'inactive zav_del_chem_actvbb summary'
2643       else if (m == 2) then
2644          tmpveca(p1st:num_chem) = zav_del_chem_aqchbb(p1st:num_chem)
2645          write(lundiag,'(/a)') 'inactive zav_del_chem_aqchbb summary'
2646       else if (m == 3) then
2647          tmpveca(p1st:num_chem) = zav_del_chem_respbb(p1st:num_chem)
2648          write(lundiag,'(/a)') 'inactive zav_del_chem_respbb summary'
2649       else if (m == 4) then
2650          tmpveca(p1st:num_chem) = zav_del_chem_totlbb(p1st:num_chem)
2651          write(lundiag,'(/a)') 'inactive zav_del_chem_totlbb summary'
2652       else if (m == 5) then
2653          tmpveca(p1st:num_chem) = zav_del_chem_resdbb(p1st:num_chem)
2654          write(lundiag,'(/a)') 'inactive zav_del_chem_resdbb summary'
2655       else if (m == 6) then
2656          do l = p1st, num_chem
2657          tmpveca(l) = sum( rhodz(kts:kte)*chem_av_old(kts:kte,l) ) / rhodzsum
2658          end do
2659          write(lundiag,'(/a)') 'inactive zav_chem_av_old summary'
2660       else if (m == 7) then
2661          do l = p1st, num_chem
2662          tmpveca(l) = sum( rhodz(kts:kte)*chem_av_old(kts:kte,l) ) / rhodzsum
2663          end do
2664          write(lundiag,'(/a)') '***final zav_chem_av_old summary'
2665       else if (m == 8) then
2666          do l = p1st, num_chem
2667          tmpveca(l) = sum( rhodz(kts:kte)*chem_av_new(kts:kte,l) ) / rhodzsum
2668          end do
2669          write(lundiag,'(/a)') '***final zav_chem_av_new summary'
2670       else if (m == 9) then
2671          do l = p1st, num_chem
2672          tmpveca(l) = sum( rhodz(kts:kte)* &
2673             (chem_av_new(kts:kte,l)-chem_av_old(kts:kte,l)) ) / rhodzsum
2674          end do
2675          write(lundiag,'(/a)') '***final zav_chem_av_del summary'
2676       else
2677          cycle
2678       end if
2679
2680       write(lundiag,'(a,1p,10e12.4)') 'h2o2                 ', &
2681          tmpveca(p_h2o2)*1.0e3
2682
2683       tmpg = tmpveca(p_so2) + tmpveca(p_sulf)
2684       write(lundiag,'(a,1p,10e12.4)') 'so2+h2so4, individ   ', &
2685          tmpg*1.0e3, &
2686          tmpveca(p_so2)*1.0e3, &
2687          tmpveca(p_sulf)*1.0e3
2688       itype = 1
2689       tmpe = 28.966/96.0 ; tmpa = 0.0_r8 ; tmpc = 0.0_r8
2690       tmpi = 0.0_r8 ; tmpj = 0.0_r8
2691       do isize = 1, nsize_aer(itype)
2692          tmpa = tmpa + tmpveca(lptr_so4_aer(isize,itype,ai_phase))
2693          tmpc = tmpc + tmpveca(lptr_so4_aer(isize,itype,cw_phase))
2694          tmpi = tmpi + zav_del_chem_wetrem(lptr_so4_aer(isize,itype,ai_phase))
2695          tmpj = tmpj + zav_del_chem_wetrem(lptr_so4_aer(isize,itype,cw_phase))
2696       end do
2697       write(lundiag,'(a,1p,10e12.4)') 'so4_a  total, individ', &
2698          tmpa*tmpe, &
2699          ( tmpveca(lptr_so4_aer(isize,itype,ai_phase))*tmpe, isize=1,nsize_aer(itype) )
2700       write(lundiag,'(a,1p,10e12.4)') 'so4_cw total, individ', &
2701          tmpc*tmpe, &
2702          ( tmpveca(lptr_so4_aer(isize,itype,cw_phase))*tmpe, isize=1,nsize_aer(itype) )
2703       if (m >= 7) then
2704          tmph = tmpg*1.0e3 + (tmpa+tmpc)*tmpe
2705          write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2706          if (m <= 8) then
2707             write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2708          else
2709             tmpk = (tmpi+tmpj)*tmpe &
2710                  + (zav_del_chem_wetrem(p_so2)+zav_del_chem_wetrem(p_sulf))*1.0e3
2711             write(lundiag,'(a,1p,10e12.4)') '   all othr, wet, tot', (tmph-tmpk), tmpk, tmph
2712          end if
2713       end if
2714
2715       write(lundiag,'(a,1p,10e12.4)') 'nh3                  ', &
2716          tmpveca(p_nh3)*1.0e3
2717       tmpe = 28.966/18.0 ; tmpa = 0.0_r8 ; tmpc = 0.0_r8
2718       tmpi = 0.0_r8 ; tmpj = 0.0_r8
2719       do isize = 1, nsize_aer(itype)
2720          tmpa = tmpa + tmpveca(lptr_nh4_aer(isize,itype,ai_phase))
2721          tmpc = tmpc + tmpveca(lptr_nh4_aer(isize,itype,cw_phase))
2722          tmpi = tmpi + zav_del_chem_wetrem(lptr_nh4_aer(isize,itype,ai_phase))
2723          tmpj = tmpj + zav_del_chem_wetrem(lptr_nh4_aer(isize,itype,cw_phase))
2724       end do
2725       write(lundiag,'(a,1p,10e12.4)') 'nh4_a  total, individ', &
2726          tmpa*tmpe, &
2727          ( tmpveca(lptr_nh4_aer(isize,itype,ai_phase))*tmpe, isize=1,nsize_aer(itype) )
2728       write(lundiag,'(a,1p,10e12.4)') 'nh4_cw total, individ', &
2729          tmpc*tmpe, &
2730          ( tmpveca(lptr_nh4_aer(isize,itype,cw_phase))*tmpe, isize=1,nsize_aer(itype) )
2731       if (m >= 7) then
2732          tmph = tmpveca(p_nh3)*1.0e3 + (tmpa+tmpc)*tmpe
2733          write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2734          if (m <= 8) then
2735             write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2736          else
2737             tmpk = (tmpi+tmpj)*tmpe + zav_del_chem_wetrem(p_nh3)*1.0e3
2738             write(lundiag,'(a,1p,10e12.4)') '   all othr, wet, tot', (tmph-tmpk), tmpk, tmph
2739          end if
2740       end if
2741
2742       write(lundiag,'(a,1p,10e12.4)') 'hno3                 ', &
2743          tmpveca(p_hno3)*1.0e3
2744       tmpe = 28.966/62.0 ; tmpa = 0.0_r8 ; tmpc = 0.0_r8
2745       tmpi = 0.0_r8 ; tmpj = 0.0_r8
2746       do isize = 1, nsize_aer(itype)
2747          tmpa = tmpa + tmpveca(lptr_no3_aer(isize,itype,ai_phase))
2748          tmpc = tmpc + tmpveca(lptr_no3_aer(isize,itype,cw_phase))
2749          tmpi = tmpi + zav_del_chem_wetrem(lptr_no3_aer(isize,itype,ai_phase))
2750          tmpj = tmpj + zav_del_chem_wetrem(lptr_no3_aer(isize,itype,cw_phase))
2751       end do
2752       write(lundiag,'(a,1p,10e12.4)') 'no3_a  total, individ', &
2753          tmpa*tmpe, &
2754          ( tmpveca(lptr_no3_aer(isize,itype,ai_phase))*tmpe, isize=1,nsize_aer(itype) )
2755       write(lundiag,'(a,1p,10e12.4)') 'no3_cw total, individ', &
2756          tmpc*tmpe, &
2757          ( tmpveca(lptr_no3_aer(isize,itype,cw_phase))*tmpe, isize=1,nsize_aer(itype) )
2758       if (m >= 7) then
2759          tmph = tmpveca(p_hno3)*1.0e3 + (tmpa+tmpc)*tmpe
2760          write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2761          if (m <= 8) then
2762             write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2763          else
2764             tmpk = (tmpi+tmpj)*tmpe + zav_del_chem_wetrem(p_hno3)*1.0e3
2765             write(lundiag,'(a,1p,10e12.4)') '   all othr, wet, tot', (tmph-tmpk), tmpk, tmph
2766          end if
2767       end if
2768
2769       write(lundiag,'(a,1p,10e12.4)') 'hcl                  ', &
2770          tmpveca(p_hcl)*1.0e3
2771       tmpe = 28.966/35.5 ; tmpa = 0.0_r8 ; tmpc = 0.0_r8
2772       tmpi = 0.0_r8 ; tmpj = 0.0_r8
2773       do isize = 1, nsize_aer(itype)
2774          tmpa = tmpa + tmpveca(lptr_cl_aer(isize,itype,ai_phase))
2775          tmpc = tmpc + tmpveca(lptr_cl_aer(isize,itype,cw_phase))
2776          tmpi = tmpi + zav_del_chem_wetrem(lptr_cl_aer(isize,itype,ai_phase))
2777          tmpj = tmpj + zav_del_chem_wetrem(lptr_cl_aer(isize,itype,cw_phase))
2778       end do
2779       write(lundiag,'(a,1p,10e12.4)') 'cl_a   total, individ', &
2780          tmpa*tmpe, &
2781          ( tmpveca(lptr_cl_aer(isize,itype,ai_phase))*tmpe, isize=1,nsize_aer(itype) )
2782       write(lundiag,'(a,1p,10e12.4)') 'cl_cw  total, individ', &
2783          tmpc*tmpe, &
2784          ( tmpveca(lptr_cl_aer(isize,itype,cw_phase))*tmpe, isize=1,nsize_aer(itype) )
2785       if (m >= 7) then
2786          tmph = tmpveca(p_hcl)*1.0e3 + (tmpa+tmpc)*tmpe
2787          if (m <= 8) then
2788             write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2789          else
2790             tmpk = (tmpi+tmpj)*tmpe + zav_del_chem_wetrem(p_hcl)*1.0e3
2791             write(lundiag,'(a,1p,10e12.4)') '   all othr, wet, tot', (tmph-tmpk), tmpk, tmph
2792          end if
2793       end if
2794
2795       tmpe = 1.0e-6 ; tmpa = 0.0_r8 ; tmpc = 0.0_r8
2796        tmpi = 0.0_r8 ; tmpj = 0.0_r8
2797       do isize = 1, nsize_aer(itype)
2798          tmpa = tmpa + tmpveca(numptr_aer(isize,itype,ai_phase))
2799          tmpc = tmpc + tmpveca(numptr_aer(isize,itype,cw_phase))
2800          tmpi = tmpi + zav_del_chem_wetrem(numptr_aer(isize,itype,ai_phase))
2801          tmpj = tmpj + zav_del_chem_wetrem(numptr_aer(isize,itype,cw_phase))
2802       end do
2803       write(lundiag,'(a,1p,10e12.4)') 'num_a  total, individ', &
2804          tmpa*tmpe, &
2805          ( tmpveca(numptr_aer(isize,itype,ai_phase))*tmpe, isize=1,nsize_aer(itype) )
2806       write(lundiag,'(a,1p,10e12.4)') 'num_cw total, individ', &
2807          tmpc*tmpe, &
2808          ( tmpveca(numptr_aer(isize,itype,cw_phase))*tmpe, isize=1,nsize_aer(itype) )
2809       if (m >= 7) then
2810          tmph = (tmpa+tmpc)*tmpe
2811          if (m <= 8) then
2812             write(lundiag,'(a,1p,10e12.4)') '   all total         ', tmph
2813          else
2814             tmpk = (tmpi+tmpj)*tmpe
2815             write(lundiag,'(a,1p,10e12.4)') '   all othr, wet, tot', (tmph-tmpk), tmpk, tmph
2816          end if
2817       end if
2818
2819       end do ! m = 1, 6
2820
2821       return
2822       end subroutine chem_cup_1d_diags_pt41
2823
2824
2825 !-------------------------------------------------------------------------------
2826       subroutine chem_cup_1d_diags_pt71( &
2827          lundiag_inp, kcldbot, kcldtop, &
2828          kts, kte, ktep1, p1st, num_chem, &
2829          ktau, grid_id, i, j, &
2830          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2831          ncomp_aer, nphase_aer, nsize_aer, ntype_aer, &
2832          ai_phase, cw_phase, &
2833          massptr_aer, numptr_aer, &
2834          lptr_cl_aer, lptr_nh4_aer, lptr_no3_aer, lptr_oin_aer, lptr_so4_aer, &
2835          af_up, af_inact, &
2836          chem_av_new, chem_up, chem_inact_new, chem_incu, &
2837          chem_name )
2838
2839       use module_configure, only:  &
2840          p_qc, &
2841          p_h2o2, p_hcl, p_hno3, p_nh3, p_so2, p_sulf
2842
2843 ! arguments
2844       integer, intent(in) :: lundiag_inp, kcldbot, kcldtop
2845       integer, intent(in) :: kts, kte, ktep1, p1st, num_chem
2846       integer, intent(in) :: ktau, grid_id, i, j
2847
2848       integer, intent(in) :: maxd_acomp, maxd_aphase, maxd_atype, maxd_asize
2849       integer, intent(in) :: &
2850          nphase_aer, ntype_aer, &
2851          ai_phase, cw_phase, &
2852          nsize_aer( maxd_atype ), &
2853          ncomp_aer( maxd_atype ), &
2854          massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
2855          numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
2856          lptr_so4_aer(maxd_asize, maxd_atype, maxd_aphase), &
2857          lptr_oin_aer(maxd_asize, maxd_atype, maxd_aphase), &
2858          lptr_no3_aer(maxd_asize, maxd_atype, maxd_aphase), &
2859          lptr_nh4_aer(maxd_asize, maxd_atype, maxd_aphase), &
2860          lptr_cl_aer(maxd_asize, maxd_atype, maxd_aphase)
2861
2862       real(r8), intent(in), dimension(kts:kte) :: af_up, af_inact
2863
2864       real(r8), intent(in), dimension(kts:kte,1:num_chem) :: &
2865          chem_av_new, chem_inact_new, chem_incu
2866
2867       real(r8), intent(in), dimension(kts:ktep1,1:num_chem) :: chem_up
2868
2869       character(len=12), intent(in) :: chem_name(num_chem)
2870
2871 ! local variables
2872       integer :: k, n, n2
2873       integer :: lundiag
2874       real(r8) :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg, tmph, tmpi, tmpj, tmpk
2875       real(r8) :: tmpveca(num_chem)
2876
2877
2878       lundiag = lundiag_inp
2879 !     lundiag = 122
2880
2881       write(lundiag,'(/a,i10,5i5)') &
2882          'chem_cup_1d_diags_pt71 - ktau, id, i, j =', ktau, grid_id, i, j
2883
2884       n2 = (nsize_aer(1)+1)/2
2885
2886       write(lundiag,'(/a,i2.2,a,i2.2,a)') &
2887          'af_up, af_inact;   so4_a01:', n2, &
2888          ' gridav/up/inact/incu;   so4_cw01:', n2, ' ...'
2889       tmpa = (28.966/96.0)
2890       do k = kcldtop+1, kts, -1
2891          tmpveca = 0.0
2892          do n = 1, (nsize_aer(1)+1)/2
2893             tmpveca(1) = tmpveca(1) + chem_av_new(   k,lptr_so4_aer(n,1,1))
2894             tmpveca(2) = tmpveca(2) + chem_up(       k,lptr_so4_aer(n,1,1))
2895             tmpveca(3) = tmpveca(3) + chem_inact_new(k,lptr_so4_aer(n,1,1))
2896             tmpveca(4) = tmpveca(4) + chem_incu(     k,lptr_so4_aer(n,1,1))
2897             tmpveca(5) = tmpveca(5) + chem_av_new(   k,lptr_so4_aer(n,1,2))
2898             tmpveca(6) = tmpveca(6) + chem_up(       k,lptr_so4_aer(n,1,2))
2899             tmpveca(7) = tmpveca(7) + chem_inact_new(k,lptr_so4_aer(n,1,2))
2900             tmpveca(8) = tmpveca(8) + chem_incu(     k,lptr_so4_aer(n,1,2))
2901          end do
2902          write(lundiag,'(i3,2x,2f8.5,1p,2(2x,4e10.2))') &
2903             k, af_up(k), af_inact(k), tmpveca(1:8)*tmpa
2904       end do
2905
2906       write(lundiag,'(/a,i2.2,a,i2.2,a)') &
2907          'af_up, af_inact;   no3_a01:', n2, &
2908          ' gridav/up/inact/incu;   no3_cw01:', n2, ' ...'
2909       tmpa = (28.966/62.0)
2910       do k = kcldtop+1, kts, -1
2911          tmpveca = 0.0
2912          do n = 1, (nsize_aer(1)+1)/2
2913             tmpveca(1) = tmpveca(1) + chem_av_new(   k,lptr_no3_aer(n,1,1))
2914             tmpveca(2) = tmpveca(2) + chem_up(       k,lptr_no3_aer(n,1,1))
2915             tmpveca(3) = tmpveca(3) + chem_inact_new(k,lptr_no3_aer(n,1,1))
2916             tmpveca(4) = tmpveca(4) + chem_incu(     k,lptr_no3_aer(n,1,1))
2917             tmpveca(5) = tmpveca(5) + chem_av_new(   k,lptr_no3_aer(n,1,2))
2918             tmpveca(6) = tmpveca(6) + chem_up(       k,lptr_no3_aer(n,1,2))
2919             tmpveca(7) = tmpveca(7) + chem_inact_new(k,lptr_no3_aer(n,1,2))
2920             tmpveca(8) = tmpveca(8) + chem_incu(     k,lptr_no3_aer(n,1,2))
2921          end do
2922          write(lundiag,'(i3,2x,2f8.5,1p,2(2x,4e10.2))') &
2923             k, af_up(k), af_inact(k), tmpveca(1:8)*tmpa
2924       end do
2925
2926       write(lundiag,'(/2a)') &
2927          'af_up, af_inact;   so2', &
2928          ' gridav/up/inact/incu;   hno3 ...'
2929       tmpa = 1.0e3
2930       do k = kcldtop+1, kts, -1
2931          tmpveca = 0.0
2932          do n = 1, (nsize_aer(1)+1)/2
2933             tmpveca(1) = tmpveca(1) + chem_av_new(   k,p_so2)
2934             tmpveca(2) = tmpveca(2) + chem_up(       k,p_so2)
2935             tmpveca(3) = tmpveca(3) + chem_inact_new(k,p_so2)
2936             tmpveca(4) = tmpveca(4) + chem_incu(     k,p_so2)
2937             tmpveca(5) = tmpveca(5) + chem_av_new(   k,p_hno3)
2938             tmpveca(6) = tmpveca(6) + chem_up(       k,p_hno3)
2939             tmpveca(7) = tmpveca(7) + chem_inact_new(k,p_hno3)
2940             tmpveca(8) = tmpveca(8) + chem_incu(     k,p_hno3)
2941          end do
2942          write(lundiag,'(i3,2x,2f8.5,1p,2(2x,4e10.2))') &
2943             k, af_up(k), af_inact(k), tmpveca(1:8)*tmpa
2944       end do
2945
2946       return
2947       end subroutine chem_cup_1d_diags_pt71
2948
2949
2950 !-------------------------------------------------------------------------------
2951       subroutine chem_cup_activate_up( &
2952          lunerr, lundiag, idiagaa, &
2953          kts, kte, p1st, num_chem, &
2954          ktau, grid_id, i, j, k, iflagaa, &
2955          pcen, tcen, rhocen, qcw, &
2956          rhodz, af_up, wact, &
2957          tmp_mf_up, tmp_mfxchem_up, dchemdt_up_activa, &
2958          maxd_acomp, maxd_aphase, maxd_atype, maxd_asize, &
2959          ncomp_aer, nsize_aer, ntype_aer, &
2960          ai_phase, cw_phase, msectional, &
2961          massptr_aer, numptr_aer, &
2962          dlo_sect, dhi_sect, dens_aer, hygro_aer, sigmag_aer )
2963 !
2964 ! compute changes to updraft mixing ratios from aerosol activation
2965 !    as updraft air moves through layer k
2966 !
2967 ! when iflagaa = 1 or 2, do traditional cloud-base activation
2968 ! when iflagaa = 10,     do secondary activation
2969 !
2970       use module_mixactivate, only:  activate
2971
2972 ! arguments
2973       integer :: lunerr, lundiag, idiagaa
2974       integer :: kts, kte, p1st, num_chem
2975       integer :: ktau, grid_id, i, j, k, iflagaa
2976
2977       integer :: maxd_acomp, maxd_aphase, maxd_atype, maxd_asize
2978       integer :: &
2979          ntype_aer, &
2980          ai_phase, cw_phase, msectional, &
2981          nsize_aer( maxd_atype ), &
2982          ncomp_aer( maxd_atype ), &
2983          massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
2984          numptr_aer( maxd_asize, maxd_atype, maxd_aphase )
2985       real(r4) :: &
2986          dlo_sect( maxd_asize, maxd_atype ), &
2987          dhi_sect( maxd_asize, maxd_atype ), &
2988          dens_aer( maxd_acomp, maxd_atype ), &
2989          hygro_aer( maxd_acomp, maxd_atype ), &
2990          sigmag_aer(maxd_asize, maxd_atype)
2991
2992       real(r8) :: pcen, tcen, rhocen, qcw
2993       real(r8) :: rhodz, af_up, wact
2994       real(r8) :: tmp_mf_up
2995       real(r8) :: tmp_mfxchem_up(1:num_chem)
2996       real(r8) :: dchemdt_up_activa(kts:kte,1:num_chem)
2997
2998 ! local variables
2999       integer :: isize, itype
3000       integer :: la, lc, l2
3001       real(r8) :: tmpb, tmpc
3002       real(r8) :: tmpvol
3003       real(r8) :: tmp_chem_up(1:num_chem)
3004       real(r8), dimension( 1:maxd_asize, 1:maxd_atype ) :: &
3005                      fn, fm, hygro, numbr, volum
3006
3007 ! single precision variable used for calls to other wrf-chem routines
3008       real(r4) :: flux_fullact_sp
3009       real(r4) :: rhocen_sp, tcen_sp, wact_sp
3010       real(r4) :: smax_prescribed_sp
3011       real(r4), dimension( 1:maxd_asize, 1:maxd_atype ) :: &
3012                      fn_sp, fs_sp, fm_sp, fluxn_sp, fluxs_sp, fluxm_sp, &
3013                      hygro_sp, numbr_sp, volum_sp
3014
3015 !     if ( do_activa )
3016 ! get current mixing ratios in updraft
3017       tmp_chem_up(p1st:num_chem) = tmp_mfxchem_up(p1st:num_chem)/tmp_mf_up
3018
3019 ! calculate activation fractions using the abdul-razak and ghan
3020 !    parameterization (subr. activate)
3021 ! for this need aerosol number, volume, and weighted hygroscopicity for
3022 !    each aerosol type and size
3023       hygro(:,:) = 0.0_r8
3024       numbr(:,:) = 0.0_r8
3025       volum(:,:) = 0.0_r8
3026
3027       do itype = 1, ntype_aer
3028       do isize = 1, nsize_aer(itype)
3029          la = numptr_aer(isize,itype,ai_phase)
3030          numbr(isize,itype) = numbr(isize,itype) + max( 0.0_r8, tmp_chem_up(la) )
3031          do l2 = 1, ncomp_aer(itype)
3032             la = massptr_aer(l2,isize,itype,ai_phase)
3033             tmpvol = max( 0.0_r8, tmp_chem_up(la) ) / dens_aer(l2,itype)
3034             volum(isize,itype) = volum(isize,itype) + tmpvol
3035             hygro(isize,itype) = hygro(isize,itype) + tmpvol*hygro_aer(l2,itype)
3036          end do
3037       end do ! isize
3038       end do ! itype
3039
3040       do itype = 1, ntype_aer
3041       do isize = 1, nsize_aer(itype)
3042          hygro(isize,itype) = hygro(isize,itype) / max( 1.0e-35_r8, volum(isize,itype) )
3043          ! convert numbr from (#/kg) to (#/m3)
3044          numbr(isize,itype) = numbr(isize,itype)*rhocen
3045          ! convert volum to (m3/m3) -- need 1e-12 factor because
3046          !    (rhocen*chem)/dens_aer = [(ugaero/m3air)/(gaero/cm3aero)]
3047          volum(isize,itype) = volum(isize,itype)*rhocen*1.0e-12_r8
3048
3049          ! avoid zero numbr or volum
3050          tmpb = 1.0e-32_r8  ! 1e-32 m3/m3 volume ~= 1e-20 ug/m3 mass
3051          tmpc = (dlo_sect(isize,itype) + dhi_sect(isize,itype))*0.5e-2_r8  ! dcen in m
3052          tmpc = 1.91_r8 * tmpb / (tmpc**3)  ! number in #/m3 when volume = tmpb
3053          if ((volum(isize,itype) < tmpb) .or. (numbr(isize,itype) < tmpc)) then
3054             volum(isize,itype) = tmpb
3055             numbr(isize,itype) = tmpc
3056             hygro(isize,itype) = 0.3_r8
3057          end if
3058       end do ! isize
3059       end do ! itype
3060
3061 !     subr. activate( wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,  &
3062 !                     msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
3063 !                     na, volc, dlo_sect, dhi_sect, sigman, hygro, &
3064 !                     fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, &
3065 !                     grid_id, ktau, ii, jj, kk )
3066
3067 !     call activate( wact, 0.0, 0.0, 0.0, 1.0, tcen, rhocen,  &
3068 !                     msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
3069 !                     numbr, volum, dlo_sect, dhi_sect, sigmag_aer, hygro, &
3070 !                     fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, &
3071 !                     grid_id, ktau, i, j, k )
3072
3073       wact_sp = wact
3074       tcen_sp = tcen
3075       rhocen_sp = rhocen
3076       numbr_sp = numbr
3077       volum_sp = volum
3078       hygro_sp = hygro
3079
3080       if (iflagaa < 10) then
3081          call activate( &
3082                       wact_sp, 0.0, 0.0, 0.0, 1.0, tcen_sp, rhocen_sp,  &
3083                       msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
3084                       numbr_sp, volum_sp, dlo_sect, dhi_sect, sigmag_aer, hygro_sp, &
3085                       fn_sp, fs_sp, fm_sp, fluxn_sp, fluxs_sp, fluxm_sp, flux_fullact_sp, &
3086                       grid_id, ktau, i, j, k )
3087       else
3088          if (qcw < qcw_inup_smallaa) return
3089          smax_prescribed_sp = 0.001  ! 0.1% supersaturation
3090          call activate( &
3091                       wact_sp, 0.0, 0.0, 0.0, 1.0, tcen_sp, rhocen_sp,  &
3092                       msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
3093                       numbr_sp, volum_sp, dlo_sect, dhi_sect, sigmag_aer, hygro_sp, &
3094                       fn_sp, fs_sp, fm_sp, fluxn_sp, fluxs_sp, fluxm_sp, flux_fullact_sp, &
3095                       grid_id, ktau, i, j, k, smax_prescribed_sp )
3096       end if
3097
3098       fn = fn_sp
3099       fm = fm_sp
3100
3101       if (idiagaa > 0) then
3102          write(lundiag,'(/a,i10,5i5,1p,e10.2)/a') &
3103             'chem_cup_activate_up - ktau, id, i, j, k, iflagaa, wact', &
3104             ktau, grid_id, i, j, k, iflagaa, wact, &
3105             'chem_cup_activate_up - fn then fm'
3106          write(lundiag,'(8f11.7)') &
3107             ((fn(isize,itype), isize=1,nsize_aer(itype)), itype=1,ntype_aer)
3108          write(lundiag,'(8f11.7)') &
3109             ((fm(isize,itype), isize=1,nsize_aer(itype)), itype=1,ntype_aer)
3110 !        do itype = 1, ntype_aer
3111 !        do isize = 1, nsize_aer(itype)
3112 !           tmpb = volum(isize,itype) / max(numbr(isize,itype),1.0e-10)
3113 !           tmpb = ( tmpb*6.0/3.14159 )**(1.0/3.0)
3114 !           write(lundiag,'(a,2i5,1p,4e10.2)') 'itype, isize, diams', itype, isize, &
3115 !              dlo_sect(isize,itype), tmpb, dhi_sect(isize,itype)
3116 !        end do ! isize
3117 !        end do ! itype
3118       end if
3119
3120 ! apply activation fractions
3121       do itype = 1, ntype_aer
3122       do isize = 1, nsize_aer(itype)
3123          do l2 = 0, ncomp_aer(itype)
3124             if (l2 == 0) then
3125                la = numptr_aer(isize,itype,ai_phase)
3126                lc = numptr_aer(isize,itype,cw_phase)
3127             else
3128                la = massptr_aer(l2,isize,itype,ai_phase)
3129                lc = massptr_aer(l2,isize,itype,cw_phase)
3130             end if
3131             if ((la < p1st) .or. (la > num_chem)) cycle
3132             if ((lc < p1st) .or. (lc > num_chem)) cycle
3133
3134             if (l2 == 0) then
3135                tmpb = tmp_chem_up(la)*fn(isize,itype)
3136             else
3137                tmpb = tmp_chem_up(la)*fm(isize,itype)
3138             end if
3139             ! update chem_up
3140             tmp_chem_up(la) = tmp_chem_up(la) - tmpb
3141             tmp_chem_up(lc) = tmp_chem_up(lc) + tmpb
3142
3143             ! update mf_up*chem_up
3144             tmp_mfxchem_up(la) = tmp_chem_up(la)*tmp_mf_up
3145             tmp_mfxchem_up(lc) = tmp_chem_up(lc)*tmp_mf_up
3146
3147             ! increment dchemdt_up_activa
3148             tmpc = tmpb*tmp_mf_up/(rhodz*af_up)
3149             dchemdt_up_activa(k,la) = dchemdt_up_activa(k,la) - tmpc
3150             dchemdt_up_activa(k,lc) = dchemdt_up_activa(k,lc) + tmpc
3151
3152          end do ! l2
3153       end do ! isize
3154       end do ! itype
3155
3156       return
3157       end subroutine chem_cup_activate_up
3158
3159
3160 !-------------------------------------------------------------------------------
3161       subroutine chem_cup_aqchem( &
3162          config_flags, aer_mech_id, &
3163          lunerr, lundiag, idiagaa, &
3164          kts, kte, p1st, num_chem, &
3165          p_qc, num_moist, &
3166          ktau, grid_id, i, j, &
3167          iflagaa, ido_aqchem, &
3168          dt_aqchem, &
3169          pcen, tcen, rhocen, rhodz, qcw, ph_no2, &
3170          af_up, tmp_gas_aqfrac_up, tmp_mf_up, tmp_mfxchem_up, &
3171          dchemdt_up_aqchem, chem_inact )
3172 !
3173 ! compute changes to updraft mixing ratios from cloud chemistry
3174 !    as updraft air moves through layer k
3175 !
3176 ! when iflagaa >= kts, do cloud chemistry
3177 !    in updraft for the layer given by k=iflagaa
3178 ! when iflagaa < kts, do cloud chemistry
3179 !    in inactive cloud for the layers having ido_aqchem(k)>0
3180 !
3181       use module_configure, only:  grid_config_rec_type
3182       use module_mosaic_cloudchem, only:  mosaic_cloudchem_driver
3183
3184 ! arguments
3185       type(grid_config_rec_type), intent(in) :: config_flags
3186
3187       integer :: aer_mech_id
3188       integer :: lunerr, lundiag, idiagaa
3189       integer :: kts, kte, p1st, num_chem
3190       integer :: p_qc, num_moist
3191       integer :: ktau, grid_id, i, j
3192       integer :: iflagaa, ido_aqchem(kts:kte)
3193
3194       real(r8) :: dt_aqchem
3195       real(r8), dimension( kts:kte ) :: pcen, tcen, rhocen, rhodz, qcw, ph_no2
3196       real(r8) :: af_up
3197       real(r8) :: tmp_gas_aqfrac_up(1:num_chem)
3198       real(r8) :: tmp_mf_up
3199       real(r8) :: tmp_mfxchem_up(1:num_chem)
3200       real(r8) :: dchemdt_up_aqchem(kts:kte,1:num_chem)
3201       real(r8) :: chem_inact(kts:kte,1:num_chem)
3202
3203 ! local variables
3204       integer :: k, k2
3205       integer :: l
3206       real(r8) :: tmpb
3207       real(r8) :: tmp_chem_up(1:num_chem), tmp_chem_up_old(1:num_chem)
3208       real(r4), dimension( 1:1, kts:kte, 1:1 ) ::  tmp_ph_cw
3209
3210 ! single precision variable used for calls to other wrf-chem routines
3211       real(r4) :: dt_aqchem_sp
3212       real(r4), dimension( 1:1, kts:kte, 1:1 ) :: &
3213          tmp_alt_sp, tmp_cldfra_sp, tmp_p_sp, tmp_ph_no2_sp, tmp_rho_sp, tmp_t_sp
3214       real(r4), dimension( 1:1, kts:kte, 1:1, 1:num_chem ) :: &
3215          tmp_chem_sp, tmp_chem_old_sp, tmp_gas_aqfrac_sp
3216       real(r4), dimension( 1:1, kts:kte, 1:1, 1:num_moist ) :: &
3217          tmp_moist_sp
3218
3219
3220       if (aer_mech_id /= 3) return
3221
3222       tmp_gas_aqfrac_up = 0.0
3223
3224       tmp_chem_old_sp = 0.0
3225       tmp_chem_sp = 0.0
3226       tmp_cldfra_sp = 0.0
3227       tmp_moist_sp = 0.0
3228       tmp_ph_no2_sp = 0.0
3229       tmp_gas_aqfrac_sp = 0.0
3230       dt_aqchem_sp = dt_aqchem
3231
3232       if (iflagaa >= kts) then
3233          ! doing aqchem for 1 layer of updraft
3234          k = iflagaa
3235          ! first get current mixing ratios in updraft
3236          tmp_chem_up_old(p1st:num_chem) = tmp_mfxchem_up(p1st:num_chem)/tmp_mf_up
3237
3238          tmp_chem_old_sp(1,k,1,p1st:num_chem) = tmp_chem_up_old(p1st:num_chem)
3239          tmp_chem_sp(1,k,1,p1st:num_chem) = tmp_chem_old_sp(1,k,1,p1st:num_chem)
3240
3241          tmp_cldfra_sp(1,k,1) = 1.0
3242          tmp_moist_sp(1,k,1,p_qc) = qcw(k)
3243          tmp_ph_no2_sp(1,k,1) = ph_no2(k)
3244
3245       else
3246          ! doing aqchem for multiple levels of inactive cloud
3247          do k = kts, kte
3248             if (ido_aqchem(k) <= 0) cycle
3249             tmp_chem_old_sp(1,k,1,p1st:num_chem) = chem_inact(k,p1st:num_chem)
3250             tmp_chem_sp(1,k,1,p1st:num_chem) = tmp_chem_old_sp(1,k,1,p1st:num_chem)
3251
3252             tmp_cldfra_sp(1,k,1) = 1.0
3253             tmp_moist_sp(1,k,1,p_qc) = qcw(k)
3254             tmp_ph_no2_sp(1,k,1) = ph_no2(k)
3255          end do
3256       end if
3257
3258       do k = kts, kte
3259          tmp_p_sp(1,k,1) = pcen(k)
3260          tmp_t_sp(1,k,1) = tcen(k)
3261          tmp_rho_sp(1,k,1) = rhocen(k)
3262          tmp_alt_sp(1,k,1) = 1.0/rhocen(k)
3263       end do
3264
3265       if (aer_mech_id == 3) then
3266 !        subr. mosaic_cloudchem_driver(   &
3267 !           id, ktau, ktauc, dtstepc, config_flags,   &
3268 !           p_phy, t_phy, rho_phy, alt,   &
3269 !           cldfra, ph_no2,   &
3270 !           moist, chem,   &
3271 !           gas_aqfrac, numgas_aqfrac,   &
3272 !           ids,ide, jds,jde, kds,kde,   &
3273 !           ims,ime, jms,jme, kms,kme,   &
3274 !           its,ite, jts,jte, kts,kte )
3275
3276          call mosaic_cloudchem_driver(   &
3277             grid_id, ktau, ktau, dt_aqchem_sp, config_flags,   &
3278             tmp_p_sp, tmp_t_sp, tmp_rho_sp, tmp_alt_sp,   &
3279             tmp_cldfra_sp, tmp_ph_no2_sp,   &
3280             tmp_moist_sp, tmp_chem_sp, tmp_ph_cw,   &
3281             tmp_gas_aqfrac_sp, num_chem,   &
3282             1,1,  1,1,  kts,kte,  &
3283             1,1,  1,1,  kts,kte,  &
3284             1,1,  1,1,  kts,kte )
3285
3286 !     else
3287 !        call cloudchem routine for other aerosol treatments
3288       end if
3289
3290
3291       if (iflagaa >= kts) then
3292          k = iflagaa
3293          do l = p1st, num_chem
3294             tmp_gas_aqfrac_up(l) = tmp_gas_aqfrac_sp(1,k,1,l)
3295
3296             if (tmp_chem_sp(1,k,1,l) == tmp_chem_old_sp(1,k,1,l)) cycle
3297
3298             tmp_chem_up(l) = tmp_chem_sp(1,k,1,l)
3299             tmpb = tmp_chem_up(l) - tmp_chem_up_old(l)
3300
3301             ! update mf_up*chem_up
3302             tmp_mfxchem_up(l) = tmp_chem_up(l)*tmp_mf_up
3303             ! increment dchemdt_up_aqchem
3304             dchemdt_up_aqchem(k,l) = dchemdt_up_aqchem(k,l) &
3305                                    + tmpb*tmp_mf_up/(rhodz(k)*af_up)
3306          end do
3307
3308       else
3309          do k = kts, kte
3310             if (ido_aqchem(k) <= 0) cycle
3311             do l = p1st, num_chem
3312                if (tmp_chem_sp(1,k,1,l) == tmp_chem_old_sp(1,k,1,l)) cycle
3313                chem_inact(k,l) = tmp_chem_sp(1,k,1,l)
3314             end do
3315          end do
3316       end if
3317
3318       return
3319       end subroutine chem_cup_aqchem
3320
3321
3322 !-------------------------------------------------------------------------------
3323       subroutine chem_cup_check_adjust_inputs( &
3324          lunerr, lundiag, idiagaa, &
3325          kts, kte, ktep1, &
3326          ktau, grid_id, i, j, &
3327          ishall, do_dndraft, &
3328          kcldbot, kcldtop, kcldbotliq, &
3329          kupdrbot, kupdrtop, kdndrbot, kdndrtop, &
3330          iok, &
3331          tau_active, tau_inactive, &
3332          dz, zcen, zbnd, pcen, tcen, rhocen, &
3333          af_lscld, af_cucld, af_up, af_dn, &
3334          qcw_incu, qci_incu, &
3335          qcw_inup, qci_inup, &
3336          mf_up, mf_up_ent, mf_up_det, &
3337          mf_dn, mf_dn_ent, mf_dn_det )
3338 !
3339 ! apply "reality checks" to inputs
3340 ! to avoid problems in the chem_cup_1d algorithm
3341 !
3342
3343 ! arguments
3344       integer, intent(in) :: lunerr, lundiag, idiagaa
3345       integer, intent(in) :: kts, kte, ktep1
3346       integer, intent(in) :: ktau, grid_id, i, j
3347       integer, intent(in) :: ishall
3348       integer, intent(inout) :: kcldbot, kcldtop, kcldbotliq
3349       integer, intent(inout) :: kupdrbot, kupdrtop, kdndrbot, kdndrtop
3350       integer, intent(inout) :: iok
3351
3352       logical, intent(in) :: do_dndraft
3353
3354       real(r8), intent(in) :: tau_active, tau_inactive
3355
3356       real(r8), intent(in), dimension(kts:kte) :: &
3357          dz, zcen, pcen, tcen, rhocen, af_lscld
3358
3359       real(r8), intent(inout), dimension(kts:kte) :: &
3360          af_cucld, af_up, af_dn, &
3361          qcw_incu, qci_incu, qcw_inup, qci_inup, &
3362          mf_up_ent, mf_up_det, mf_dn_ent, mf_dn_det
3363
3364       real(r8), intent(in), dimension(kts:ktep1) :: &
3365          zbnd
3366
3367       real(r8), intent(inout), dimension(kts:ktep1) :: &
3368          mf_up, mf_dn
3369
3370 ! local variables
3371       integer :: k, ktmpa, ktmpb, ktmpc, ktmpd
3372
3373       real(r8) :: tmpa
3374
3375
3376       iok = -1
3377
3378 ! check for tau_active too small
3379       if (tau_active < tau_active_smallaa) then
3380          write(lunerr,'(2a,i10,3i5)') &
3381             'chem_cup_check_adjust_inputs - ', &
3382             'tau_active < tau_active_smallaa', ktau, grid_id, i, j
3383          return
3384       end if
3385
3386
3387 ! identify bottom and top of the updraft
3388       mf_up(kts) = 0.0_r8
3389       mf_up(ktep1) = 0.0_r8
3390       kupdrbot = kts-1
3391       kupdrtop = kts-1
3392       tmpa = 1.0e20
3393       do k = kts, kte
3394          if (mf_up(k) >= aw_up_smallaa*rhocen(k)) then
3395             if (kupdrbot < kts) kupdrbot = k-1
3396             kupdrtop = k
3397             tmpa = min( tmpa, mf_up(k) )
3398          else
3399             mf_up(k) = 0.0_r8
3400          end if
3401       end do
3402       if (kupdrbot < kts) then
3403          write(lunerr,'(2a,i10,3i5)') &
3404             'chem_cup_check_adjust_inputs - ', &
3405             'no mf_up > aw_up_smallaa*rho', ktau, grid_id, i, j
3406          return
3407       end if
3408
3409 ! force mf_up to have "reasonable" positive values throughout the updraft layers
3410       do k = kupdrbot+1, kupdrtop
3411          mf_up(k) = max( mf_up(k), tmpa )
3412       end do
3413
3414
3415       if ( do_dndraft ) then
3416 ! identify bottom and top of the dndraft
3417       mf_dn(kts) = 0.0_r8
3418       mf_dn(ktep1) = 0.0_r8
3419       kdndrbot = kts-1
3420       kdndrtop = kts-1
3421       tmpa = -1.0e20
3422       do k = kts, kte
3423          if (k > kupdrtop) then
3424             mf_dn(k) = 0.0_r8 ! do not allow downdraft to extend higher than updraft
3425          else if (mf_dn(k) <= -aw_up_smallaa*rhocen(k)) then
3426             if (kdndrbot < kts) kdndrbot = k-1
3427             kdndrtop = k
3428             tmpa = max( tmpa, mf_dn(k) )
3429          else
3430             mf_dn(k) = 0.0_r8
3431          end if
3432       end do
3433       if (kdndrbot < kts) then
3434          write(lunerr,'(2a,i10,3i5)') &
3435             'chem_cup_check_adjust_inputs - ', &
3436             'no mf_dn > aw_dn_smallaa*rho', ktau, grid_id, i, j
3437          return
3438       end if
3439
3440 ! force mf_dn to have "reasonable" negative values throughout the dndraft layers
3441       do k = kdndrbot+1, kdndrtop
3442          mf_dn(k) = min( mf_dn(k), tmpa )
3443       end do
3444       end if ! ( do_dndraft )
3445
3446
3447 ! identify layers with cloudwater and cloud-ice
3448       ktmpa = kts-1 ; ktmpb = kts-1 ; ktmpc = kts-1 ; ktmpd = kts-1
3449       do k = kts, kte
3450          if ((k < kupdrbot) .or. (k > kupdrtop)) then
3451             qcw_inup(k) = 0.0_r8
3452             qci_inup(k) = 0.0_r8
3453             qcw_incu(k) = 0.0_r8
3454             qci_incu(k) = 0.0_r8
3455          else
3456             if (qcw_inup(k) < qcw_inup_smallaa) then
3457                qcw_inup(k) = 0.0_r8
3458             else
3459                if (ktmpa < kts) ktmpa = k
3460                ktmpb = k
3461             end if
3462             if (qci_inup(k) < qci_inup_smallaa) then
3463                qci_inup(k) = 0.0_r8
3464             else
3465                if (ktmpc < kts) ktmpc = k
3466                ktmpd = k
3467             end if
3468
3469             if (qcw_incu(k) < qcw_incu_smallaa) then
3470                qcw_incu(k) = 0.0_r8
3471             end if
3472             if (qci_incu(k) < qci_incu_smallaa) then
3473                qci_incu(k) = 0.0_r8
3474             end if
3475          end if
3476       end do
3477       if ((ktmpa < kts) .and. (ktmpc < kts)) then
3478          write(lunerr,'(2a,i10,3i5)') &
3479             'chem_cup_check_adjust_inputs - ', &
3480             'no qcw/qci_inup > qcw/i_inup_smallaa', ktau, grid_id, i, j
3481          return
3482       end if
3483 ! at this point,  kupdrbot <= ktmpa <= ktmpb <= kupdrtop
3484 !    and/or       kupdrbot <= ktmpc <= ktmpd <= kupdrtop
3485
3486       kcldbotliq = ktmpa
3487       if (ktmpa < kts ) then
3488          ! in this case, ktmpc >= kts
3489          kcldbot = ktmpc
3490          kcldtop = ktmpd
3491          kcldbotliq = kts-1
3492       else if (ktmpc < kts ) then
3493          ! in this case, ktmpa >= kts
3494          kcldbot = ktmpa
3495          kcldtop = ktmpb
3496       else
3497          kcldbot = min( ktmpa, ktmpc )
3498          kcldtop = max( ktmpb, ktmpd )
3499       end if
3500
3501 ! check/adjust af_cucld and af_up
3502 !
3503 ! *** this should be done better
3504 !     check the ecpp code, which does something like
3505 !        af_up = max( af_up, mf_up/rho*wup_maxaa
3506       do k = kts, kte
3507          if ((k < kupdrbot) .or. (k > kupdrtop)) then
3508             af_up(k) = 0.0_r8
3509          else
3510             af_up(k) = max( af_up(k), af_up_smallaa )
3511             af_up(k) = min( af_up(k), af_up_maxaa )
3512          end if
3513       end do
3514
3515       do k = kts, kte
3516          if ((k < kcldbot) .or. (k > kcldtop)) then
3517             af_cucld(k) = 0.0_r8
3518          else
3519             af_cucld(k) = max( af_cucld(k), af_up(k), af_cucld_smallaa )
3520             af_cucld(k) = min( af_cucld(k), af_cucld_maxaa )
3521          end if
3522       end do
3523
3524 ! calculate mf_up_det and adjust mf_up_ent if needed
3525       do k = kts, kte
3526          if ((k < kupdrbot) .or. (k > kupdrtop)) then
3527             mf_up_ent(k) = 0.0_r8
3528             mf_up_det(k) = 0.0_r8
3529          else
3530             tmpa = mf_up(k+1) - mf_up(k)
3531             mf_up_det(k) = mf_up_ent(k) - tmpa
3532             if (mf_up_det(k) < 0.0_r8) then
3533                mf_up_ent(k) = tmpa
3534                mf_up_det(k) = 0.0_r8
3535             end if
3536         end if
3537       end do
3538
3539 ! calculate mf_dn_det and adjust mf_dn_ent if needed
3540       do k = kts, kte
3541          if ((k < kdndrbot) .or. (k > kdndrtop)) then
3542             mf_dn_ent(k) = 0.0_r8
3543             mf_dn_det(k) = 0.0_r8
3544          else
3545             tmpa = mf_dn(k+1) - mf_dn(k)
3546             mf_dn_det(k) = mf_dn_ent(k) - tmpa
3547             if (mf_dn_det(k) < 0.0_r8) then
3548                mf_dn_ent(k) = tmpa
3549                mf_dn_det(k) = 0.0_r8
3550             end if
3551         end if
3552       end do
3553
3554       iok = 0
3555       return
3556       end subroutine chem_cup_check_adjust_inputs
3557
3558
3559 !-----------------------------------------------------------------
3560       end module module_chem_cup