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