chem/module_sorgam_vbs_cloudchem.F

   1 !**********************************************************************************
   2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
   3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
   4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
   5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
   6 !
   7 ! MOSAIC module: see module_mosaic_driver.F for references and terms of use
   8 !**********************************************************************************
   9
  10         module module_sorgam_vbs_cloudchem
  11
  12
  13
  14         integer, parameter :: l_so4_aqyy = 1
  15         integer, parameter :: l_no3_aqyy = 2
  16         integer, parameter :: l_cl_aqyy  = 3
  17         integer, parameter :: l_nh4_aqyy = 4
  18         integer, parameter :: l_na_aqyy  = 5
  19         integer, parameter :: l_oin_aqyy = 6
  20         integer, parameter :: l_bc_aqyy  = 7
  21         integer, parameter :: l_oc_aqyy  = 8
  22
  23         integer, parameter :: nyyy = 8
  24
  25 ! "negligible volume" = (1 #/kg ~= 1e-6 #/cm3) * ((dp=1e-6 cm)**3 * pi/6)
  26         real, parameter :: smallvolaa = 0.5e-18
  27
  28
  29
  30         contains
  31
  32
  33
  34 !-----------------------------------------------------------------------
  35         subroutine sorgam_vbs_cloudchem_driver(   &
  36             id, ktau, ktauc, dtstepc, config_flags,   &
  37             p_phy, t_phy, rho_phy, alt,   &
  38             cldfra, ph_no2,   &
  39             moist, chem,   &
  40             gas_aqfrac, numgas_aqfrac,   &
  41             ids,ide, jds,jde, kds,kde,   &
  42             ims,ime, jms,jme, kms,kme,   &
  43             its,ite, jts,jte, kts,kte )
  44
  45         use module_state_description, only:   &
  46                 num_moist, num_chem, p_qc
  47
  48         use module_configure, only:  grid_config_rec_type
  49
  50         use module_data_sorgam_vbs, only:  cw_phase, nphase_aer
  51
  52
  53         implicit none
  54
  55 !   subr arguments
  56         integer, intent(in) ::   &
  57                 id, ktau, ktauc,   &
  58                 numgas_aqfrac,   &
  59                 ids, ide, jds, jde, kds, kde,   &
  60                 ims, ime, jms, jme, kms, kme,   &
  61                 its, ite, jts, jte, kts, kte
  62 !   id - domain index
  63 !   ktau - time step number
  64 !   ktauc - gas and aerosol chemistry time step number
  65 !   numgas_aqfrac - last dimension of gas_aqfrac
  66
  67 !   [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for 'domain'
  68 !   [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for 'memory'
  69 !       Most arrays that are arguments to chem_driver
  70 !       are dimensioned with these spatial indices.
  71 !   [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for 'tile'
  72 !       chem_driver and routines under it do calculations
  73 !       over these spatial indices.
  74
  75         type(grid_config_rec_type), intent(in) :: config_flags
  76 !   config_flags - configuration and control parameters
  77
  78         real, intent(in) ::   &
  79             dtstepc
  80 !   dtstepc - time step for gas and aerosol chemistry(s)
  81
  82         real, intent(in),   &
  83                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  84                 p_phy, t_phy, rho_phy, alt, cldfra, ph_no2
  85 !   p_phy - air pressure (Pa)
  86 !   t_phy - temperature (K)
  87 !   rho_phy - moist air density (kg/m^3)
  88 !   alt - dry air specific volume (m^3/kg)
  89 !   cldfra - cloud fractional area (0-1)
  90 !   ph_no2 - no2 photolysis rate (1/min)
  91
  92         real, intent(in),   &
  93                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
  94                 moist
  95 !   moist - mixing ratios of moisture species (water vapor,
  96 !       cloud water, ...) (kg/kg for mass species, #/kg for number species)
  97
  98         real, intent(inout),   &
  99                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
 100                 chem
 101 !   chem - mixing ratios of trace gas and aerosol species (ppm for gases,
 102 !       ug/kg for aerosol mass species, #/kg for aerosol number species)
 103
 104         real, intent(inout),   &
 105                 dimension( ims:ime, kms:kme, jms:jme, numgas_aqfrac ) :: &
 106                 gas_aqfrac
 107 !   gas_aqfrac - fraction (0-1) of gas that is dissolved in cloud water
 108
 109
 110 !   local variables
 111         integer :: it, jt, kt, l
 112         integer :: icase
 113         integer :: igaschem_onoff, iphotol_onoff, iradical_onoff
 114
 115         real :: rbox(num_chem)
 116         real :: gas_aqfrac_box(numgas_aqfrac)
 117         real :: ph_aq_box
 118         real, parameter :: qcldwtr_cutoff = 1.0e-6
 119         real :: qcldwtr
 120
 121
 122 !   check that cw_phase is active
 123         if ((cw_phase .le. 0) .or. (cw_phase .gt. nphase_aer)) then
 124             print *, '*** sorgam_vbs_cloudchem_driver - cw_phase not active'
 125             return
 126         end if
 127
 128         print 93010, 'entering sorgam_vbs_cloudchem_driver - ktau =', ktau
 129
 130         icase = 0
 131
 132 !   iphotol_onoff = 1 if photolysis rate calcs are on; 0 if off
 133         iphotol_onoff = 0
 134         if (config_flags%phot_opt .gt. 0) iphotol_onoff = 1
 135 !   igaschem_onoff = 1 if gas-phase chemistry is on; 0 if off
 136         igaschem_onoff = 0
 137         if (config_flags%gaschem_onoff .gt. 0) igaschem_onoff = 1
 138
 139 !   iradical_onoff turns aqueous radical chemistry on/off
 140 !   set iradical_onoff=0 if either photolysis or gas-phase chem are off
 141         if ((igaschem_onoff .le. 0) .or. (iphotol_onoff  .le. 0)) then
 142             iradical_onoff = 0
 143         else
 144             iradical_onoff = 1
 145         end if
 146 !   following line turns aqueous radical chem off unconditionally
 147         iradical_onoff = 0
 148
 149
 150         do 3920 jt = jts, jte
 151         do 3910 it = its, ite
 152
 153         do 3800 kt = kts, kte
 154
 155         qcldwtr = moist(it,kt,jt,p_qc)
 156         if (qcldwtr .le. qcldwtr_cutoff) goto 3800
 157
 158
 159         icase = icase + 1
 160
 161 !   detailed dump for debugging
 162         if (ktau .eq. -13579) then
 163 !       if ((ktau .eq. 30) .and. (it .eq. 23) .and.   &
 164 !             (jt .eq.  1) .and. (kt .eq. 11)) then
 165           call sorgam_cloudchem_dumpaa(   &
 166             id, ktau, ktauc, dtstepc, config_flags,   &
 167             p_phy, t_phy, rho_phy, alt,   &
 168             cldfra, ph_no2,   &
 169             moist, chem,   &
 170             gas_aqfrac, numgas_aqfrac,   &
 171             ids,ide, jds,jde, kds,kde,   &
 172             ims,ime, jms,jme, kms,kme,   &
 173             its,ite, jts,jte, kts,kte,   &
 174             qcldwtr_cutoff,   &
 175             it, jt, kt )
 176         end if
 177
 178 !   map from wrf-chem 3d arrays to pegasus clm & sub arrays
 179         rbox(1:num_chem) = chem(it,kt,jt,1:num_chem)
 180
 181 !   make call '1box' cloudchem routine
 182 !       print 93010, 'calling sorgam_cloudchem_1 at ijk =', it, jt, kt
 183         call sorgam_cloudchem_1box(   &
 184             id, ktau, ktauc, dtstepc,   &
 185             iphotol_onoff, iradical_onoff,   &
 186             ph_no2(it,kt,jt),   &
 187             ph_aq_box, gas_aqfrac_box,   &
 188             numgas_aqfrac, it, jt, kt, icase,   &
 189             rbox, qcldwtr,   &
 190             t_phy(it,kt,jt), p_phy(it,kt,jt), rho_phy(it,kt,jt),&
 191             config_flags)
 192
 193 !   map back to wrf-chem 3d arrays
 194         chem(it,kt,jt,1:num_chem) = rbox(1:num_chem)
 195         gas_aqfrac(it,kt,jt,:) = gas_aqfrac_box(:)
 196
 197
 198 3800    continue
 199
 200 3910    continue
 201 3920    continue
 202
 203         print 93010, 'leaving  sorgam_vbs_cloudchem_driver - ktau =', ktau, icase
 204 93010   format( a, 8(1x,i6) )
 205
 206         return
 207         end subroutine sorgam_vbs_cloudchem_driver
 208
 209
 210
 211 !-----------------------------------------------------------------------
 212         subroutine sorgam_cloudchem_1box(   &
 213             id, ktau, ktauc, dtstepc,   &
 214             iphotol_onoff, iradical_onoff,   &
 215             photol_no2_box,   &
 216             ph_aq_box, gas_aqfrac_box,   &
 217             numgas_aqfrac, it, jt, kt, icase,   &
 218             rbox, qcw_box, temp_box, pres_box, rho_box, &
 219             config_flags )
 220
 221         use module_configure, only: grid_config_rec_type
 222
 223         use module_state_description, only:   &
 224                 num_moist, num_chem
 225
 226         use module_data_sorgam_vbs, only:   &
 227                 msectional, maxd_asize, maxd_atype,   &
 228                 cw_phase, nsize_aer, ntype_aer, do_cloudchem_aer,   &
 229                 lptr_so4_aer, lptr_no3_aer, lptr_nh4_aer,   &
 230                 lptr_orgpa_aer, lptr_ec_aer, lptr_p25_aer, &
 231                 lptr_cl_aer, lptr_na_aer
 232
 233         use module_data_cmu_bulkaqchem, only:   &
 234                 meqn1max
 235
 236
 237         implicit none
 238
 239 !   subr arguments
 240
 241         type(grid_config_rec_type), intent(in) :: config_flags
 242
 243         integer, intent(in) ::   &
 244                 id, ktau, ktauc,   &
 245                 numgas_aqfrac, it, jt, kt,   &
 246                 icase, iphotol_onoff, iradical_onoff
 247
 248         real, intent(in) ::   &
 249             dtstepc, photol_no2_box,   &
 250             qcw_box,   &                ! cloud water (kg/kg)
 251             temp_box,   &               ! air temp (K)
 252             pres_box,   &               ! air pres (Pa)
 253             rho_box                     ! air dens (kg/m3)
 254
 255         real, intent(inout) :: ph_aq_box
 256
 257         real, intent(inout), dimension( num_chem ) :: rbox
 258 !   rbox has same units as chem [gas = ppmv, AP mass = ug/kg, AP number = #/kg]
 259
 260         real, intent(inout), dimension( numgas_aqfrac ) :: gas_aqfrac_box
 261
 262 !   local variables
 263         integer :: iphase
 264         integer :: icase_in, idecomp_hmsa_hso5,   &
 265                 iradical_in, istat_aqop
 266
 267         integer :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
 268
 269         real :: co2_mixrat_in
 270         real :: ph_cmuaq_cur
 271         real :: photol_no2_in
 272         real :: xprescribe_ph
 273
 274         real :: yaq_beg(meqn1max), yaq_end(meqn1max)
 275         real :: rbox_sv1(num_chem)
 276         real :: rbulk_cwaer(nyyy,2)
 277
 278         real, dimension( maxd_asize, maxd_atype ) :: fr_partit_cw
 279         real, dimension( 2, 3 ) :: xvol_old
 280
 281
 282 !
 283 !   set the lptr_yyy_cwaer
 284 !
 285         iphase = cw_phase
 286         lptr_yyy_cwaer(:,:,l_so4_aqyy) = lptr_so4_aer(:,:,iphase)
 287         lptr_yyy_cwaer(:,:,l_no3_aqyy) = lptr_no3_aer(:,:,iphase)
 288         lptr_yyy_cwaer(:,:,l_nh4_aqyy) = lptr_nh4_aer(:,:,iphase)
 289         lptr_yyy_cwaer(:,:,l_oin_aqyy) = lptr_p25_aer(:,:,iphase)
 290         lptr_yyy_cwaer(:,:,l_bc_aqyy ) = lptr_ec_aer( :,:,iphase)
 291         lptr_yyy_cwaer(:,:,l_oc_aqyy ) = lptr_orgpa_aer( :,:,iphase)
 292 ! na and cl added
 293         lptr_yyy_cwaer(:,:,l_cl_aqyy ) = lptr_cl_aer(:,:,iphase)
 294         lptr_yyy_cwaer(:,:,l_na_aqyy ) = lptr_na_aer(:,:,iphase)
 295 !       lptr_yyy_cwaer(:,:,l_cl_aqyy ) = -999888777
 296 !       lptr_yyy_cwaer(:,:,l_na_aqyy ) = -999888777
 297
 298 !
 299 !
 300 !   do bulk cloud-water chemistry
 301 !
 302 !
 303         icase_in = icase
 304         iradical_in = 1
 305         idecomp_hmsa_hso5 = 1
 306
 307         co2_mixrat_in = 350.0
 308
 309         photol_no2_in = photol_no2_box
 310
 311 !   when xprescribe_ph >= 0, ph is forced to its value
 312         xprescribe_ph = -1.0e31
 313
 314 !   turn off aqueous phase photolytic and radical chemistry
 315 !   if either of the iphotol_onoff and iradical_onoff flags are 0
 316         if ((iphotol_onoff .le. 0) .or. (iradical_onoff .le. 0)) then
 317             photol_no2_in = 0.0
 318             iradical_in = 0
 319         end if
 320
 321 #if defined ( ccboxtest_box_testing_active)
 322 !   following is for off-line box testing only
 323         call ccboxtest_extra_args_aa( 'get',   &
 324                 co2_mixrat_in, iradical_in,   &
 325                 idecomp_hmsa_hso5, icase_in,   &
 326                 xprescribe_ph )
 327 #endif
 328
 329         rbox_sv1(:) = rbox(:)
 330         gas_aqfrac_box(:) = 0.0
 331
 332
 333 !   make call to interface_to_aqoperator1
 334         call sorgam_interface_to_aqoperator1(   &
 335             istat_aqop,   &
 336             dtstepc,   &
 337             rbox, gas_aqfrac_box,   &
 338             qcw_box, temp_box, pres_box, rho_box,   &
 339             rbulk_cwaer, lptr_yyy_cwaer,   &
 340             co2_mixrat_in, photol_no2_in, xprescribe_ph,   &
 341             iradical_in, idecomp_hmsa_hso5,   &
 342             yaq_beg, yaq_end, ph_cmuaq_cur,   &
 343             numgas_aqfrac, id, it, jt, kt, ktau, icase_in, &
 344             config_flags )
 345
 346         ph_aq_box = ph_cmuaq_cur
 347
 348
 349 #if defined ( ccboxtest_box_testing_active)
 350 !   following is for off-line box testing only
 351         call ccboxtest_extra_args_bb( 'put',   &
 352                 yaq_beg, yaq_end, ph_cmuaq_cur )
 353 #endif
 354
 355
 356 !
 357 !
 358 !   calculate fraction of cloud-water associated with each activated aerosol bin
 359 !
 360 !
 361         call sorgam_partition_cldwtr(   &
 362             rbox, fr_partit_cw, xvol_old,   &
 363             id, it, jt, kt, icase_in )
 364
 365 !
 366 !
 367 !   distribute changes in bulk cloud-water composition among size bins
 368 !
 369 !
 370         call sorgam_distribute_bulk_changes(   &
 371             rbox, rbox_sv1, fr_partit_cw,   &
 372             rbulk_cwaer, lptr_yyy_cwaer,   &
 373             id, it, jt, kt, icase_in )
 374
 375
 376 !
 377 !   do move-sections
 378 !
 379         if (msectional .lt. 1000000000) then
 380             call sorgam_cloudchem_apply_mode_transfer(   &
 381                 rbox, rbox_sv1, xvol_old,   &
 382                 id, it, jt, kt, icase_in )
 383         end if
 384
 385
 386
 387         return
 388         end subroutine sorgam_cloudchem_1box
 389
 390
 391
 392 !-----------------------------------------------------------------------
 393         subroutine sorgam_interface_to_aqoperator1(   &
 394             istat_aqop,   &
 395             dtstepc,   &
 396             rbox, gas_aqfrac_box,   &
 397             qcw_box, temp_box, pres_box, rho_box,   &
 398             rbulk_cwaer, lptr_yyy_cwaer,   &
 399             co2_mixrat_in, photol_no2_in, xprescribe_ph,   &
 400             iradical_in, idecomp_hmsa_hso5,   &
 401             yaq_beg, yaq_end, ph_cmuaq_cur,   &
 402             numgas_aqfrac, id, it, jt, kt, ktau, icase, &
 403             config_flags )
 404
 405         use module_configure, only: grid_config_rec_type
 406
 407         use module_state_description, only:   &
 408                 num_chem, param_first_scalar, p_qc,   &
 409                 p_nh3, p_hno3, p_hcl, p_sulf, p_hcho,   &
 410                 p_ora1, p_so2, p_h2o2, p_o3, p_ho,   &
 411                 p_ho2, p_no3, p_no, p_no2, p_hono,   &
 412                 p_pan, p_ch3o2, p_ch3oh, p_op1, &
 413                 p_form, p_facd, p_oh, p_meo2, p_meoh, p_mepx, &
 414                 CB05_SORG_VBS_AQ_KPP
 415
 416         use module_data_cmu_bulkaqchem, only:   &
 417                 meqn1max, naers, ngas,   &
 418                 na4, naa, nac, nae, nah, nahmsa, nahso5,   &
 419                 nan, nao, nar, nas, naw,   &
 420                 ng4, nga, ngc, ngch3co3h, ngch3o2, ngch3o2h, ngch3oh,   &
 421                 ngh2o2, nghcho, nghcooh, nghno2, ngho2,   &
 422                 ngn, ngno, ngno2, ngno3, ngo3, ngoh, ngpan, ngso2
 423
 424         use module_cmu_bulkaqchem, only:  aqoperator1
 425
 426         use module_data_sorgam_vbs, only:   &
 427                 maxd_asize, maxd_atype,   &
 428                 cw_phase, nsize_aer, ntype_aer, do_cloudchem_aer,   &
 429                 lptr_so4_aer, lptr_no3_aer, lptr_nh4_aer,   &
 430                 lptr_orgpa_aer, lptr_ec_aer, lptr_p25_aer, &
 431                 lptr_cl_aer, lptr_na_aer
 432
 433
 434         implicit none
 435
 436 !   subr arguments
 437
 438         type(grid_config_rec_type), intent(in) :: config_flags
 439
 440         integer, intent(in) ::   &
 441                 iradical_in, idecomp_hmsa_hso5,   &
 442                 numgas_aqfrac, id, it, jt, kt, ktau, icase
 443         integer, intent(inout) ::   &
 444                 istat_aqop
 445
 446         integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
 447
 448         real, intent(in) ::   &
 449             dtstepc, co2_mixrat_in,   &
 450             photol_no2_in, xprescribe_ph,   &
 451             qcw_box, temp_box, pres_box, rho_box
 452
 453         real, intent(inout) :: ph_cmuaq_cur
 454
 455         real, intent(inout), dimension( num_chem ) :: rbox      ! ppm or ug/kg
 456
 457         real, intent(inout), dimension( numgas_aqfrac ) :: gas_aqfrac_box
 458
 459         real, intent(inout), dimension( nyyy, 2 ) :: rbulk_cwaer
 460
 461         real, intent(inout), dimension( meqn1max ) :: yaq_beg, yaq_end
 462
 463
 464 !   local variables
 465         integer :: i, iphase, isize, itype
 466         integer :: iaq, istat_fatal, istat_warn
 467         integer :: l, lyyy
 468         integer :: p1st
 469 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
 470         integer :: lunxx
 471 #endif
 472
 473         real, parameter :: eps=0.622   ! (mw h2o)/(mw air)
 474
 475
 476         real :: cair_moleperm3
 477         real :: dum, dumb
 478         real :: factgas, factlwc, factpatm, factphoto
 479         real :: factaerbc, factaercl, factaerna, factaernh4,   &
 480                 factaerno3, factaeroc, factaeroin, factaerso4
 481         real :: lwc
 482         real :: p_atm, photo_in
 483         real :: rh
 484         real :: temp, tstep_beg_sec, tstep_end_sec
 485         real :: totsulf_beg, totsulf_end
 486         real :: gas(ngas), aerosol(naers)
 487         real :: gas_aqfrac_cmu(ngas)
 488
 489         double precision tstep_beg_sec_dp, tstep_end_sec_dp,   &
 490           temp_dp, p_atm_dp, lwc_dp, rh_dp,   &
 491           co2_mixrat_in_dp, photo_in_dp, ph_cmuaq_cur_dp,   &
 492           xprescribe_ph_dp
 493         double precision gas_dp(ngas), gas_aqfrac_cmu_dp(ngas),   &
 494           aerosol_dp(naers), yaq_beg_dp(meqn1max), yaq_end_dp(meqn1max)
 495
 496
 497
 498         p1st = param_first_scalar
 499
 500 !
 501 !   units conversion factors
 502 !   'cmuaq-bulk' value = pegasus value X factor
 503 !
 504 !   [pres in atmospheres] = [pres in pascals] * factpatm
 505         factpatm = 1.0/1.01325e5
 506 !   [cldwtr in g-h2o/m3-air] = [cldwtr in kg-h2o/kg-air] * factlwc
 507         factlwc = 1.0e3*rho_box
 508 !   [aq photolysis rate scaling factor in --] = [jno2 in 1/min] * factphoto
 509         factphoto = 1.6
 510
 511 !   [gas in ppm] = [gas in ppm] * factgas
 512         factgas = 1.0
 513
 514 !   [aerosol in ug/m3-air] = [aerosol in ug/kg-air] * factaer
 515         dum = rho_box
 516         factaerso4   = dum
 517         factaerno3   = dum
 518         factaercl    = dum
 519         factaernh4   = dum
 520         factaerna    = dum
 521         factaeroin   = dum
 522         factaeroc    = dum
 523         factaerbc    = dum
 524
 525 #if defined ( ccboxtest_box_testing_active)
 526 !   If aboxtest_units_convert=10, turn off units conversions both here
 527 !   and in module_mosaic.  This is for testing, to allow exact agreements.
 528         if (aboxtest_units_convert .eq. 10) then
 529             factpatm = 1.0
 530             factlwc = 1.0
 531             factphoto = 1.0
 532             factgas = 1.0
 533             factaerso4   = 1.0
 534             factaerno3   = 1.0
 535             factaercl    = 1.0
 536             factaernh4   = 1.0
 537             factaerna    = 1.0
 538             factaeroin   = 1.0
 539             factaeroc    = 1.0
 540             factaerbc    = 1.0
 541         end if
 542 #endif
 543
 544 !
 545 !   map from rbox to gas,aerosol
 546 !
 547         temp = temp_box
 548
 549         lwc = qcw_box * factlwc
 550         p_atm = pres_box * factpatm
 551
 552 !   rce 2005-jul-11 - set p_atm so that cmu code's cair will match cairclm
 553 !       p_atm = cairclm(kpeg)*1.0e3*0.082058e0*temp
 554 !   for made-sorgam, set p_atm so that cmu code's (cair*28.966e3)
 555 !       will match rho_box
 556         p_atm = (rho_box/28.966)*0.082058e0*temp
 557
 558         photo_in = photol_no2_in * factphoto
 559
 560         rh = 1.0
 561         iaq = 1
 562
 563         tstep_beg_sec = 0.0
 564         tstep_end_sec = dtstepc
 565
 566 !   map gases and convert to ppm
 567         gas(:) = 0.0
 568
 569         if (p_nh3   >= p1st) gas(nga     ) = rbox(p_nh3  )*factgas
 570         if (p_hno3  >= p1st) gas(ngn     ) = rbox(p_hno3 )*factgas
 571         if (p_hcl   >= p1st) gas(ngc     ) = rbox(p_hcl  )*factgas
 572         if (p_sulf  >= p1st) gas(ng4     ) = rbox(p_sulf )*factgas
 573
 574 !       if (p_hcho  >= p1st) gas(nghcho  ) = rbox(p_hcho )*factgas
 575 !       if (p_ora1  >= p1st) gas(nghcooh ) = rbox(p_ora1 )*factgas
 576         if (p_so2   >= p1st) gas(ngso2   ) = rbox(p_so2  )*factgas
 577         if (p_h2o2  >= p1st) gas(ngh2o2  ) = rbox(p_h2o2 )*factgas
 578         if (p_o3    >= p1st) gas(ngo3    ) = rbox(p_o3   )*factgas
 579 !       if (p_ho    >= p1st) gas(ngoh    ) = rbox(p_ho   )*factgas
 580         if (p_ho2   >= p1st) gas(ngho2   ) = rbox(p_ho2  )*factgas
 581         if (p_no3   >= p1st) gas(ngno3   ) = rbox(p_no3  )*factgas
 582
 583         if (p_no    >= p1st) gas(ngno    ) = rbox(p_no   )*factgas
 584         if (p_no2   >= p1st) gas(ngno2   ) = rbox(p_no2  )*factgas
 585         if (p_hono  >= p1st) gas(nghno2  ) = rbox(p_hono )*factgas
 586         if (p_pan   >= p1st) gas(ngpan   ) = rbox(p_pan  )*factgas
 587 !       if (p_ch3o2 >= p1st) gas(ngch3o2 ) = rbox(p_ch3o2)*factgas
 588 !       if (p_ch3oh >= p1st) gas(ngch3oh ) = rbox(p_ch3oh)*factgas
 589 !       if (p_op1   >= p1st) gas(ngch3o2h) = rbox(p_op1  )*factgas
 590
 591 ! CB05 case added below
 592         if ((config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP)) then
 593
 594           if (p_form  >= p1st) gas(nghcho  ) = rbox(p_form )*factgas
 595           if (p_facd  >= p1st) gas(nghcooh ) = rbox(p_facd )*factgas
 596           if (p_oh    >= p1st) gas(ngoh    ) = rbox(p_oh   )*factgas
 597           if (p_meo2  >= p1st) gas(ngch3o2 ) = rbox(p_meo2 )*factgas
 598           if (p_meoh  >= p1st) gas(ngch3oh ) = rbox(p_meoh )*factgas
 599           if (p_mepx  >= p1st) gas(ngch3o2h) = rbox(p_mepx )*factgas
 600
 601         else
 602
 603           if (p_hcho  >= p1st) gas(nghcho  ) = rbox(p_hcho )*factgas
 604           if (p_ora1  >= p1st) gas(nghcooh ) = rbox(p_ora1 )*factgas
 605           if (p_ho    >= p1st) gas(ngoh    ) = rbox(p_ho   )*factgas
 606           if (p_ch3o2 >= p1st) gas(ngch3o2 ) = rbox(p_ch3o2)*factgas
 607           if (p_ch3oh >= p1st) gas(ngch3oh ) = rbox(p_ch3oh)*factgas
 608           if (p_op1   >= p1st) gas(ngch3o2h) = rbox(p_op1  )*factgas
 609
 610         endif
 611
 612
 613 !   compute bulk activated-aerosol mixing ratios
 614         aerosol(:) = 0.0
 615         rbulk_cwaer(:,:) = 0.0
 616
 617         iphase = cw_phase
 618         do itype = 1, ntype_aer
 619         do isize = 1, nsize_aer(itype)
 620         if ( .not. do_cloudchem_aer(isize,itype) ) cycle
 621
 622         do lyyy = 1, nyyy
 623
 624         l = lptr_yyy_cwaer(isize,itype,lyyy)
 625         if (l .ge. p1st) rbulk_cwaer(lyyy,1) = rbulk_cwaer(lyyy,1) + rbox(l)
 626
 627         end do
 628
 629         end do
 630         end do
 631
 632 !   map them to 'aerosol' array and convert to ug/m3
 633         aerosol(na4) = rbulk_cwaer(l_so4_aqyy,1) * factaerso4
 634         aerosol(nan) = rbulk_cwaer(l_no3_aqyy,1) * factaerno3
 635         aerosol(nac) = rbulk_cwaer(l_cl_aqyy, 1) * factaercl
 636         aerosol(naa) = rbulk_cwaer(l_nh4_aqyy,1) * factaernh4
 637         aerosol(nas) = rbulk_cwaer(l_na_aqyy, 1) * factaerna
 638         aerosol(nar) = rbulk_cwaer(l_oin_aqyy,1) * factaeroin
 639         aerosol(nae) = rbulk_cwaer(l_bc_aqyy, 1) * factaerbc
 640         aerosol(nao) = rbulk_cwaer(l_oc_aqyy, 1) * factaeroc
 641
 642
 643 !
 644 !   make call to aqoperator1
 645 !
 646 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
 647         lunxx = 87
 648         lunxx = -1
 649         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
 650         if (lunxx .gt. 0) then
 651             write(lunxx,*)
 652             write(lunxx,*)
 653             write(lunxx,*) 'interface_to_aqoperator1 - icase, irad, idecomp'
 654             write(lunxx,9870) icase, iradical_in, idecomp_hmsa_hso5
 655             write(lunxx,*) 'it, jt, kt, ktau'
 656             write(lunxx,9870) it, jt, kt, ktau
 657             write(lunxx,*) 'temp, p_atm, lwc, photo, co2, xprescribe_ph'
 658             write(lunxx,9875) temp, p_atm, lwc, photo_in, co2_mixrat_in, xprescribe_ph
 659             write(lunxx,*) 'pres_box, rho_box, qcw_box, dt_sec'
 660             write(lunxx,9875) pres_box, rho_box, qcw_box,   &
 661                 (tstep_end_sec-tstep_beg_sec)
 662             write(lunxx,*) 'gas (1=nh3, 2=hno3, 3=hcl, 4=h2so4, 11=so2, 12=h2o2, 18=o3)'
 663             write(lunxx,9875) gas
 664             write(lunxx,*) 'rbox(nh3, hno3, hcl, h2so4, so2, h2o2, o3)'
 665             write(lunxx,9875) rbox(p_nh3), rbox(p_hno3), rbox(p_hcl),   &
 666                 rbox(p_sulf), rbox(p_so2), rbox(p_h2o2), rbox(p_o3)
 667             write(lunxx,*) 'aerosol (1=na, 3=nh4, 4=no3, 5=cl, 6=so4, 8=ec, 9=oc, 10=crus)'
 668             write(lunxx,9875) aerosol
 669             write(lunxx,*) 'rbulk_cwaer (1=so4, 2=no3, 3-cl, 4=nh4, 5=na, 6=oin, 7=bc, 8=oc)'
 670             write(lunxx,9875) rbulk_cwaer(:,1)
 671             if (icase .le. -5) then
 672                 write(*,*)   &
 673                  '*** stopping in interface_to_aqop1 at icase =', icase
 674                 call wrf_error_fatal('*** stopping in interface_to_aqop1')
 675             end if
 676         end if
 677 9870    format( 8i5 )
 678 9875    format( 5(1pe14.6) )
 679 #endif
 680
 681 #if 0
 682 ! Print outs for debugging of aqoperator1... wig, 26-Oct-2005
 683 !!$    if( (id == 1 .and. ktau >=  207 ) .or. &
 684 !!$        (id == 2 .and. ktau >=  610 ) .or. &
 685 !!$        (id == 3 .and. ktau >= 1830 ) ) then
 686        write(6,'(a)') '---Begin input for aqoperator1---'
 687        write(6,'(a,4i)') 'id, it, jt, kt =', id, it, jt, kt
 688        write(6,'(a,1p,2e20.12)') 'tstep_beg_sec, tstep_end_sec = ',  &
 689            tstep_beg_sec, tstep_end_sec
 690        do l=1,ngas
 691           write(6,'("gas(",i2,") = ",1p,1e20.12)') l, gas(l)
 692        end do
 693        do l=1,naers
 694           write(6,'("aerosol(",i2,") = ",1p,1e20.12)') l, aerosol(l)
 695        end do
 696        write(6,'(a,1p,4e20.12)') "temp, p_atm, lwc, rh = ", temp, p_atm, lwc, rh
 697        write(6,'(a,1p,3e20.12)') "co2_mixrat_in, photo_in, xprescribe_ph = ",   &
 698            co2_mixrat_in, photo_in, xprescribe_ph
 699        write(6,'(a,3i)') " iradical_in, idecomp_hmsa_hso5, iaq = ",   &
 700            iradical_in, idecomp_hmsa_hso5, iaq
 701        write(6,'(a)') "---End input for aqoperator1---"
 702 !!$    end if
 703 #endif
 704
 705
 706 !   convert arguments to double prec
 707         tstep_beg_sec_dp = 0.0d0
 708         if (tstep_beg_sec .ne. 0.0) tstep_beg_sec_dp = tstep_beg_sec
 709         tstep_end_sec_dp = 0.0d0
 710         if (tstep_end_sec .ne. 0.0) tstep_end_sec_dp = tstep_end_sec
 711         temp_dp = 0.0d0
 712         if (temp .ne. 0.0) temp_dp = temp
 713         p_atm_dp = 0.0d0
 714         if (p_atm .ne. 0.0) p_atm_dp = p_atm
 715         lwc_dp = 0.0d0
 716         if (lwc .ne. 0.0) lwc_dp = lwc
 717         rh_dp = 0.0d0
 718         if (rh .ne. 0.0) rh_dp = rh
 719         co2_mixrat_in_dp = 0.0d0
 720         if (co2_mixrat_in .ne. 0.0) co2_mixrat_in_dp = co2_mixrat_in
 721         photo_in_dp = 0.0d0
 722         if (photo_in .ne. 0.0) photo_in_dp = photo_in
 723         xprescribe_ph_dp = 0.0d0
 724         if (xprescribe_ph .ne. 0.0) xprescribe_ph_dp = xprescribe_ph
 725         ph_cmuaq_cur_dp = 0.0d0
 726         if (ph_cmuaq_cur .ne. 0.0) ph_cmuaq_cur_dp = ph_cmuaq_cur
 727
 728         do i = 1, ngas
 729             gas_dp(i) = 0.0d0
 730             if (gas(i) .ne. 0.0) gas_dp(i) = gas(i)
 731         end do
 732         do i = 1, naers
 733             aerosol_dp(i) = 0.0d0
 734             if (aerosol(i) .ne. 0.0) aerosol_dp(i) = aerosol(i)
 735         end do
 736         do i = 1, ngas
 737             gas_aqfrac_cmu_dp(i) = 0.0d0
 738             if (gas_aqfrac_cmu(i) .ne. 0.0) gas_aqfrac_cmu_dp(i) = gas_aqfrac_cmu(i)
 739         end do
 740         do i = 1, meqn1max
 741             yaq_beg_dp(i) = 0.0d0
 742             if (yaq_beg(i) .ne. 0.0) yaq_beg_dp(i) = yaq_beg(i)
 743         end do
 744         do i = 1, meqn1max
 745             yaq_end_dp(i) = 0.0d0
 746             if (yaq_end(i) .ne. 0.0) yaq_end_dp(i) = yaq_end(i)
 747         end do
 748
 749
 750 !   total sulfur species conc as sulfate (ug/m3)
 751         cair_moleperm3 = 1.0e3*p_atm_dp/(0.082058e0*temp_dp)
 752         totsulf_beg = ( aerosol_dp(na4)/96.   &
 753                       + aerosol_dp(nahso5)/113. + aerosol_dp(nahmsa)/111.   &
 754                       + (gas_dp(ngso2) + gas_dp(ng4))*cair_moleperm3 )*96.0
 755
 756 !       call aqoperator1(   &
 757 !           istat_fatal, istat_warn,   &
 758 !           tstep_beg_sec, tstep_end_sec,   &
 759 !           gas, aerosol, gas_aqfrac_cmu,   &
 760 !           temp, p_atm, lwc, rh,   &
 761 !           co2_mixrat_in, photo_in, xprescribe_ph,   &
 762 !           iradical_in, idecomp_hmsa_hso5, iaq,   &
 763 !           yaq_beg, yaq_end, ph_cmuaq_cur )
 764
 765         call aqoperator1(   &
 766             istat_fatal, istat_warn,   &
 767             tstep_beg_sec_dp, tstep_end_sec_dp,   &
 768             gas_dp, aerosol_dp, gas_aqfrac_cmu_dp,   &
 769             temp_dp, p_atm_dp, lwc_dp, rh_dp,   &
 770             co2_mixrat_in_dp, photo_in_dp, xprescribe_ph_dp,   &
 771             iradical_in, idecomp_hmsa_hso5, iaq,   &
 772             yaq_beg_dp, yaq_end_dp, ph_cmuaq_cur_dp )
 773
 774         totsulf_end = ( aerosol_dp(na4)/96.   &
 775                       + aerosol_dp(nahso5)/113. + aerosol_dp(nahmsa)/111.   &
 776                       + (gas_dp(ngso2) + gas_dp(ng4))*cair_moleperm3 )*96.0
 777
 778
 779 !   convert arguments back to single prec
 780         tstep_beg_sec = tstep_beg_sec_dp
 781         tstep_end_sec = tstep_end_sec_dp
 782         temp = temp_dp
 783         p_atm = p_atm_dp
 784         lwc = lwc_dp
 785         rh = rh_dp
 786 !       co2_mixrat_in = co2_mixrat_in_dp        ! this has intent(in)
 787 !       photo_in = photo_in_dp                  ! this has intent(in)
 788 !       xprescribe_ph = xprescribe_ph_dp        ! this has intent(in)
 789         ph_cmuaq_cur = ph_cmuaq_cur_dp
 790
 791         do i = 1, ngas
 792             gas(i) = gas_dp(i)
 793         end do
 794         do i = 1, naers
 795             aerosol(i) = aerosol_dp(i)
 796         end do
 797         do i = 1, ngas
 798             gas_aqfrac_cmu(i) = gas_aqfrac_cmu_dp(i)
 799         end do
 800         do i = 1, meqn1max
 801             yaq_beg(i) = yaq_beg_dp(i)
 802         end do
 803         do i = 1, meqn1max
 804             yaq_end(i) = yaq_end_dp(i)
 805         end do
 806
 807
 808 !
 809 !   warning message when status flags are non-zero
 810 !
 811         istat_aqop = 0
 812         if (istat_fatal .ne. 0) then
 813             write(6,*)   &
 814                 '*** sorgam_cloudchem_driver, subr interface_to_aqoperator1'
 815             write(6,'(a,4i5,2i10)')   &
 816                 '    id,it,jt,kt, istat_fatal, warn =',   &
 817                 id, it, jt, kt, istat_fatal, istat_warn
 818             istat_aqop = -10
 819         end if
 820
 821 !
 822 !   warning message when sulfur mass balance error exceeds the greater
 823 !       of (1.0e-3 ug/m3) OR (1.0e-3 X total sulfur mixing ratio)
 824 !
 825         dum = totsulf_end - totsulf_beg
 826         dumb = max( totsulf_beg, totsulf_end )
 827         if (abs(dum) .gt. max(1.0e-3,1.0e-3*dumb)) then
 828             write(6,*)   &
 829                 '*** sorgam_cloudchem_driver, sulfur balance warning'
 830             write(6,'(a,4i5,1p,3e12.4)')   &
 831                 '    id,it,jt,kt, total_sulfur_beg, _end, _error =',   &
 832                 id, it, jt, kt, totsulf_beg, totsulf_end, dum
 833         end if
 834
 835 !
 836 !   map from [gas,aerosol,gas_aqfrac_box] to [rbox,gas_aqfrac_box]
 837 !
 838         gas_aqfrac_box(:) = 0.0
 839
 840         if (p_nh3   >= p1st) then
 841             rbox(p_nh3  ) = gas(nga     )/factgas
 842             if (p_nh3   <= numgas_aqfrac)   &
 843                 gas_aqfrac_box(p_nh3   ) = gas_aqfrac_cmu(nga     )
 844         end if
 845         if (p_hno3  >= p1st) then
 846             rbox(p_hno3 ) = gas(ngn     )/factgas
 847             if (p_hno3  <= numgas_aqfrac)   &
 848                 gas_aqfrac_box(p_hno3  ) = gas_aqfrac_cmu(ngn     )
 849         end if
 850         if (p_hcl   >= p1st) then
 851             rbox(p_hcl  ) = gas(ngc     )/factgas
 852             if (p_hcl   <= numgas_aqfrac)   &
 853                 gas_aqfrac_box(p_hcl   ) = gas_aqfrac_cmu(ngc     )
 854         end if
 855         if (p_sulf  >= p1st) then
 856             rbox(p_sulf ) = gas(ng4     )/factgas
 857             if (p_sulf  <= numgas_aqfrac)   &
 858                 gas_aqfrac_box(p_sulf  ) = gas_aqfrac_cmu(ng4     )
 859         end if
 860
 861 !       if (p_hcho  >= p1st) then
 862 !           rbox(p_hcho ) = gas(nghcho  )/factgas
 863 !           if (p_hcho  <= numgas_aqfrac)   &
 864 !               gas_aqfrac_box(p_hcho  ) = gas_aqfrac_cmu(nghcho  )
 865 !       end if
 866 !       if (p_ora1  >= p1st) then
 867 !           rbox(p_ora1 ) = gas(nghcooh )/factgas
 868 !           if (p_ora1  <= numgas_aqfrac)   &
 869 !               gas_aqfrac_box(p_ora1  ) = gas_aqfrac_cmu(nghcooh )
 870 !       end if
 871         if (p_so2   >= p1st) then
 872             rbox(p_so2  ) = gas(ngso2   )/factgas
 873             if (p_so2   <= numgas_aqfrac)   &
 874                 gas_aqfrac_box(p_so2   ) = gas_aqfrac_cmu(ngso2   )
 875         end if
 876         if (p_h2o2  >= p1st) then
 877             rbox(p_h2o2 ) = gas(ngh2o2  )/factgas
 878             if (p_h2o2  <= numgas_aqfrac)   &
 879                 gas_aqfrac_box(p_h2o2  ) = gas_aqfrac_cmu(ngh2o2  )
 880         end if
 881         if (p_o3    >= p1st) then
 882             rbox(p_o3   ) = gas(ngo3    )/factgas
 883             if (p_o3    <= numgas_aqfrac)   &
 884                 gas_aqfrac_box(p_o3    ) = gas_aqfrac_cmu(ngo3    )
 885         end if
 886 !       if (p_ho    >= p1st) then
 887 !           rbox(p_ho   ) = gas(ngoh    )/factgas
 888 !           if (p_ho    <= numgas_aqfrac)   &
 889 !               gas_aqfrac_box(p_ho    ) = gas_aqfrac_cmu(ngoh    )
 890 !       end if
 891         if (p_ho2   >= p1st) then
 892             rbox(p_ho2  ) = gas(ngho2   )/factgas
 893             if (p_ho2   <= numgas_aqfrac)   &
 894                 gas_aqfrac_box(p_ho2   ) = gas_aqfrac_cmu(ngho2   )
 895         end if
 896         if (p_no3   >= p1st) then
 897             rbox(p_no3  ) = gas(ngno3   )/factgas
 898             if (p_no3   <= numgas_aqfrac)   &
 899                 gas_aqfrac_box(p_no3   ) = gas_aqfrac_cmu(ngno3   )
 900         end if
 901
 902         if (p_no    >= p1st) then
 903             rbox(p_no   ) = gas(ngno    )/factgas
 904             if (p_no    <= numgas_aqfrac)   &
 905                 gas_aqfrac_box(p_no    ) = gas_aqfrac_cmu(ngno    )
 906         end if
 907         if (p_no2   >= p1st) then
 908             rbox(p_no2  ) = gas(ngno2   )/factgas
 909             if (p_no2   <= numgas_aqfrac)   &
 910                 gas_aqfrac_box(p_no2   ) = gas_aqfrac_cmu(ngno2   )
 911         end if
 912         if (p_hono  >= p1st) then
 913             rbox(p_hono ) = gas(nghno2  )/factgas
 914             if (p_hono  <= numgas_aqfrac)   &
 915                 gas_aqfrac_box(p_hono  ) = gas_aqfrac_cmu(nghno2  )
 916         end if
 917         if (p_pan   >= p1st) then
 918             rbox(p_pan  ) = gas(ngpan   )/factgas
 919             if (p_pan   <= numgas_aqfrac)   &
 920                 gas_aqfrac_box(p_pan   ) = gas_aqfrac_cmu(ngpan   )
 921         end if
 922 !       if (p_ch3o2 >= p1st) then
 923 !           rbox(p_ch3o2) = gas(ngch3o2 )/factgas
 924 !           if (p_ch3o2 <= numgas_aqfrac)   &
 925 !               gas_aqfrac_box(p_ch3o2 ) = gas_aqfrac_cmu(ngch3o2 )
 926 !       end if
 927 !       if (p_ch3oh >= p1st) then
 928 !           rbox(p_ch3oh) = gas(ngch3oh )/factgas
 929 !           if (p_ch3oh <= numgas_aqfrac)   &
 930 !               gas_aqfrac_box(p_ch3oh ) = gas_aqfrac_cmu(ngch3oh )
 931 !       end if
 932 !       if (p_op1   >= p1st) then
 933 !           rbox(p_op1  ) = gas(ngch3o2h)/factgas
 934 !           if (p_op1   <= numgas_aqfrac)   &
 935 !               gas_aqfrac_box(p_op1   ) = gas_aqfrac_cmu(ngch3o2h)
 936 !       end if
 937
 938 ! CB05 case added
 939        if ( (config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP) ) then
 940
 941           if (p_form  >= p1st) then
 942               rbox(p_form ) = gas(nghcho  )/factgas
 943               if (p_form  <= numgas_aqfrac)   &
 944                   gas_aqfrac_box(p_form  ) = gas_aqfrac_cmu(nghcho  )
 945           end if
 946           if (p_facd  >= p1st) then
 947               rbox(p_facd ) = gas(nghcooh )/factgas
 948               if (p_facd  <= numgas_aqfrac)   &
 949                   gas_aqfrac_box(p_facd  ) = gas_aqfrac_cmu(nghcooh )
 950           end if
 951           if (p_oh    >= p1st) then
 952               rbox(p_oh   ) = gas(ngoh    )/factgas
 953               if (p_oh    <= numgas_aqfrac)   &
 954                   gas_aqfrac_box(p_oh    ) = gas_aqfrac_cmu(ngoh    )
 955           end if
 956           if (p_meo2  >= p1st) then
 957               rbox(p_meo2 ) = gas(ngch3o2 )/factgas
 958               if (p_meo2  <= numgas_aqfrac)   &
 959                   gas_aqfrac_box(p_meo2  ) = gas_aqfrac_cmu(ngch3o2 )
 960           end if
 961
 962           if (p_meoh  >= p1st) then
 963               rbox(p_meoh ) = gas(ngch3oh )/factgas
 964               if (p_meoh  <= numgas_aqfrac)   &
 965                   gas_aqfrac_box(p_meoh  ) = gas_aqfrac_cmu(ngch3oh )
 966           end if
 967           if (p_mepx  >= p1st) then
 968               rbox(p_mepx ) = gas(ngch3o2h)/factgas
 969               if (p_mepx  <= numgas_aqfrac)   &
 970                   gas_aqfrac_box(p_mepx  ) = gas_aqfrac_cmu(ngch3o2h)
 971           end if
 972
 973        else
 974           if (p_hcho  >= p1st) then
 975               rbox(p_hcho ) = gas(nghcho  )/factgas
 976               if (p_hcho  <= numgas_aqfrac)   &
 977                   gas_aqfrac_box(p_hcho  ) = gas_aqfrac_cmu(nghcho  )
 978           end if
 979           if (p_ora1  >= p1st) then
 980               rbox(p_ora1 ) = gas(nghcooh )/factgas
 981               if (p_ora1  <= numgas_aqfrac)   &
 982                   gas_aqfrac_box(p_ora1  ) = gas_aqfrac_cmu(nghcooh )
 983           end if
 984           if (p_ho    >= p1st) then
 985               rbox(p_ho   ) = gas(ngoh    )/factgas
 986               if (p_ho    <= numgas_aqfrac)   &
 987                   gas_aqfrac_box(p_ho    ) = gas_aqfrac_cmu(ngoh    )
 988           end if
 989           if (p_ch3o2 >= p1st) then
 990               rbox(p_ch3o2) = gas(ngch3o2 )/factgas
 991               if (p_ch3o2 <= numgas_aqfrac)   &
 992                   gas_aqfrac_box(p_ch3o2 ) = gas_aqfrac_cmu(ngch3o2 )
 993           end if
 994           if (p_ch3oh >= p1st) then
 995               rbox(p_ch3oh) = gas(ngch3oh )/factgas
 996               if (p_ch3oh <= numgas_aqfrac)   &
 997                   gas_aqfrac_box(p_ch3oh ) = gas_aqfrac_cmu(ngch3oh )
 998           end if
 999           if (p_op1   >= p1st) then
1000               rbox(p_op1  ) = gas(ngch3o2h)/factgas
1001               if (p_op1   <= numgas_aqfrac)   &
1002                   gas_aqfrac_box(p_op1   ) = gas_aqfrac_cmu(ngch3o2h)
1003           end if
1004
1005        end if
1006
1007 ! end of addition
1008
1009         rbulk_cwaer(l_so4_aqyy,2) = aerosol(na4)/factaerso4
1010         rbulk_cwaer(l_no3_aqyy,2) = aerosol(nan)/factaerno3
1011         rbulk_cwaer(l_cl_aqyy, 2) = aerosol(nac)/factaercl
1012         rbulk_cwaer(l_nh4_aqyy,2) = aerosol(naa)/factaernh4
1013         rbulk_cwaer(l_na_aqyy, 2) = aerosol(nas)/factaerna
1014         rbulk_cwaer(l_oin_aqyy,2) = aerosol(nar)/factaeroin
1015         rbulk_cwaer(l_bc_aqyy, 2) = aerosol(nae)/factaerbc
1016         rbulk_cwaer(l_oc_aqyy, 2) = aerosol(nao)/factaeroc
1017
1018
1019 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1020         lunxx = 87
1021         lunxx = -1
1022         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1023         if (lunxx .gt. 0) then
1024             write(lunxx,*)
1025             write(lunxx,*) 'interface_to_aqoperator1 - after call'
1026             write(lunxx,*) 'gas (1=nh3, 2=hno3, 3=hcl, 4=h2so4, 11=so2, 12=h2o2, 18=o3)'
1027             write(lunxx,9875) gas
1028             write(lunxx,*) 'rbox(nh3, hno3, hcl, h2so4, so2, h2o2, o3)'
1029             write(lunxx,9875) rbox(p_nh3), rbox(p_hno3), rbox(p_hcl),   &
1030                 rbox(p_sulf), rbox(p_so2), rbox(p_h2o2), rbox(p_o3)
1031             write(lunxx,*) 'aerosol (1=na, 3=nh4, 4=no3, 5=cl, 6=so4, 8=ec, 9=oc, 10=crus)'
1032             write(lunxx,9875) aerosol
1033             write(lunxx,*) 'rbulk_cwaer (1=so4, 2=no3, 3-cl, 4=nh4, 5=na, 6=oin, 7=bc, 8=oc)'
1034             write(lunxx,9875) rbulk_cwaer(:,2)
1035             write(lunxx,*) 'ph_cmuaq_cur'
1036             write(lunxx,9875) ph_cmuaq_cur
1037             if (icase .le. -5) then
1038                 write(*,*)   &
1039                  '*** stopping in interface_to_aqop1 at icase =', icase
1040                 call wrf_error_fatal('*** stopping in interface_to_aqop1')
1041             end if
1042         end if
1043 #endif
1044
1045
1046         return
1047         end subroutine sorgam_interface_to_aqoperator1
1048
1049
1050
1051 !-----------------------------------------------------------------------
1052         subroutine sorgam_partition_cldwtr(   &
1053             rbox, fr_partit_cw, xvol_old,   &
1054             id, it, jt, kt, icase )
1055
1056         use module_state_description, only:   &
1057                 param_first_scalar, num_chem
1058
1059         use module_data_sorgam_vbs, only:   &
1060                 maxd_asize, maxd_atype,   &
1061                 ai_phase, cw_phase, nsize_aer, ntype_aer, ncomp_aer,   &
1062                 do_cloudchem_aer, massptr_aer, numptr_aer,   &
1063                 dens_aer, sigmag_aer,   &
1064                 dcen_sect, dlo_sect, dhi_sect,   &
1065                 volumcen_sect, volumlo_sect, volumhi_sect
1066
1067
1068         implicit none
1069
1070 !   subr arguments
1071         integer, intent(in) :: id, it, jt, kt, icase
1072
1073         real, intent(inout), dimension( 1:num_chem ) :: rbox
1074
1075         real, intent(inout), dimension( maxd_asize, maxd_atype ) ::   &
1076                 fr_partit_cw
1077
1078         real, intent(inout), dimension( 2, 3 ) :: xvol_old
1079
1080 !   local variables
1081         integer :: isize, itype
1082         integer :: jdone_mass, jdone_numb, jpos, jpos_mass, jpos_numb
1083         integer :: l, ll
1084         integer :: p1st
1085 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1086         integer :: lunxx
1087 #endif
1088
1089         real, parameter :: partit_wght_mass = 0.5
1090
1091         real :: tmpa, tmpb, tmpc
1092         real :: tmp_cwvolfrac, tmp_lnsg
1093         real :: tmass, tnumb, umass, unumb, wmass, wnumb
1094         real :: xmass_c, xmass_a, xmass_t, xvolu_c, xvolu_a, xvolu_t
1095         real :: xnumb_c1, xnumb_a1, xnumb_t1, xnumb_c2, xnumb_a2, xnumb_t2
1096         real, dimension( maxd_asize, maxd_atype ) :: fmass, fnumb, xmass, xnumb, xnumbsv
1097
1098
1099         p1st = PARAM_FIRST_SCALAR
1100
1101         tmass = 0.0
1102         tnumb = 0.0
1103         umass = 0.0
1104         unumb = 0.0
1105
1106 !   compute
1107 !     xmass, xnumb = mass, number mixing ratio for a bin
1108 !     tmass, tnumb = sum over all bins of xmass, xnumb
1109 !     umass, unumb = max over all bins of xmass, xnumb
1110 !   set xmass, xnumb = 0.0 if bin mass, numb < 1.0e-37
1111 !   constrain xnumb so that mean particle volume is
1112 !     within bin boundaries
1113 !   for made-sorgam, x/t/umass are g/kg-air, x/t/unumb are #/kg-air
1114         do itype = 1, ntype_aer
1115         do isize = 1, nsize_aer(itype)
1116             if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1117             xmass_c = 0.0
1118             xvolu_c = 0.0
1119             xvolu_a = 0.0
1120             do ll = 1, ncomp_aer(itype)
1121                 l = massptr_aer(ll,isize,itype,cw_phase)
1122                 if (l .ge. p1st) then
1123                     tmpa = max( 0.0, rbox(l) )*1.0e-6
1124                     xmass_c = xmass_c + tmpa
1125                     xvolu_c = xvolu_c + tmpa/dens_aer(ll,itype)
1126                 end if
1127                 l = massptr_aer(ll,isize,itype,ai_phase)
1128                 if (l .ge. p1st) then
1129                     tmpa = max( 0.0, rbox(l) )*1.0e-6
1130                     xvolu_a = xvolu_a + tmpa/dens_aer(ll,itype)
1131                 end if
1132             end do
1133
1134             xnumb_c1 = max( 0.0, rbox(numptr_aer(isize,itype,cw_phase)) )
1135             xnumb_a1 = max( 0.0, rbox(numptr_aer(isize,itype,ai_phase)) )
1136             xnumbsv(isize,itype) = xnumb_c1
1137             xnumb_t1 = xnumb_a1 + xnumb_c1
1138             xvolu_t = xvolu_a + xvolu_c
1139
1140 !   do "bounding" activated+interstitial combined number
1141 !   and calculate dgnum for activated+interstitial combined
1142             if (xvolu_t < smallvolaa) then
1143                 xnumb_t2 = xvolu_t/volumcen_sect(isize,itype)
1144             else if (xnumb_t1 < xvolu_t/volumhi_sect(isize,itype)) then
1145                 xnumb_t2 = xvolu_t/volumhi_sect(isize,itype)
1146             else if (xnumb_t1 > xvolu_t/volumlo_sect(isize,itype)) then
1147                 xnumb_t2 = xvolu_t/volumlo_sect(isize,itype)
1148             else
1149                 xnumb_t2 = xnumb_t1
1150             end if
1151
1152 !   do "bounding" of activated number
1153 !   tmp_cwvolfrac = (cw volume)/(ai + cw volume)
1154             tmp_cwvolfrac = xvolu_c/max(xvolu_t,1.e-30)
1155             tmp_lnsg = log(sigmag_aer(isize,itype))
1156             if ((xvolu_c < smallvolaa) .or. (tmp_cwvolfrac < 1.0e-10)) then
1157 !   for very small cw volume or volume fraction,
1158 !   use (ai+cw number)*(cw volume fraction)
1159                 xnumb_c2 = xnumb_t2*tmp_cwvolfrac
1160                 tmpa = -7.0 ; tmpb = -7.0 ; tmpc = -7.0
1161             else
1162 !   tmpa is value of (ln(dpcut)-ln(dgvol))/ln(sigmag) for which
1163 !   "norm01 upper tail" cummulative pdf is equal to tmp_cwvolfrac
1164                 tmpa = norm01_uptail_inv( tmp_cwvolfrac )
1165 !   tmpb is corresponding value of (ln(dp)-ln(dgnum))/ln(sigmag)
1166                 tmpb = tmpa + 3.0*tmp_lnsg
1167 !   tmpc is corresponding "norm01 upper tail" cummulative pdf
1168                 tmpc = norm01_uptail( tmpb )
1169 !   minimum number of activated particles occurs when
1170 !   activated particles are dp>dpcut and interstitial are dp<dpcut,
1171 !   giving (cw number) == (ai+cw number)*tmpc
1172                 xnumb_c2 = max( xnumb_c1,  xnumb_t2*tmpc )
1173 !   assuming activated particles are at least as big as interstitial ones,
1174 !   the minimum number of activated particles occurs when ai & cw have same sizes,
1175 !   giving (cw number) == (ai+cw number)*(cw volume fraction)
1176                 xnumb_c2 = min( xnumb_c2,  xnumb_t2*tmp_cwvolfrac )
1177             end if
1178             xnumb_a2 = xnumb_t2 - xnumb_c2
1179             rbox(numptr_aer(isize,itype,cw_phase)) = xnumb_c2
1180             rbox(numptr_aer(isize,itype,ai_phase)) = xnumb_a2
1181
1182 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1183 !   diagnostics when lunxx > 0
1184         lunxx = 86
1185         lunxx = -1
1186         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1187         if (lunxx > 0) then
1188             if ((isize == 1) .and. (itype == 1)) write(lunxx,'(a)')
1189             write(lunxx,'(a,3i4,4x,2i4,2x,l2)') 'partition-cw  i,j,k,  is,it', &
1190                 it, jt, kt, isize, itype
1191             write(lunxx,'(a,1p,5e12.4)') 'cw  vol, num old/adj      ', &
1192                 xvolu_c, xnumb_c1, xnumb_c2
1193             write(lunxx,'(a,1p,5e12.4)') 'ai  vol, num old/adj      ', &
1194                 xvolu_a, xnumb_a1, xnumb_a2
1195             write(lunxx,'(a,1p,5e12.4)') 'a+c vol, num old/adj      ', &
1196                 xvolu_t, xnumb_t1, xnumb_t2
1197             write(lunxx,'(a,1p,5e12.4)') 'lnsg, cwvolfr, cwnumfr 1/2', &
1198                 tmp_lnsg, tmp_cwvolfrac, &
1199                 xnumb_c1/max(xnumb_t1,1.0e-30), &
1200                 xnumb_c2/max(xnumb_t2,1.0e-30)
1201             write(lunxx,'(a,1p,5e12.4)') 'tmpa/b/c                  ', &
1202                 tmpa, tmpb, tmpc
1203             write(lunxx,'(a,1p,5e12.4)') 'dlo, dcen, dhi_sect       ', &
1204                 dlo_sect(isize,itype), dcen_sect(isize,itype), &
1205                 dhi_sect(isize,itype)
1206             write(lunxx,'(a,1p,5e12.4)') 'vlo, vcen, vhi_sect       ', &
1207                 volumlo_sect(isize,itype), volumcen_sect(isize,itype), &
1208                 volumhi_sect(isize,itype)
1209         end if
1210 #endif
1211
1212             if (xmass_c .lt. 1.0e-37) xmass_c = 0.0
1213             xmass(isize,itype) = xmass_c
1214             if (xnumb_c2 .lt. 1.0e-37) xnumb_c2 = 0.0
1215             xnumb(isize,itype) = xnumb_c2
1216             xnumbsv(isize,itype) = xnumb_c1
1217
1218             tmass = tmass + xmass(isize,itype)
1219             tnumb = tnumb + xnumb(isize,itype)
1220             umass = max( umass, xmass(isize,itype) )
1221             unumb = max( unumb, xnumb(isize,itype) )
1222
1223             if ((itype == 1) .and. (isize <= 2)) then
1224                 xvol_old(isize,1) = xvolu_c
1225                 xvol_old(isize,2) = xvolu_a
1226                 xvol_old(isize,3) = xvolu_t
1227             end if
1228         end do
1229         end do
1230
1231 !   compute
1232 !     fmass, fnumb = fraction of total mass, number that is in a bin
1233 !   if tmass<1e-35 and umass>0, set fmass=1 for bin with largest xmass
1234 !   if tmass<1e-35 and umass=0, set fmass=0 for all
1235         jdone_mass = 0
1236         jdone_numb = 0
1237         jpos_mass = 0
1238         jpos_numb = 0
1239         do itype = 1, ntype_aer
1240         do isize = 1, nsize_aer(itype)
1241             if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1242             fmass(isize,itype) = 0.0
1243             if (tmass .ge. 1.0e-35) then
1244                 fmass(isize,itype) = xmass(isize,itype)/tmass
1245             else if (umass .gt. 0.0) then
1246                 if ( (jdone_mass .eq. 0) .and.   &
1247                      (xmass(isize,itype) .eq. umass) ) then
1248                     jdone_mass = 1
1249                     fmass(isize,itype) = 1.0
1250                 end if
1251             end if
1252             if (fmass(isize,itype) .gt. 0) jpos_mass = jpos_mass + 1
1253
1254             fnumb(isize,itype) = 0.0
1255             if (tnumb .ge. 1.0e-35) then
1256                 fnumb(isize,itype) = xnumb(isize,itype)/tnumb
1257             else if (unumb .gt. 0.0) then
1258                 if ( (jdone_numb .eq. 0) .and.   &
1259                      (xnumb(isize,itype) .eq. unumb) ) then
1260                     jdone_numb = 1
1261                     fnumb(isize,itype) = 1.0
1262                 end if
1263             end if
1264             if (fnumb(isize,itype) .gt. 0) jpos_numb = jpos_numb + 1
1265         end do
1266         end do
1267
1268 !   if only 1 bin has fmass or fnumb > 0, set value to 1.0 exactly
1269         if ((jpos_mass .eq. 1) .or. (jpos_numb .eq. 1)) then
1270         do itype = 1, ntype_aer
1271         do isize = 1, nsize_aer(itype)
1272             if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1273             if (jpos_mass .eq. 1) then
1274                 if (fmass(isize,itype) .gt. 0) fmass(isize,itype) = 1.0
1275             end if
1276             if (jpos_numb .eq. 1) then
1277                 if (fnumb(isize,itype) .gt. 0) fnumb(isize,itype) = 1.0
1278             end if
1279         end do
1280         end do
1281         end if
1282
1283 !
1284 !   compute fr_partit_cw as weighted average of fmass & fnumb, except
1285 !   if tmass<1e-35 and umass=0, use only fnumb
1286 !   if tnumb<1e-35 and unumb=0, use only fmass
1287 !   if tmass,tnumb<1e-35 and umass,unumb=0,
1288 !               set fr_partit_cw=1 for center bin of itype=1
1289 !
1290         fr_partit_cw(:,:) = 0.0
1291         if ((jpos_mass .eq. 0) .and. (jpos_numb .eq. 0)) then
1292             itype = 1
1293             isize = (nsize_aer(itype)+1)/2
1294             fr_partit_cw(isize,itype) = 1.0
1295
1296         else if (jpos_mass .eq. 0) then
1297             fr_partit_cw(:,:) = fnumb(:,:)
1298
1299         else if (jpos_numb .eq. 0) then
1300             fr_partit_cw(:,:) = fmass(:,:)
1301
1302         else
1303             wmass = max( 0.0, min( 1.0, partit_wght_mass ) )
1304             wnumb = 1.0 - wmass
1305             fr_partit_cw(:,:) =  wmass*fmass(:,:) + wnumb*fnumb(:,:)
1306
1307             jpos = 0
1308             do itype = 1, ntype_aer
1309             do isize = 1, nsize_aer(itype)
1310                 if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1311                 if (fr_partit_cw(isize,itype) .gt. 0.0) jpos = jpos + 1
1312             end do
1313             end do
1314
1315 !   if only 1 bin has fr_partit_cw > 0, set value to 1.0 exactly
1316             if (jpos .eq. 1) then
1317             do itype = 1, ntype_aer
1318             do isize = 1, nsize_aer(itype)
1319                 if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1320                 if (fr_partit_cw(isize,itype) .gt. 0.0)   &
1321                         fr_partit_cw(isize,itype) = 1.0
1322             end do
1323             end do
1324             end if
1325         end if
1326
1327
1328 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1329 !   diagnostics when lunxx > 0
1330         lunxx = 86
1331         lunxx = -1
1332         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1333 !       if (icase .gt. 9) lunxx = -1
1334         if (lunxx .gt. 0) then
1335             write(lunxx,9800)
1336             write(lunxx,9800)   &
1337                 'partition_cldwtr - icase, jpos, jpos_mass, jpos_numb'
1338             write(lunxx,9810)  icase, jpos, jpos_mass, jpos_numb
1339             write(lunxx,9800) 'tmass, umass, wmass'
1340             write(lunxx,9820)  tmass, umass, wmass
1341             write(lunxx,9800) 'tnumb, unumb, wnumb'
1342             write(lunxx,9820)  tnumb, unumb, wnumb
1343             write(lunxx,9800) 'xmass, fmass, xnumb_orig/adj, fnumb, fr_partit_cw'
1344             tmpa = 0.0
1345             tmpb = 0.0
1346             tmpc = 0.0
1347             do itype = 1, ntype_aer
1348             do isize = 1, nsize_aer(itype)
1349                 if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1350                 write(lunxx,9820)  xmass(isize,itype), fmass(isize,itype),   &
1351                         xnumbsv(isize,itype), xnumb(isize,itype),   &
1352                         fnumb(isize,itype), fr_partit_cw(isize,itype)
1353                 tmpa = tmpa + fmass(isize,itype)
1354                 tmpb = tmpb + fnumb(isize,itype)
1355                 tmpc = tmpc + fr_partit_cw(isize,itype)
1356             end do
1357             end do
1358             write(lunxx,9800)   &
1359                 'sum_fmass-1.0,  sum_fnumb-1.0,  sum_fr_partit-1.0'
1360             write(lunxx,9820)  (tmpa-1.0), (tmpb-1.0), (tmpc-1.0)
1361             if (icase .le. -5) then
1362                 write(*,*) '*** stopping in partition_cldwtr at icase =', icase
1363                 call wrf_error_fatal('*** stopping in partition_cldwtr')
1364             end if
1365 9800        format( a )
1366 9810        format( 5i10 )
1367 9820        format( 6(1pe10.2) )
1368         end if
1369 #endif
1370
1371
1372         return
1373         end subroutine sorgam_partition_cldwtr
1374
1375
1376
1377 !-----------------------------------------------------------------------
1378         subroutine sorgam_distribute_bulk_changes(   &
1379                 rbox, rbox_sv1, fr_partit_cw,   &
1380                 rbulk_cwaer, lptr_yyy_cwaer,   &
1381                 id, it, jt, kt, icase )
1382
1383         use module_state_description, only:   &
1384                 param_first_scalar, num_chem
1385
1386         use module_scalar_tables, only:  chem_dname_table
1387
1388         use module_data_sorgam_vbs, only:   &
1389                 maxd_asize, maxd_atype,   &
1390                 cw_phase, nsize_aer, ntype_aer, do_cloudchem_aer,   &
1391                 lptr_so4_aer, lptr_no3_aer, lptr_nh4_aer,   &
1392                 lptr_orgpa_aer, lptr_ec_aer, lptr_p25_aer
1393
1394
1395         implicit none
1396
1397 !   subr arguments
1398         integer, intent(in) :: id, it, jt, kt, icase
1399
1400         integer, intent(in) :: lptr_yyy_cwaer(maxd_asize,maxd_atype,nyyy)
1401
1402         real, intent(inout), dimension( 1:num_chem ) :: rbox, rbox_sv1
1403
1404         real, intent(in), dimension( maxd_asize, maxd_atype ) ::   &
1405                 fr_partit_cw
1406
1407         real, intent(in), dimension( nyyy, 2 ) :: rbulk_cwaer
1408
1409
1410 !   local variables
1411         integer :: iphase, isize, itype
1412         integer :: idone, icount, ncount
1413         integer :: jpos, jpos_sv
1414         integer :: l, lunxxaa, lunxxbb, lyyy
1415         integer :: p1st
1416 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1417         integer :: lunxx
1418 #endif
1419
1420         real :: duma, dumb, dumc
1421         real :: fr, frsum_cur
1422         real :: fr_cur(maxd_asize,maxd_atype)
1423         real :: del_r_current, del_r_remain
1424         real :: del_rbulk_cwaer(nyyy)
1425
1426
1427         p1st = param_first_scalar
1428
1429         do lyyy = 1, nyyy
1430             del_rbulk_cwaer(lyyy) = rbulk_cwaer(lyyy,2) - rbulk_cwaer(lyyy,1)
1431         end do
1432
1433         iphase = cw_phase
1434
1435
1436         jpos = 0
1437         do itype = 1, ntype_aer
1438         do isize = 1, nsize_aer(itype)
1439             if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1440             if (fr_partit_cw(isize,itype) .gt. 0) jpos = jpos + 1
1441         end do
1442         end do
1443         jpos_sv = jpos
1444
1445 !
1446 !   distribution is trivial when only 1 bin has fr_partit_cw > 0
1447 !
1448         if (jpos_sv .eq. 1) then
1449             do lyyy = 1, nyyy
1450
1451             do itype = 1, ntype_aer
1452             do isize = 1, nsize_aer(itype)
1453                 if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1454                 fr = fr_partit_cw(isize,itype)
1455                 if (fr .eq. 1.0) then
1456                     l = lptr_yyy_cwaer(isize,itype,lyyy)
1457                     if (l .ge. p1st) rbox(l) = rbulk_cwaer(lyyy,2)
1458                 end if
1459             end do
1460             end do
1461
1462             end do
1463             goto 7900
1464         end if
1465
1466
1467         do 3900 lyyy = 1, nyyy
1468
1469 !
1470 !   distribution is simple when del_rbulk_cwaer(lyyy) >= 0
1471 !
1472         if (del_rbulk_cwaer(lyyy) .eq. 0.0) then
1473             goto 3900
1474         else if (del_rbulk_cwaer(lyyy) .gt. 0.0) then
1475             do itype = 1, ntype_aer
1476             do isize = 1, nsize_aer(itype)
1477                 if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1478                 fr = fr_partit_cw(isize,itype)
1479                 if (fr .gt. 0.0) then
1480                     l = lptr_yyy_cwaer(isize,itype,lyyy)
1481                     if (l .ge. p1st) then
1482                         rbox(l) = rbox(l) + fr*del_rbulk_cwaer(lyyy)
1483                     end if
1484                 end if
1485             end do
1486             end do
1487
1488             goto 3900
1489         end if
1490
1491 !
1492 !   distribution is complicated when del_rbulk_cwaer(lyyy) < 0,
1493 !   because you cannot produce any negative mixrats
1494 !
1495         del_r_remain = del_rbulk_cwaer(lyyy)
1496         fr_cur(:,:) = fr_partit_cw(:,:)
1497
1498         ncount = max( 1, jpos_sv*2 )
1499         icount = 0
1500
1501 !   iteration loop
1502         do while (icount .le. ncount)
1503
1504         icount = icount + 1
1505         del_r_current = del_r_remain
1506         jpos = 0
1507         frsum_cur = 0.0
1508
1509         do itype = 1, ntype_aer
1510         do isize = 1, nsize_aer(itype)
1511         if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1512
1513         fr = fr_cur(isize,itype)
1514
1515         if (fr .gt. 0.0) then
1516             l = lptr_yyy_cwaer(isize,itype,lyyy)
1517             if (l .ge. p1st) then
1518                 duma = fr*del_r_current
1519                 dumb = rbox(l) + duma
1520                 if (dumb .gt. 0.0) then
1521                     jpos = jpos + 1
1522                 else if (dumb .eq. 0.0) then
1523                     fr_cur(isize,itype) = 0.0
1524                 else
1525                     duma = -rbox(l)
1526                     dumb = 0.0
1527                     fr_cur(isize,itype) = 0.0
1528                 end if
1529                 del_r_remain = del_r_remain - duma
1530                 rbox(l) = dumb
1531                 frsum_cur = frsum_cur + fr_cur(isize,itype)
1532             else
1533                 fr_cur(isize,itype) = 0.0
1534             end if
1535         end if
1536
1537         end do  ! isize = 1, nsize_aer
1538         end do  ! itype = 1, ntype_aer
1539
1540 !   done if jpos = jpos_sv, because bins reached zero mixrat
1541         if (jpos .eq. jpos_sv) then
1542             idone = 1
1543 !   del_r_remain starts as negative, so done if non-negative
1544         else if (del_r_remain .ge. 0.0) then
1545             idone = 2
1546 !   del_r_remain starts as negative, so done if non-negative
1547         else if (abs(del_r_remain) .le. 1.0e-7*abs(del_rbulk_cwaer(lyyy))) then
1548             idone = 3
1549 !   done if all bins have fr_cur = 0
1550         else if (frsum_cur .le. 0.0) then
1551             idone = 4
1552 !   same thing basically
1553         else if (jpos .le. 0) then
1554             idone = 5
1555         else
1556             idone = 0
1557         end if
1558
1559 !   check for done, and (conditionally) print message
1560         if (idone .gt. 0) then
1561             lunxxaa = 6
1562 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1563 !           lunxxaa = 86
1564             if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxxaa = 171
1565 #endif
1566             if ((lunxxaa .gt. 0) .and. (icount .gt. (1+jpos_sv)/2)) then
1567                 write(lunxxaa,9800)   &
1568                 'distribute_bulk_changes - icount>jpos_sv/2 - i,j,k'
1569                 write(lunxxaa,9810)  it, jt, kt
1570                 write(lunxxaa,9800) 'icase, lyyy, idone, icount, jpos, jpos_sv'
1571                 write(lunxxaa,9810)  icase, lyyy, idone, icount, jpos, jpos_sv
1572             end if
1573             goto 3900
1574         end if
1575
1576 !   rescale fr_cur for next iteration
1577         fr_cur(:,:) = fr_cur(:,:)/frsum_cur
1578
1579         end do  ! while (icount .le. ncount)
1580
1581
1582 !   icount > ncount, so print message
1583         lunxxbb = 6
1584 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1585 !       lunxxbb = 86
1586         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxxbb = 171
1587 #endif
1588         if (lunxxbb .gt. 0) then
1589             write(lunxxbb,9800)
1590             write(lunxxbb,9800)   &
1591                 'distribute_bulk_changes - icount>ncount - i,j,k'
1592             write(lunxxbb,9810)  it, jt, kt
1593             write(lunxxbb,9800) 'icase, lyyy, icount, ncount, jpos_sv, jpos'
1594             write(lunxxbb,9810)  icase, lyyy, icount, ncount, jpos_sv, jpos
1595             write(lunxxbb,9800) 'rbulk_cwaer(1), del_rbulk_cwaer, del_r_remain, frsum_cur, (frsum_cur-1.0)'
1596             write(lunxxbb,9820)  rbulk_cwaer(lyyy,1), del_rbulk_cwaer(lyyy),   &
1597                 del_r_remain, frsum_cur, (frsum_cur-1.0)
1598         end if
1599 9800    format( a )
1600 9801    format( 3a )
1601 9810    format( 7i10 )
1602 9820    format( 7(1pe10.2) )
1603 9840    format( 2i3, 5(1pe14.6) )
1604
1605
1606 3900    continue
1607
1608 7900    continue
1609
1610
1611 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1612 !   diagnostics for testing
1613         lunxx = 88
1614         lunxx = -1
1615         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1616         if (lunxx .gt. 0) then
1617             icount = 0
1618             do lyyy = 1, nyyy
1619             duma = del_rbulk_cwaer(lyyy)
1620             if ( abs(duma) .gt.   &
1621                     max( 1.0e-35, 1.0e-5*abs(rbulk_cwaer(lyyy,1)) ) ) then
1622                 icount = icount + 1
1623                 if (icount .eq. 1) write(lunxx,9800)
1624                 if (icount .eq. 1) write(lunxx,9800)
1625                 write(lunxx,9800)
1626                 l = lptr_yyy_cwaer(1,1,lyyy)
1627                 if (l .ge. p1st) then
1628                     write(lunxx,9801) 'distribute_bulk_changes - ',   &
1629                         chem_dname_table(id,l)(1:12), ' - icase, lyyy, l11'
1630                 else
1631                     write(lunxx,9801) 'distribute_bulk_changes - ',   &
1632                         'name = ?????', ' - icase, lyyy, l11'
1633                 end if
1634                 write(lunxx,9810) icase, lyyy, l
1635                 write(lunxx,9800) ' tp sz  rbox_sv1, rbox, del_rbox' //   &
1636                         ', del_rbox/del_rbulk_cwaer, (...-fr_partit_cw)'
1637                 write(lunxx,9840) 0, 0, rbulk_cwaer(lyyy,1:2), duma
1638                 do itype = 1, ntype_aer
1639                 do isize = 1, nsize_aer(itype)
1640                     if ( .not. do_cloudchem_aer(isize,itype) ) cycle
1641                     l = lptr_yyy_cwaer(isize,itype,lyyy)
1642                     if (l .lt. p1st) cycle
1643                     dumb = rbox(l) - rbox_sv1(l)
1644                     dumc = dumb/max( abs(duma), 1.0e-35 )
1645                     if (duma .lt. 0.0) dumc = -dumc
1646                     write(lunxx,9840) itype, isize, rbox_sv1(l), rbox(l),   &
1647                         dumb, dumc, (dumc-fr_partit_cw(isize,itype))
1648                 end do
1649                 end do
1650             end if
1651             end do
1652             if (icase .le. -5) then
1653                 write(*,*)   &
1654                 '*** stop in distribute_bulk_changes diags, icase =', icase
1655                 call wrf_error_fatal('*** stop in distribute_bulk_changes diags')
1656             end if
1657         end if
1658 #endif
1659
1660
1661         return
1662         end subroutine sorgam_distribute_bulk_changes
1663
1664
1665
1666 !-----------------------------------------------------------------------
1667         subroutine sorgam_cloudchem_apply_mode_transfer(   &
1668                 rbox, rbox_sv1, xvol_old,   &
1669                 id, it, jt, kt, icase )
1670
1671         use module_state_description, only:   &
1672                 param_first_scalar, num_chem
1673
1674         use module_scalar_tables, only:  chem_dname_table
1675
1676         use module_data_sorgam_vbs, only:   &
1677                 pirs,   &
1678                 msectional,   &
1679                 maxd_asize, maxd_atype,   &
1680                 ai_phase, cw_phase, nsize_aer, ntype_aer, ncomp_aer,   &
1681                 do_cloudchem_aer, massptr_aer, numptr_aer, dens_aer,   &
1682                 sigmag_aer, dcen_sect, dlo_sect, dhi_sect,   &
1683                 volumcen_sect, volumlo_sect, volumhi_sect,   &
1684                 lptr_so4_aer, lptr_nh4_aer, lptr_p25_aer
1685
1686         use module_aerosols_sorgam_vbs, only:  getaf
1687
1688
1689         implicit none
1690
1691 !   subr arguments
1692         integer, intent(in) :: id, it, jt, kt, icase
1693
1694         real, intent(inout), dimension( 1:num_chem ) :: rbox, rbox_sv1
1695
1696         real, intent(in), dimension( 2, 3 ) :: xvol_old
1697
1698
1699 !   local variables
1700         integer :: idum_msect
1701         integer :: ii, isize, isize_ait, isize_acc, itype
1702         integer :: jj
1703         integer :: l, lfrm, ltoo, ll
1704         integer :: lptr_dum(maxd_asize,maxd_atype)
1705         integer :: p1st
1706 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1707         integer :: lunxx
1708 #endif
1709
1710         logical :: skip_xfer
1711
1712         real :: delvol(2)
1713         real :: fracrem_num, fracrem_vol, fracxfr_num, fracxfr_vol
1714         real :: rbox_sv2(1:num_chem)
1715         real :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf
1716         real :: tmp_cwnumfrac, tmp_cwvolfrac
1717         real :: tmp_dpmeanvol, tmp_lnsg
1718         real :: xcut_num, xcut_vol
1719         real :: xdgnum_aaa(2), xlnsg(2)
1720         real :: xnum_aaa(2,3)
1721         real :: xvol_aaa(2,3)
1722
1723
1724 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1725 !   diagnostics for testing
1726         lunxx = -1
1727         if ((it==23) .and. (jt==1) .and. (kt<=20)) lunxx = 171
1728 #endif
1729
1730         p1st = param_first_scalar
1731
1732 !
1733 !   initial calculations for aitken (ii=1) and accum (ii=2) modes
1734 !
1735         skip_xfer = .false.
1736
1737         do ii = 1, 2
1738         itype = 1
1739         isize = ii
1740
1741 !   calculate new volumes for activated (jj=1), interstitial (jj=2), and combined (jj=3)
1742         tmpa = 0.0
1743         do ll = 1, ncomp_aer(itype)
1744             l = massptr_aer(ll,isize,itype,cw_phase)
1745             if (l >= p1st) tmpa = tmpa + rbox(l)/dens_aer(ll,itype)
1746         end do
1747         xvol_aaa(ii,1) = tmpa*1.0e-6
1748         xvol_aaa(ii,2) = xvol_old(ii,2)
1749         xnum_aaa(ii,1) = rbox(numptr_aer(isize,itype,cw_phase))
1750         xnum_aaa(ii,2) = rbox(numptr_aer(isize,itype,ai_phase))
1751
1752         xvol_aaa(ii,3) = xvol_aaa(ii,1) + xvol_aaa(ii,2)
1753         xnum_aaa(ii,3) = xnum_aaa(ii,1) + xnum_aaa(ii,2)
1754         delvol(ii) = xvol_aaa(ii,1) - xvol_old(ii,1)
1755
1756 !   check for negligible number or volume in aitken mode
1757         if (ii == 1) then
1758             if (xvol_aaa(ii,3) < smallvolaa) then
1759                 skip_xfer = .true.
1760                 exit
1761             end if
1762         end if
1763
1764 !   calculate dgnum for activated+interstitial combined using new volume
1765         if (xvol_aaa(ii,3) < smallvolaa) then
1766             tmp_dpmeanvol = dcen_sect(isize,itype)
1767         else
1768             tmp_dpmeanvol = xvol_aaa(ii,3)/xnum_aaa(ii,3)
1769             tmp_dpmeanvol = (tmp_dpmeanvol*6.0/pirs)**0.33333333
1770         end if
1771         xlnsg(ii) = log(sigmag_aer(isize,itype))
1772         xdgnum_aaa(ii) = tmp_dpmeanvol*exp(-1.5*xlnsg(ii)*xlnsg(ii))
1773 !   tmp_cwvolfrac = (cw volume)/(ai + cw volume)
1774         tmp_cwvolfrac = xvol_aaa(ii,1)/max(xvol_aaa(ii,3),1.e-30)
1775
1776 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1777 !   diagnostics for testing
1778         if (lunxx > 0) then
1779             if (ii == 1) write(lunxx,'(a)')
1780             write(lunxx,'(a,3i4,i8,2x,l2)') 'merge i,j,k,  ii,skip', &
1781                 it, jt, kt, ii, skip_xfer
1782             write(lunxx,'(a,1p,5e12.4)') 'cw  vol old/aaa, num aaa        ', &
1783                 xvol_old(ii,1), xvol_aaa(ii,1), xnum_aaa(ii,1)
1784             write(lunxx,'(a,1p,5e12.4)') 'ai  vol old/aaa, num aaa        ', &
1785                 xvol_old(ii,2), xvol_aaa(ii,2), xnum_aaa(ii,2)
1786             write(lunxx,'(a,1p,5e12.4)') 'a+c vol old/aaa, num aaa        ', &
1787                 xvol_old(ii,3), xvol_aaa(ii,3), xnum_aaa(ii,3)
1788             write(lunxx,'(a,1p,5e12.4)') 'cwnum/volfrac, dpmeanvol, dgnum ', &
1789                 (xnum_aaa(ii,1)/max(xnum_aaa(ii,3),1.e-30)), &
1790                 tmp_cwvolfrac, tmp_dpmeanvol, xdgnum_aaa(ii)
1791         end if
1792 #endif
1793
1794         end do  ! ii = 1, 2
1795
1796
1797 !   check for mode merging
1798         if ( skip_xfer ) return
1799         if (delvol(1) > delvol(2)) then
1800             continue
1801         else if ( (xdgnum_aaa(1) > 0.03e-4) .and. (xnum_aaa(1,3) > xnum_aaa(2,3)) ) then
1802             continue
1803         else
1804             return
1805         end if
1806 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1807         if (lunxx > 0) then
1808           if (delvol(1) > delvol(2)) then
1809             write(lunxx,'(a)') 'do merging - criterion 1'
1810           else if ( (xdgnum_aaa(1) > 0.03e-4) .and. (xnum_aaa(1,3) > xnum_aaa(2,3)) ) then
1811             write(lunxx,'(a)') 'do merging - criterion 2'
1812           end if
1813         end if
1814 #endif
1815
1816 !
1817 !   calc transfer fractions for volume/mass and number
1818 !   approach follows that in module_aerosols_sorgam (subr aerostep) except
1819 !   >>  the first steps of the calculation are done using the total
1820 !       (interstitial+activated) size distributions
1821 !   >>  the number and volume (moment-3) transfer amounts are then limited
1822 !       to the number/volume of the aitken activated distribution
1823 !
1824 !   xcut_num = [ln(dintsect/dgnuc)/xxlsgn], where dintsect is the diameter
1825 !   at which the aitken-mode and accum-mode number distribs intersect (overlap).
1826         tmpa = sqrt(2.0)
1827 !       aaa = getaf( nu0, ac0, dgnuc, dgacc, xxlsgn, xxlsga, sqrt2 )
1828         xcut_num = tmpa * getaf( xnum_aaa(1,3), xnum_aaa(2,3), &
1829                 xdgnum_aaa(1), xdgnum_aaa(2), xlnsg(1), xlnsg(2), tmpa )
1830
1831 !   forcing xcut_vol>0 means that no more than half of the aitken volume
1832 !   will be transferred
1833         tmpd = xcut_num
1834         tmpc = 3.0*xlnsg(1)
1835         xcut_vol = max( xcut_num-tmpc, 0.0 )
1836         xcut_num = xcut_vol + tmpc
1837         fracxfr_vol = norm01_uptail( xcut_vol )
1838         fracxfr_num = norm01_uptail( xcut_num )
1839         tmpe = fracxfr_vol ; tmpf = fracxfr_num
1840
1841         tmp_cwvolfrac = xvol_aaa(1,1)/max(xvol_aaa(1,3),1.e-30)
1842         tmp_cwnumfrac = xnum_aaa(1,1)/max(xnum_aaa(1,3),1.e-30)
1843         if ( (fracxfr_vol >= tmp_cwvolfrac) .or. &
1844              (fracxfr_num >= tmp_cwnumfrac) ) then
1845 !   limit volume fraction transferred to tmp_cwvolfrac
1846 !   limit number fraction transferred to tmp_cwnumfrac
1847             fracxfr_num = 1.0
1848             fracxfr_vol = 1.0
1849         else
1850 !   at this point, fracxfr_num/vol are fraction of
1851 !       interstitial+activated num/vol to be transferred
1852 !   convert them to fraction of activated num/vol to be transferred
1853             fracxfr_vol = fracxfr_vol/max(1.0e-10,tmp_cwvolfrac)
1854             fracxfr_num = fracxfr_num/max(1.0e-10,tmp_cwnumfrac)
1855 !   number fraction transferred cannot exceed volume fraction
1856             fracxfr_num = min( fracxfr_num, fracxfr_vol )
1857         end if
1858         fracrem_vol = 1.0 - fracxfr_vol
1859         fracrem_num = 1.0 - fracxfr_num
1860 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1861         if (lunxx > 0) then
1862             write(lunxx,'(a,1p,5e12.4)') 'xcut_num1/num2/vol              ', &
1863                 tmpd, xcut_num, xcut_vol
1864             write(lunxx,'(a,1p,5e12.4)') 'fracxfr_num3/vol3/num1,vol1     ', &
1865                 tmpf, tmpe, fracxfr_num, fracxfr_vol
1866         end if
1867 #endif
1868         if ( skip_xfer ) return
1869
1870 !   do the transfer
1871         rbox_sv2(:) = rbox(:)
1872         itype = 1
1873         isize_ait = 1
1874         isize_acc = 2
1875
1876         lfrm = numptr_aer(isize_ait,itype,cw_phase)
1877         ltoo = numptr_aer(isize_acc,itype,cw_phase)
1878         rbox(ltoo) = rbox(ltoo) + rbox(lfrm)*fracxfr_num
1879         rbox(lfrm) = rbox(lfrm)*fracrem_num
1880
1881         do ll = 1, ncomp_aer(itype)
1882             lfrm = massptr_aer(ll,isize_ait,itype,cw_phase)
1883             ltoo = massptr_aer(ll,isize_acc,itype,cw_phase)
1884             if (lfrm >= p1st) then
1885                 if (ltoo >= p1st) rbox(ltoo) = rbox(ltoo) &
1886                                              + rbox(lfrm)*fracxfr_vol
1887                 rbox(lfrm) = rbox(lfrm)*fracrem_vol
1888             end if
1889         end do
1890
1891
1892 #if defined ( ccboxtest_box_testing_active ) || defined ( cctemp_testing_active )
1893 !   more diagnostics for testing
1894         if (lunxx .gt. 0) then
1895             do ll = 1, 4
1896                 if (ll .eq. 1) then
1897                     lptr_dum(:,:) = lptr_so4_aer(:,:,cw_phase)
1898                 else if (ll .eq. 2) then
1899                     lptr_dum(:,:) = lptr_nh4_aer(:,:,cw_phase)
1900                 else if (ll .eq. 3) then
1901                     lptr_dum(:,:) = lptr_p25_aer(:,:,cw_phase)
1902                 else if (ll .eq. 4) then
1903                     lptr_dum(:,:) = numptr_aer(:,:,cw_phase)
1904                 end if
1905
1906                 if (ll .eq. 1) write(lunxx,'(a)')
1907                 write(lunxx,'(2a,i6,i3,2x,a)') 'sorgam_cloudchem_apply_mode_transfer',   &
1908                     ' - icase, ll', icase, ll,   &
1909                     chem_dname_table(id,lptr_dum(1,1))(1:12)
1910                 write(lunxx,'(a)') ' ty sz  rbox_sv1, rbox, rsub'
1911
1912                 itype = 1
1913                 do ii = 1, 2
1914                     if (ii == 1) then
1915                         isize = isize_ait
1916                     else
1917                         isize = isize_acc
1918                     end if
1919                     l = lptr_dum(isize,itype)
1920                     write(lunxx,'(2i3,1p,5e14.6)') &
1921                         itype, isize, rbox_sv1(l), rbox_sv2(l), rbox(l)
1922                 end do
1923             end do
1924
1925             if (icase .le. -5) then
1926                 write(*,*)   &
1927                 '*** stop in sorgam_cloudchem_apply_mode_transfer diags, icase =',   &
1928                 icase
1929                 call wrf_error_fatal('*** stop in sorgam_cloudchem_apply_mode_transfer diags')
1930             end if
1931         end if
1932 #endif
1933
1934
1935         return
1936         end subroutine sorgam_cloudchem_apply_mode_transfer
1937
1938
1939
1940 !-----------------------------------------------------------------------
1941         real function norm01_uptail( x )
1942 !
1943 !   norm01_uptail = cummulative pdf complement of normal(0,1) pdf
1944 !            = integral from x to +infinity of [normal(0,1) pdf]
1945 !
1946 !   erfc_num_recipes is from press et al, numerical recipes, 1990, page 164
1947 !
1948         implicit none
1949         real x, xabs
1950         real*8 erfc_approx, tmpa, t, z
1951
1952         xabs = abs(x)
1953         if (xabs >= 12.962359) then
1954             if (x > 0.0) then
1955                 norm01_uptail = 0.0
1956             else
1957                 norm01_uptail = 1.0
1958             end if
1959             return
1960         end if
1961
1962         z = xabs / sqrt(2.0_8)
1963         t = 1.0_8/(1.0_8 + 0.5_8*z)
1964
1965 !       erfc_approx =
1966 !     &   t*exp( -z*z - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
1967 !     &   t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
1968 !     &                                    t*(-1.13520398 +
1969 !     &   t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
1970
1971         tmpa =  ( -z*z - 1.26551223_8 + t*(1.00002368_8 + t*(0.37409196_8 +   &
1972           t*(0.09678418_8 + t*(-0.18628806_8 + t*(0.27886807_8 +   &
1973                                            t*(-1.13520398_8 +   &
1974           t*(1.48851587_8 + t*(-0.82215223_8 + t*0.17087277_8 )))))))))
1975
1976         erfc_approx = t * exp(tmpa)
1977         if (x .lt. 0.0) erfc_approx = 2.0_8 - erfc_approx
1978
1979         norm01_uptail = 0.5_8 * erfc_approx
1980
1981         return
1982         end function norm01_uptail
1983
1984 !-----------------------------------------------------------------------
1985         real function norm01_uptail_inv( x )
1986 !
1987 !   norm01_uptail_inv = inverse of norm01_uptail
1988 !       if y = norm01_uptail_inv( x ), then
1989 !       {integral from y to +infinity of [normal(0,1) pdf]} = x
1990 !   y is computed using newton's method
1991 !
1992         implicit none
1993
1994 !   fn parameters
1995         real x
1996
1997 !   local variables
1998         integer niter
1999         real dfdyinv, f, pi, sqrt2pi, tmpa, y, ynew
2000
2001         parameter (pi = 3.1415926535897932384626434)
2002
2003         if (x .le. 1.0e-38) then
2004             norm01_uptail_inv = 12.962359
2005             return
2006         else if (x .ge. 1.0) then
2007             norm01_uptail_inv = -12.962359
2008             return
2009         end if
2010
2011         sqrt2pi = sqrt( 2.0*pi )
2012
2013 !   initial guess
2014 !   crude
2015 !       y = 3.0*(0.5 - x)
2016 !   better
2017         tmpa = x
2018         tmpa = max( 0.0, min( 1.0, tmpa ) )
2019         tmpa = 4.0*tmpa*(1.0 - tmpa)
2020         tmpa = max( 1.0e-38, min( 1.0, tmpa ) )
2021         y = sqrt( -(pi/2.0)*log(tmpa) )
2022         if (x > 0.5) y = -y
2023
2024         f = norm01_uptail(y) - x
2025         do niter = 1, 100
2026 !   iterate - dfdy is computed analytically
2027             dfdyinv = -sqrt2pi * exp( 0.5*y*y )
2028             ynew = y - f*dfdyinv
2029             f = norm01_uptail(ynew) - x
2030
2031             if ( (ynew == y) .or.   &
2032                  (abs(f) <= abs(x)*1.0e-5) ) then
2033                 exit
2034             end if
2035             y = ynew
2036         end do
2037 9100    format( 'niter/x/f/y/ynew', i5, 4(1pe16.8) )
2038
2039         norm01_uptail_inv = ynew
2040         return
2041         end function norm01_uptail_inv
2042
2043
2044
2045 !-----------------------------------------------------------------------
2046         subroutine sorgam_cloudchem_dumpaa(   &
2047             id, ktau, ktauc, dtstepc, config_flags,   &
2048             p_phy, t_phy, rho_phy, alt,   &
2049             cldfra, ph_no2,   &
2050             moist, chem,   &
2051             gas_aqfrac, numgas_aqfrac,   &
2052             ids,ide, jds,jde, kds,kde,   &
2053             ims,ime, jms,jme, kms,kme,   &
2054             its,ite, jts,jte, kts,kte,   &
2055             qcldwtr_cutoff,   &
2056             itcur, jtcur, ktcur )
2057
2058         use module_state_description, only:   &
2059                 num_moist, num_chem, p_qc
2060         use module_scalar_tables, only:  chem_dname_table
2061         use module_configure, only:  grid_config_rec_type
2062         use module_data_sorgam_vbs
2063
2064
2065         implicit none
2066
2067 !   subr arguments
2068         integer, intent(in) ::   &
2069                 id, ktau, ktauc,   &
2070                 numgas_aqfrac,   &
2071                 ids, ide, jds, jde, kds, kde,   &
2072                 ims, ime, jms, jme, kms, kme,   &
2073                 its, ite, jts, jte, kts, kte,   &
2074                 itcur, jtcur, ktcur
2075 !   id - domain index
2076 !   ktau - time step number
2077 !   ktauc - gas and aerosol chemistry time step number
2078 !   numgas_aqfrac - last dimension of gas_aqfrac
2079
2080 !   [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for 'domain'
2081 !   [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for 'memory'
2082 !       Most arrays that are arguments to chem_driver
2083 !       are dimensioned with these spatial indices.
2084 !   [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for 'tile'
2085 !       chem_driver and routines under it do calculations
2086 !       over these spatial indices.
2087
2088         type(grid_config_rec_type), intent(in) :: config_flags
2089 !   config_flags - configuration and control parameters
2090
2091         real, intent(in) ::   &
2092             dtstepc, qcldwtr_cutoff
2093 !   dtstepc - time step for gas and aerosol chemistry(s)
2094
2095         real, intent(in),   &
2096                 dimension( ims:ime, kms:kme, jms:jme ) :: &
2097                 p_phy, t_phy, rho_phy, alt, cldfra, ph_no2
2098 !   p_phy - air pressure (Pa)
2099 !   t_phy - temperature (K)
2100 !   rho_phy - moist air density (kg/m^3)
2101 !   alt - dry air specific volume (m^3/kg)
2102 !   cldfra - cloud fractional area (0-1)
2103 !   ph_no2 - no2 photolysis rate (1/min)
2104
2105         real, intent(in),   &
2106                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
2107                 moist
2108 !   moist - mixing ratios of moisture species (water vapor,
2109 !       cloud water, ...) (kg/kg for mass species, #/kg for number species)
2110
2111         real, intent(inout),   &
2112                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
2113                 chem
2114 !   chem - mixing ratios of trace gas and aerosol species (ppm for gases,
2115 !       ug/kg for aerosol mass species, #/kg for aerosol number species)
2116
2117         real, intent(inout),   &
2118                 dimension( ims:ime, kms:kme, jms:jme, numgas_aqfrac ) :: &
2119                 gas_aqfrac
2120 !   gas_aqfrac - fraction (0-1) of gas that is dissolved in cloud water
2121
2122
2123 !   local variables
2124         integer :: it, jt, kt, l, ll, n
2125         integer :: isize, itype
2126
2127         real :: dumai, dumcw
2128         real :: qcldwtr
2129
2130
2131         it = itcur
2132         jt = jtcur
2133         kt = ktcur
2134
2135         write(*,*)
2136         write(*,*)
2137         write(*,*)
2138         write(*,9100)
2139         write(*,9102) ktau, it, jt, kt
2140 9100    format( 7('----------') )
2141 9102    format(   &
2142         'sorgam_cloudchem_dumpaa - ktau, i, j, k =', 4i5 )
2143
2144         do 2900 itype = 1, ntype_aer
2145         do 2900 isize = 1, nsize_aer(itype)
2146         if ( .not. do_cloudchem_aer(isize,itype) ) goto 2900
2147
2148         write(*,9110) isize
2149 9110    format( / 'isize, itype =', 2i3 /   &
2150         '  k    cldwtr      mass-ai   numb-ai      mass-cw   numb-cw' )
2151
2152         do 2800 kt = kte, kts, -1
2153
2154         dumai = 0.0
2155         dumcw = 0.0
2156         do ll = 1, ncomp_aer(itype)
2157             l = massptr_aer(ll,isize,itype,1)
2158             dumai = dumai + chem(it,kt,jt,l)
2159             l = massptr_aer(ll,isize,itype,2)
2160             dumcw = dumcw + chem(it,kt,jt,l)
2161         end do
2162         write(*,9120) kt,   &
2163                 moist(it,kt,jt,p_qc),   &
2164                 dumai, chem(it,kt,jt,numptr_aer(isize,itype,1)),   &
2165                 dumcw, chem(it,kt,jt,numptr_aer(isize,itype,2))
2166 9120    format( i3, 1p, e10.2, 2(3x, 2e10.2) )
2167
2168 2800    continue
2169 2900    continue
2170
2171         write(*,*)
2172         write(*,9100)
2173         write(*,*)
2174
2175 !   map from wrf-chem 3d arrays to pegasus clm & sub arrays
2176         kt = ktcur
2177         if ((ktau .eq. 30) .and. (it .eq. 23) .and.   &
2178               (jt .eq.  1) .and. (kt .eq. 11)) then
2179             qcldwtr = moist(it,kt,jt,p_qc)
2180             write(*,*)
2181             write(*,*)
2182             write(*,9102) ktau, it, jt, kt
2183             write(*,*)
2184             write( *, '(3(1pe10.2,3x,a))' )   &
2185                 (chem(it,kt,jt,l), chem_dname_table(id,l)(1:12), l=1,num_chem)
2186             write(*,*)
2187             write( *, '(3(1pe10.2,3x,a))' )   &
2188                   p_phy(it,kt,jt), 'p_phy     ',   &
2189                   t_phy(it,kt,jt), 't_phy     ',   &
2190                 rho_phy(it,kt,jt), 'rho_phy   ',   &
2191                     alt(it,kt,jt), 'alt       ',   &
2192                           qcldwtr, 'qcldwtr   ',   &
2193                    qcldwtr_cutoff, 'qcldwtrcut'
2194             write(*,*)
2195             write(*,9100)
2196             write(*,*)
2197         end if
2198
2199
2200         return
2201         end subroutine sorgam_cloudchem_dumpaa
2202
2203
2204
2205 !-----------------------------------------------------------------------
2206         end module module_sorgam_vbs_cloudchem