chem/module_mozcart_wetscav.F

   1
   2 MODULE module_mozcart_wetscav
   3
   4    USE module_configure
   5    USE module_state_description
   6
   7 IMPLICIT NONE
   8
   9 private
  10 public :: wetscav_mozcart_init
  11 public :: wetscav_mozcart
  12
  13 save
  14
  15    real, parameter :: zero = 0.
  16    real, parameter :: one  = 1.
  17    real, parameter :: four = 4.
  18    real(8), parameter :: oner8  = 1._8
  19    real(8), parameter :: fourr8 = 4._8
  20    real, parameter :: adj_factor = one + 10.*epsilon( one )
  21    real, parameter :: TICE = 273.
  22    real, parameter :: TMIX = 258.
  23    integer, parameter :: idiag = 0
  24    integer, parameter :: jdiag = 0
  25    integer, parameter :: kdiag = 18
  26
  27    integer :: hetcnt
  28    integer :: hno3_ndx = 0
  29    integer, allocatable :: wrf2tab(:)
  30    REAL, allocatable    :: mol_wght(:)
  31    logical, allocatable :: ice_uptake(:)
  32
  33    type wet_scav
  34      character(len=12) :: name
  35      integer :: p_ndx
  36      real    :: heff(6)
  37      real    :: molecw
  38      logical :: ice_uptake
  39      real    :: reteff
  40    end type wet_scav
  41
  42    type(wet_scav), allocatable :: wet_scav_tab(:)
  43
  44    LOGICAL, EXTERNAL  :: wrf_dm_on_monitor
  45
  46 CONTAINS
  47
  48    subroutine wetscav_mozcart_init( id, numgas, config_flags )
  49 !----------------------------------------------------------------------
  50 !  Initialize the mozart, mozcart wet scavenging module
  51 !----------------------------------------------------------------------
  52    use module_scalar_tables, only : chem_dname_table
  53    use module_chem_utilities, only : upcase
  54    use module_HLawConst, only : nHLC, HLC
  55
  56 !----------------------------------------------------------------------
  57 !  dummy arguments
  58 !----------------------------------------------------------------------
  59    integer, intent(in)                           :: id
  60    integer, intent(in)                           :: numgas
  61    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  62
  63 !----------------------------------------------------------------------
  64 !  local variables
  65 !----------------------------------------------------------------------
  66    integer :: m, m1, m2
  67    integer :: astat
  68    character(len=12)  :: wrf_spc_name
  69    character(len=64)  :: HL_tbl_name
  70    character(len=256) :: message
  71
  72 is_allocated : &
  73    if( .not. allocated( wet_scav_tab ) ) then
  74      call wrf_message( ' ' )
  75      call wrf_message( 'wetscav_mozcart_init: species names' )
  76      do m = param_first_scalar,numgas,6
  77        write(message,'(6a12)') chem_dname_table(id,m:min(m+5,numgas))
  78        call wrf_message( trim(message) )
  79      enddo
  80      call wrf_message( ' ' )
  81 !----------------------------------------------------------------------
  82 !  scan HLawConst table for match with chem_opt scheme gas phase species
  83 !  first just count matches
  84 !----------------------------------------------------------------------
  85      hetcnt = 0
  86      do m = 1,nHLC
  87        HL_tbl_name = HLC(m)%name
  88        call upcase( HL_tbl_name )
  89        do m1 = param_first_scalar,numgas
  90          wrf_spc_name = chem_dname_table(id,m1)
  91          call upcase( wrf_spc_name )
  92          if( trim(HL_tbl_name) == trim(wrf_spc_name) ) then
  93            hetcnt = hetcnt + 1
  94            exit
  95          endif
  96        end do
  97      end do
  98
  99      if( hetcnt > 0 ) then
 100        allocate( wet_scav_tab(hetcnt),stat=astat )
 101        if( astat /= 0 ) then
 102          call wrf_error_fatal("mozcart_wetscav_init: failed to allocate wet_scav_tab")
 103        endif
 104      else
 105        call wrf_message( ' ' )
 106        call wrf_message( 'wetscav_mozcart_init: There are no wet scavenging mozcart species' )
 107        call wrf_message( ' ' )
 108        return
 109      endif
 110
 111 !----------------------------------------------------------------------
 112 !  now put matches in local table
 113 !----------------------------------------------------------------------
 114      m2 = 0
 115      do m = 1,nHLC
 116        HL_tbl_name = HLC(m)%name
 117        call upcase( HL_tbl_name )
 118        do m1 = param_first_scalar,numgas
 119          wrf_spc_name = chem_dname_table(id,m1)
 120          call upcase( wrf_spc_name )
 121          if( trim(HL_tbl_name) == trim(wrf_spc_name) ) then
 122            m2 = m2 + 1
 123            wet_scav_tab(m2)%name  = chem_dname_table(id,m1)
 124            wet_scav_tab(m2)%p_ndx = m1
 125            wet_scav_tab(m2)%molecw = HLC(m)%mw
 126            wet_scav_tab(m2)%heff(1:6) = HLC(m)%hcnst(1:6)
 127            exit
 128          endif
 129        end do
 130      end do
 131
 132      wet_scav_tab(:)%ice_uptake = .false.
 133      wet_scav_tab(:)%reteff = 0.
 134      do m = 1,hetcnt
 135        if( wet_scav_tab(m)%p_ndx == p_h2o2 ) then
 136          wet_scav_tab(m)%ice_uptake = .true.
 137          wet_scav_tab(m)%reteff     = .64
 138        elseif( wet_scav_tab(m)%p_ndx == p_hno3 ) then
 139          wet_scav_tab(m)%ice_uptake = .true.
 140          wet_scav_tab(m)%reteff     = 1.
 141          hno3_ndx = m
 142        elseif( wet_scav_tab(m)%p_ndx == p_hcooh ) then
 143          wet_scav_tab(m)%ice_uptake = .true.
 144          wet_scav_tab(m)%reteff     = .64
 145        elseif( wet_scav_tab(m)%p_ndx == p_ch3ooh ) then
 146          wet_scav_tab(m)%ice_uptake = .true.
 147          wet_scav_tab(m)%reteff     = .02
 148        elseif( wet_scav_tab(m)%p_ndx == p_so2 ) then
 149          wet_scav_tab(m)%ice_uptake = .true.
 150          wet_scav_tab(m)%reteff     = .02
 151        elseif( wet_scav_tab(m)%p_ndx == p_hcooh ) then
 152          wet_scav_tab(m)%ice_uptake = .true.
 153          wet_scav_tab(m)%reteff     = .68
 154        endif
 155      end do
 156
 157      allocate( wrf2tab(hetcnt),mol_wght(hetcnt),ice_uptake(hetcnt),stat=astat )
 158      if( astat /= 0 ) then
 159        call wrf_error_fatal("mozcart_wetscav_init: failed to allocate wrf2tab,mol_wght,ice_uptake")
 160      endif
 161      do m = 1,hetcnt
 162        wrf2tab(m)    = m
 163        mol_wght(m)   = wet_scav_tab(m)%molecw
 164        ice_uptake(m) = wet_scav_tab(m)%ice_uptake
 165      end do
 166
 167      call wrf_message( 'wetscav_mozcart_init: Wet scavenging mozcart species' )
 168      do m = 1,hetcnt
 169        write(message,*) m,wrf2tab(m),trim(wet_scav_tab(wrf2tab(m))%name),mol_wght(m),ice_uptake(m)
 170        call wrf_message( trim(message) )
 171      end do
 172
 173      if( wrf_dm_on_monitor() ) then
 174        call wrf_debug( 100,' ' )
 175        write(message,*) 'wetscav_mozcart_init: hetcnt = ',hetcnt
 176        call wrf_debug( 100, trim(message) )
 177        write(message,*) 'wetscav_mozcart_init: hno3_ndx = ',hno3_ndx
 178        call wrf_debug( 100, trim(message) )
 179        call wrf_debug( 100,' ' )
 180      endif
 181    endif is_allocated
 182
 183    end subroutine wetscav_mozcart_init
 184
 185    subroutine wetscav_mozcart( id, ktau, dtstep, ktauc, config_flags,                      &
 186                                dtstepc, t_phy, p8w, t8w, p_phy,                            &
 187                                chem, rho_phy, cldfra, rainprod, evapprod,                  &
 188                                qv_b4mp, qc_b4mp, qi_b4mp, qs_b4mp,                         &
 189                                gas_aqfrac, numgas_aqfrac, dz8w, dx, dy,                    &
 190                                qv, qc, qi, qs,                                             &
 191 ! 20131125 acd_ck_washout start
 192 !                               hno3_col_mdel,                                              &
 193                                delta_mass_col,                                             &
 194 ! 20131125 acd_ck_washout end
 195                                ids,ide, jds,jde, kds,kde,                                  &
 196                                ims,ime, jms,jme, kms,kme,                                  &
 197                                its,ite, jts,jte, kts,kte                                   )
 198
 199    USE module_configure, only: grid_config_rec_type
 200    USE module_state_description
 201    USE module_model_constants, only: g, mwdry
 202
 203 !----------------------------------------------------------------------
 204 !  dummy arguments
 205 !----------------------------------------------------------------------
 206    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
 207
 208    INTEGER,      INTENT(IN   )    ::                                &
 209                                       ids,ide, jds,jde, kds,kde,    &
 210                                       ims,ime, jms,jme, kms,kme,    &
 211                                       its,ite, jts,jte, kts,kte,    &
 212                                       id, ktau, ktauc, numgas_aqfrac
 213    REAL, INTENT(IN   ) :: dtstep, dtstepc
 214    REAL, INTENT(IN   ) :: dx, dy
 215 !----------------------------------------------------------------------
 216 ! all advected chemical species
 217 !----------------------------------------------------------------------
 218    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),          &
 219          INTENT(INOUT ) ::                                chem
 220 !----------------------------------------------------------------------
 221 ! fraction of gas species in cloud water
 222 !----------------------------------------------------------------------
 223    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, numgas_aqfrac ),     &
 224          INTENT(IN ) ::                                   gas_aqfrac
 225 !----------------------------------------------------------------------
 226 ! input from meteorology
 227 !----------------------------------------------------------------------
 228    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,        &
 229          INTENT(IN   ) ::                                           &
 230                                                       t_phy,        &
 231                                                       p_phy,        &
 232                                                       t8w,          &
 233                                                       p8w,          &
 234                                                       dz8w,         &
 235                                                     rho_phy
 236
 237    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::    &
 238                                                                  QV, &
 239                                                                  QI, &
 240                                                                  QC, &
 241                                                                  QS
 242    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::    &
 243                                                                  rainprod, &
 244                                                                  evapprod
 245    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::    &
 246                                                                  qv_b4mp, &
 247                                                                  qc_b4mp, &
 248                                                                  qi_b4mp, &
 249                                                                  qs_b4mp
 250
 251    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,        &
 252          INTENT(INOUT ) ::                                cldfra
 253
 254 ! 20131125 acd_ck_washout start
 255 !   REAL,  DIMENSION( ims:ime , jms:jme )         ,                  &
 256 !         INTENT(INOUT ) ::                                hno3_col_mdel
 257
 258    REAL,  DIMENSION( ims:ime , jms:jme, num_chem )         ,           &
 259          INTENT(OUT ) ::                                delta_mass_col
 260 ! 20131125 acd_ck_washout end
 261
 262 !----------------------------------------------------------------------
 263 !  local variables
 264 !----------------------------------------------------------------------
 265    real, parameter    :: t0 = 298.
 266 !  real, parameter    :: ph = 1.e-5
 267 !  real, parameter    :: hp_inv = 1./ph
 268    real, parameter    :: hion = 1.e-5
 269    real, parameter    :: hion_inv = 1./hion
 270    real, parameter    :: henry_thres = 1.e4
 271
 272    integer :: i, j, k, ktem1, m, m1
 273    integer :: pndx
 274    integer :: cld_col_cnt
 275    integer :: precip_col_cnt
 276    integer :: max_ndx(3)
 277    integer :: wrk_ndx(3)
 278    REAL :: area
 279    REAL :: e298, dhr
 280    REAL :: percent_cld
 281    REAL :: percent_precip
 282    REAL :: max_rls
 283    REAL :: layer_mass(kts:kte)
 284    REAL :: delp(kts:kte)
 285    REAL :: p(kts:kte)
 286    REAL :: t(kts:kte)
 287    REAL :: rls(kts:kte)
 288    REAL :: evaprate(kts:kte)
 289    REAL :: totprec(kts:kte)
 290    REAL :: totevap(kts:kte)
 291    REAL :: wrk(kts:kte)
 292    REAL :: tfac(kts:kte)
 293    REAL :: dk1s(kts:kte)
 294    REAL :: dk2s(kts:kte)
 295    REAL :: kh(kts:kte)
 296    REAL :: diff(its:ite,kts:kte,jts:jte)
 297    REAL :: hno3(kts:kte)
 298    REAL :: wrk_mass(hetcnt)
 299    REAL :: trc_mass(kts:kte,hetcnt)
 300    REAL :: heff(kts:kte,hetcnt)
 301
 302    logical :: is_hno3
 303    logical :: tckaqb(hetcnt)
 304
 305 !++mmb
 306    REAL :: reteff(hetcnt)
 307 !--mmb
 308
 309    character(len=128) :: message
 310
 311 has_wet_scav : &
 312    if( hetcnt > 0 ) then
 313 !----------------------------------------------------------------------
 314 !  form cloud fraction
 315 !----------------------------------------------------------------------
 316      CALL cal_cldfra3( CLDFRA, qc_b4mp, qi_b4mp, qs_b4mp,         &
 317                       ids,ide, jds,jde, kds,kde,                  &
 318                       ims,ime, jms,jme, kms,kme,                  &
 319                       its,ite, jts,jte, kts,kte                   )
 320
 321 !----------------------------------------------------------------------
 322 !  washout soluble species
 323 !----------------------------------------------------------------------
 324      ktem1 = kte - 1
 325      area = dx * dy
 326      diff(:,:,:) = 0.
 327      cld_col_cnt = 0
 328      precip_col_cnt = 0
 329 ! 20131125 acd_ck_washout start
 330 !     hno3_col_mdel(:,:) = 0.
 331      delta_mass_col(:,:,:) = 0.
 332 ! 20131125 acd_ck_washout end
 333      max_rls       = 0.
 334 jloop : &
 335      do j = jts,jte
 336 iloop : &
 337        do i = its,ite
 338            t(kts:kte)      = t_phy(i,kts:kte,j)
 339            tfac(kts:ktem1) = (t0 - t(kts:ktem1))/(t0*t(kts:ktem1))
 340            p(kts:kte)      = p_phy(i,kts:kte,j)*.01
 341            delp(kts:ktem1) = p8w(i,kts:ktem1,j) - p8w(i,kts+1:kte,j)
 342            layer_mass(kts:ktem1) = area*delp(kts:ktem1)/g
 343            totprec(kts:ktem1) = rainprod(i,kts:ktem1,j)*layer_mass(kts:ktem1)
 344            totevap(kts:ktem1) = evapprod(i,kts:ktem1,j)*layer_mass(kts:ktem1)
 345            rls(kte)      = 0.
 346            evaprate(kte) = 0.
 347            do k = ktem1,kts,-1
 348              rls(k) = max( 0.,totprec(k)-totevap(k)+rls(k+1) )
 349              evaprate(k) = min( 1.,totevap(k)/(rls(k+1)+1.e-20) )
 350            end do
 351 column_has_precip : &
 352          if( any( rls(kts:ktem1) > 0. ) ) then
 353            if( maxval(rls(kts:ktem1)) >= max_rls ) then
 354              max_rls = max( max_rls,maxval(rls(kts:ktem1)) )
 355              max_ndx(3:3) = maxloc(rls(kts:ktem1))
 356              max_ndx(1:2) = (/ i,j /)
 357              if( max_ndx(3) /= kts ) then
 358                write(message,'(''wetscav: max rls not at srf; time,i,j,k = '',4i6)') ktau,max_ndx(:)
 359                call wrf_debug( 100,trim(message) )
 360              endif
 361            endif
 362            precip_col_cnt = precip_col_cnt + 1
 363 species_loop : &
 364            do m = 1,hetcnt
 365              m1 = wrf2tab(m)
 366              pndx = wet_scav_tab(m1)%p_ndx
 367              if( pndx == p_hno3 ) then
 368                hno3(kts:kte)   = chem(i,kts:kte,j,p_hno3)
 369              endif
 370              wrk(kts:ktem1)        = 1.e-6*mol_wght(m)*chem(i,kts:ktem1,j,pndx)/mwdry
 371              trc_mass(kts:ktem1,m) = wrk(kts:ktem1)*layer_mass(kts:ktem1)
 372              wrk_mass(m)    = sum( trc_mass(kts:ktem1,m) )
 373              e298 = wet_scav_tab(m1)%heff(1)
 374              dhr  = wet_scav_tab(m1)%heff(2)
 375              kh(kts:ktem1) = e298 * exp( dhr*tfac(kts:ktem1) )
 376              e298 = wet_scav_tab(m1)%heff(3)
 377              dhr  = wet_scav_tab(m1)%heff(4)
 378              if( e298 /= 0. ) then
 379                dk1s(kts:ktem1) = e298*exp( dhr*tfac(kts:ktem1) )
 380              else
 381                dk1s(kts:ktem1) = 0.
 382              endif
 383              e298 = wet_scav_tab(m1)%heff(5)
 384              dhr  = wet_scav_tab(m1)%heff(6)
 385              if( e298 /= 0. ) then
 386                dk2s(kts:ktem1) = e298*exp( dhr*tfac(kts:ktem1) )
 387              else
 388                dk2s(kts:ktem1) = 0.
 389              endif
 390              if( pndx /= p_nh3 ) then
 391                heff(kts:ktem1,m) = kh(kts:ktem1)*(1. + dk1s(kts:ktem1)*hion_inv*(1. + dk2s(kts:ktem1)*hion_inv))
 392              else
 393 !----------------------------------------------------------------------
 394 !  special handling for nh3
 395 !----------------------------------------------------------------------
 396                heff(kts:ktem1,m) = kh(kts:ktem1)*(1. + dk1s(kts:ktem1)*hion/dk2s(kts:ktem1))
 397              endif
 398              tckaqb(m) = any( heff(kts:ktem1,m) > henry_thres )
 399 !++mmb
 400              reteff(m) = wet_scav_tab(m1)%reteff
 401 !--mmb
 402            end do species_loop
 403
 404 !----------------------------------------------------------------------
 405 !  jneu washout
 406 !----------------------------------------------------------------------
 407            CALL washout( kte-kts+1, hetcnt, dtstep, trc_mass, layer_mass, &
 408                          p, dz8w(i,kts:kte,j), rls, qc_b4mp(i,kts:kte,j), qi_b4mp(i,kts:kte,j), &
 409                          cldfra(i,kts:kte,j), t, evaprate, area, heff, &
 410 !++mmb
 411 !                        mol_wght, tckaqb, ice_uptake, i, j )
 412                          mol_wght, tckaqb, ice_uptake, i, j, reteff )
 413 !--mmb
 414
 415 species_loop1 : &
 416            do m = 1,hetcnt
 417              m1 = wrf2tab(m)
 418              pndx = wet_scav_tab(m1)%p_ndx
 419              is_hno3 = pndx == p_hno3
 420 ! 20131125 acd_ck_washout start
 421 !             if( is_hno3 ) then
 422 !              hno3_col_mdel(i,j) = sum( trc_mass(kts:ktem1,m) ) - wrk_mass(m)
 423 !            endif
 424              delta_mass_col(i,j,pndx) = sum( trc_mass(kts:ktem1,m) ) - wrk_mass(m)
 425 ! 20131125 acd_ck_washout end
 426
 427              wrk(kts:ktem1) = 1.e6*mwdry*trc_mass(kts:ktem1,m)/mol_wght(m)
 428              chem(i,kts:ktem1,j,pndx) = wrk(kts:ktem1)/layer_mass(kts:ktem1)
 429              if( is_hno3 ) then
 430                diff(i,kts:ktem1,j) = 100.*(chem(i,kts:ktem1,j,p_hno3) - hno3(kts:ktem1))/hno3(kts:ktem1)
 431              endif
 432            end do species_loop1
 433          endif column_has_precip
 434        end do iloop
 435      end do jloop
 436
 437      write(message,'(''washout: max rls @ (i,j,k) '',3i4,'' = '',1pg15.7)') max_ndx(:),max_rls
 438      call wrf_debug( 100,trim(message) )
 439
 440      percent_precip = 100.*real(precip_col_cnt)/real((ite-its+1)*(jte-jts+1))
 441      write(*,*) 'wetscav_mozcart: percent columns with precip = ',percent_precip,'%'
 442    endif has_wet_scav
 443
 444    end subroutine wetscav_mozcart
 445
 446    SUBROUTINE cal_cldfra( CLDFRA,QC,QI,                               &
 447                           ids,ide, jds,jde, kds,kde,                  &
 448                           ims,ime, jms,jme, kms,kme,                  &
 449                           its,ite, jts,jte, kts,kte                   )
 450 !---------------------------------------------------------------------
 451 ! !DESCRIPTION:
 452 ! Compute cloud fraction from input ice and cloud water fields
 453 ! if provided.
 454 !
 455 ! Whether QI or QC is active or not is determined from the indices of
 456 ! the fields into the 4D scalar arrays in WRF. These indices are
 457 ! P_QI and P_QC, respectively, and they are passed in to the routine
 458 ! to enable testing to see if QI and QC represent active fields in
 459 ! the moisture 4D scalar array carried by WRF.
 460 !
 461 ! If a field is active its index will have a value greater than or
 462 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
 463 ! this routine.
 464 !EOP
 465 !---------------------------------------------------------------------
 466
 467    IMPLICIT NONE
 468
 469    INTEGER,  INTENT(IN   )   ::           ids,ide, jds,jde, kds,kde, &
 470                                           ims,ime, jms,jme, kms,kme, &
 471                                           its,ite, jts,jte, kts,kte
 472
 473    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT  ) ::    &
 474                                                              CLDFRA
 475
 476    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::    &
 477                                                                  QI, &
 478                                                                  QC
 479 !----------------------------------------------------------------------
 480 !  local variables
 481 !----------------------------------------------------------------------
 482 !  REAL, parameter :: thresh = 1.e-6
 483    REAL, parameter :: thresh = 1.e-9
 484
 485    INTEGER :: j,k
 486
 487    DO j = jts,jte
 488      DO k = kts,kte-1
 489        where( (qc(its:ite,k,j) + qi(its:ite,k,j)) > thresh )
 490          cldfra(its:ite,k,j) = one
 491        elsewhere
 492          cldfra(its:ite,k,j) = zero
 493        endwhere
 494      ENDDO
 495    ENDDO
 496
 497    END SUBROUTINE cal_cldfra
 498
 499    SUBROUTINE cal_cldfra3( CLDFRA,QC,QI, QS,                          &
 500                           ids,ide, jds,jde, kds,kde,                  &
 501                           ims,ime, jms,jme, kms,kme,                  &
 502                           its,ite, jts,jte, kts,kte                   )
 503 !---------------------------------------------------------------------
 504 ! !DESCRIPTION:
 505 ! Compute cloud fraction from input ice and cloud water fields
 506 ! if provided.
 507 !
 508 ! Whether QI or QC is active or not is determined from the indices of
 509 ! the fields into the 4D scalar arrays in WRF. These indices are
 510 ! P_QI and P_QC, respectively, and they are passed in to the routine
 511 ! to enable testing to see if QI and QC represent active fields in
 512 ! the moisture 4D scalar array carried by WRF.
 513 !
 514 ! If a field is active its index will have a value greater than or
 515 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
 516 ! this routine.
 517 !EOP
 518 !---------------------------------------------------------------------
 519
 520    IMPLICIT NONE
 521
 522    INTEGER,  INTENT(IN   )   ::           ids,ide, jds,jde, kds,kde, &
 523                                           ims,ime, jms,jme, kms,kme, &
 524                                           its,ite, jts,jte, kts,kte
 525
 526    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT  ) ::    &
 527                                                              CLDFRA
 528
 529    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::    &
 530                                                                  QI, &
 531                                                                  QC, &
 532                                                                  QS
 533 !----------------------------------------------------------------------
 534 !  local variables
 535 !----------------------------------------------------------------------
 536    REAL, parameter :: thresh = 1.e-9
 537
 538    INTEGER :: j
 539
 540    DO j = jts,jte
 541      where( (qc(its:ite,kts:kte,j) + qi(its:ite,kts:kte,j) + qs(its:ite,kts:kte,j)) > thresh )
 542        cldfra(its:ite,kts:kte,j) = one
 543      elsewhere
 544        cldfra(its:ite,kts:kte,j) = zero
 545      endwhere
 546    ENDDO
 547
 548    END SUBROUTINE cal_cldfra3
 549
 550    SUBROUTINE cal_cldfra2( CLDFRA, QV, QC, QI, QS,                     &
 551                            t_phy, p_phy,                               &
 552                            ids,ide, jds,jde, kds,kde,                  &
 553                            ims,ime, jms,jme, kms,kme,                  &
 554                            its,ite, jts,jte, kts,kte                   )
 555
 556 !----------------------------------------------------------------------
 557 !  dummy arguments
 558 !----------------------------------------------------------------------
 559    INTEGER,  INTENT(IN   )   ::           ids,ide, jds,jde, kds,kde, &
 560                                           ims,ime, jms,jme, kms,kme, &
 561                                           its,ite, jts,jte, kts,kte
 562
 563    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT  ) ::    &
 564                                                              CLDFRA
 565
 566    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::    &
 567                                                                  QV, &
 568                                                                  QI, &
 569                                                                  QC, &
 570                                                                  QS, &
 571                                                               t_phy, &
 572                                                               p_phy
 573
 574
 575 !----------------------------------------------------------------------
 576 !  local variables
 577 !----------------------------------------------------------------------
 578    integer, parameter  :: idbg = 144, kdbg = 15, jdbg = 147
 579    REAL    , PARAMETER :: ALPHA0 = 100.
 580    REAL    , PARAMETER :: GAMMA = 0.49
 581    REAL    , PARAMETER :: QCLDMIN = 1.E-12
 582    REAL    , PARAMETER :: PEXP = 0.25
 583    REAL    , PARAMETER :: RHGRID =1.0
 584    REAL    , PARAMETER :: SVP1 = 0.61078
 585    REAL    , PARAMETER :: SVP2 = 17.2693882
 586    REAL    , PARAMETER :: SVPI2 = 21.8745584
 587    REAL    , PARAMETER :: SVP3 = 35.86
 588    REAL    , PARAMETER :: SVPI3 = 7.66
 589    REAL    , PARAMETER :: SVPT0 = 273.15
 590    REAL    , PARAMETER :: r_d = 287.
 591    REAL    , PARAMETER :: r_v = 461.6
 592    REAL    , PARAMETER :: ep_2 = r_d/r_v
 593
 594    INTEGER :: i,j,k
 595    INTEGER :: imax, jmax, kmax
 596    REAL    :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight
 597    REAL    :: QCLD, DENOM, ARG, SUBSAT, wrk
 598    REAL    :: relhum_max, wrk_max
 599
 600
 601 ! !DESCRIPTION:
 602 !----------------------------------------------------------------------
 603 ! Compute cloud fraction from input ice and cloud water fields
 604 ! if provided.
 605 !
 606 ! Whether QI or QC is active or not is determined from the indices of
 607 ! the fields into the 4D scalar arrays in WRF. These indices are
 608 ! P_QI and P_QC, respectively, and they are passed in to the routine
 609 ! to enable testing to see if QI and QC represent active fields in
 610 ! the moisture 4D scalar array carried by WRF.
 611 !
 612 ! If a field is active its index will have a value greater than or
 613 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
 614 ! this routine.
 615 !----------------------------------------------------------------------
 616 !EOP
 617
 618
 619 !-----------------------------------------------------------------------
 620 !     COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
 621 !     (modified by Ferrier, Feb '02)
 622 !
 623 !     Cloud fraction parameterization follows Randall, 1994
 624 !     (see Hong et al., 1998)
 625 !-----------------------------------------------------------------------
 626 ! Note: ep_2=287./461.6 Rd/Rv
 627 ! Note: R_D=287.
 628
 629 ! Alternative calculation for critical RH for grid saturation
 630 !     RHGRID=0.90+.08*((100.-DX)/95.)**.5
 631
 632 ! Calculate saturation mixing ratio weighted according to the fractions of
 633 ! water and ice.
 634 ! Following:
 635 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure''  J. Appl. Meteor.  6 p.204
 636 !    es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
 637 !
 638 !       over ice        over water
 639 ! a =   21.8745584      17.2693882
 640 ! b =   7.66            35.86
 641 !---------------------------------------------------------------------
 642
 643     CLDFRA(its:ite,kts:kte,jts:jte) = 0.
 644     relhum_max = -100.
 645     wrk_max    = -10000.
 646     imax = 0; kmax = 0; jmax = 0
 647
 648     DO j = jts,jte
 649       DO k = kts,kte
 650         DO i = its,ite
 651 !---------------------------------------------------------------------
 652 !    Determine cloud fraction (modified from original algorithm)
 653 !---------------------------------------------------------------------
 654           QCLD = QI(i,k,j) + QC(i,k,j) + QS(i,k,j)
 655 has_cloud : &
 656           IF( QCLD >= QCLDMIN ) THEN
 657             tc   = t_phy(i,k,j) - SVPT0
 658             esw  = 1000.0 * SVP1 * EXP( SVP2  * tc / ( t_phy(i,k,j) - SVP3  ) )
 659             esi  = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
 660             QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
 661             QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
 662
 663             weight     = (QI(i,k,j) + QS(i,k,j)) / QCLD
 664             QVS_WEIGHT = (1. - weight)*QVSW + weight*QVSI
 665             RHUM       = QV(i,k,j)/QVS_WEIGHT              !--- Relative humidity
 666 !---------------------------------------------------------------------
 667 !    Assume zero cloud fraction if there is no cloud mixing ratio
 668 !---------------------------------------------------------------------
 669             IF( RHUM >= RHGRID )THEN
 670 !---------------------------------------------------------------------
 671 !    Assume cloud fraction of unity if near saturation and the cloud
 672 !    mixing ratio is at or above the minimum threshold
 673 !---------------------------------------------------------------------
 674               CLDFRA(i,k,j) = 1.
 675             ELSE
 676 !---------------------------------------------------------------------
 677 !    Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
 678 !    modified based on assumed grid-scale saturation at RH=RHgrid.
 679 !---------------------------------------------------------------------
 680               SUBSAT = MAX( 1.E-10,RHGRID*QVS_WEIGHT - QV(i,k,j) )
 681               DENOM  = SUBSAT**GAMMA
 682               ARG    = MAX( -6.9,-ALPHA0*QCLD/DENOM )    ! <-- EXP(-6.9)=.001
 683 !---------------------------------------------------------------------
 684 !   prevent negative values  (new)
 685 !---------------------------------------------------------------------
 686               RHUM = MAX( 1.E-10, RHUM )
 687               wrk  = (RHUM/RHGRID)**PEXP*(1. - EXP( ARG ))
 688               if( rhum >= relhum_max ) then
 689                 relhum_max = rhum
 690                 imax = i
 691                 kmax = k
 692                 jmax = j
 693               endif
 694               IF( wrk >= .01 ) then
 695                 CLDFRA(i,k,j) = wrk
 696                 if( wrk >= wrk_max ) then
 697                   wrk_max = wrk
 698                 endif
 699               ENDIF
 700             ENDIF
 701           ENDIF has_cloud
 702         END DO
 703       END DO
 704     END DO
 705
 706    END SUBROUTINE cal_cldfra2
 707
 708    subroutine WASHOUT( LPAR, NTRACE, DTSCAV, QTTJFL, QM, &
 709                        POFL, DELZ, RLS, CLWC, CIWC, &
 710                        CFR, TEM, EVAPRATE, GAREA, HSTAR, &
 711 !++mmb
 712 !                      TCMASS, TCKAQB, TCNION, ii, jj )
 713                        TCMASS, TCKAQB, TCNION, ii, jj, RETEFF )
 714 !--mmb
 715 !-----------------------------------------------------------------------
 716 !---p-conde 5.4 (2007)   -----called from main-----
 717 !---called from pmain to calculate rainout and washout of tracers
 718 !---revised by JNEU 8/2007
 719 !---
 720 !-LAER has been removed - no scavenging for aerosols
 721 !-LAER could be used as LWASHTYP
 722 !---WILL THIS WORK FOR T42->T21???????????
 723 !-----------------------------------------------------------------------
 724
 725 !-----------------------------------------------------------------------
 726 !  dummy arguments
 727 !-----------------------------------------------------------------------
 728       integer, intent(in)  :: LPAR
 729       integer, intent(in)  :: NTRACE
 730       integer, intent(in)  :: ii, jj
 731       real,  intent(in)    :: DTSCAV
 732       real,  intent(in)    :: GAREA
 733       real,  intent(in)    :: QM(LPAR)
 734       real,  intent(in)    :: POFL(LPAR)
 735       real,  intent(in)    :: DELZ(LPAR)
 736       real,  intent(in)    :: RLS(LPAR)
 737       real,  intent(in)    :: CLWC(LPAR)
 738       real,  intent(in)    :: CIWC(LPAR)
 739       real,  intent(in)    :: CFR(LPAR)
 740       real,  intent(in)    :: TEM(LPAR)
 741       real,  intent(in)    :: EVAPRATE(LPAR)
 742       real,  intent(in)    :: TCMASS(NTRACE)
 743       real,  intent(in)    :: HSTAR(LPAR,NTRACE)
 744       real,  intent(inout) :: QTTJFL(LPAR,NTRACE)
 745       logical, intent(in)  :: TCKAQB(NTRACE)
 746       logical, intent(in)  :: TCNION(NTRACE)
 747 !++mmb
 748       real,  intent(in)    :: RETEFF(NTRACE)
 749 !--mmb
 750
 751 !-----------------------------------------------------------------------
 752 !  local variables
 753 !-----------------------------------------------------------------------
 754       integer  :: I, J, L, LE, LM1, N
 755       integer  :: LWASHTYP, LICETYP
 756
 757       real :: WRK, RNEW_TST
 758       real :: CLWX
 759       real :: RNEW, RPRECIP, DELTARIMEMASS, DELTARIME, RAMPCT
 760       real :: MASSLOSS
 761       real :: DOR, DNEW, DEMP, COLEFFSNOW, RHOSNOW
 762       real :: WEMP, REMP, RRAIN, RWASH
 763       real :: QTPRECIP, QTRAIN, QTCXA, QTAX
 764
 765       real :: FAMA, RAMA, DAMA, FCA, RCA, DCA
 766       real :: FAX, RAX, DAX, FCXA, RCXA, DCXA, FCXB, RCXB, DCXB
 767       real :: RAXADJ, FAXADJ, RAXADJF
 768       real :: QTDISCF, QTDISRIME, QTDISCXA
 769       real :: QTEVAPAXP, QTEVAPAXW, QTEVAPAX
 770       real :: QTWASHAX
 771       real :: QTEVAPCXAP, QTEVAPCXAW, QTEVAPCXA
 772       real :: QTWASHCXA, QTRIMECXA
 773       real :: QTRAINCXA, QTRAINCXB
 774       real :: QTTOPCA, QTTOPAA, QTTOPCAX, QTTOPAAX
 775
 776       real :: AMPCT, AMCLPCT, CLNEWPCT, CLNEWAMPCT, CLOLDPCT, CLOLDAMPCT
 777       real :: RAXLOC, RCXALOC, RCXBLOC, RCALOC, RAMALOC, RCXPCT
 778
 779       real :: QTNETLCXA, QTNETLCXB, QTNETLAX, QTNETL
 780       real :: QTDISSTAR
 781
 782       real :: CFXX(lpar)
 783       real :: QTT(lpar)
 784       real :: QTTNEW(lpar)
 785       real :: rls_wrk(lpar)
 786       real :: rnew_wrk(lpar)
 787       real :: rca_wrk(lpar)
 788       real :: fca_wrk(lpar)
 789       real :: rcxa_wrk(lpar)
 790       real :: fcxa_wrk(lpar)
 791       real :: rcxb_wrk(lpar)
 792       real :: fcxb_wrk(lpar)
 793       real :: rax_wrk(lpar,2)
 794       real :: fax_wrk(lpar,2)
 795       real :: rama_wrk(lpar)
 796       real :: fama_wrk(lpar)
 797       real :: deltarime_wrk(lpar)
 798       real :: clwx_wrk(lpar)
 799       real :: frc(lpar,3)
 800       real :: rlsog(lpar)
 801
 802       logical :: is_hno3
 803       logical :: rls_flag(lpar)
 804       logical :: rnew_flag(lpar)
 805       logical :: cf_trigger(lpar)
 806       logical :: freezing(lpar)
 807
 808       character(len=132) :: message
 809
 810
 811       real, parameter  :: TFROZ      = 240.
 812       real, parameter  :: CFMIN      = 1.0
 813       real, parameter  :: CWMIN      = 1.0e-9
 814       real, parameter  :: DMIN       = 1.0e-1    !mm
 815       real, parameter  :: VOLPOW     = 1./3.
 816       real, parameter  :: RHORAIN    = 1.0e3     !kg/m3
 817       real, parameter  :: RHOSNOWFIX = 1.0e2     !kg/m3
 818       real, parameter  :: COLEFFRAIN = 0.7
 819       real, parameter  :: COLEFFAER  = 0.05
 820
 821 !-----------------------------------------------------------------------
 822 !  Note: LE must be less than LPAR
 823 !-----------------------------------------------------------------------
 824       LE = LPAR - 1
 825       rls_flag(1:le) = rls(1:le) > zero
 826       freezing(1:le) = tem(1:le) < tice
 827       rlsog(1:le) = rls(1:le)/garea
 828
 829 species_loop : &
 830      do N = 1,NTRACE
 831        QTT(:lpar)    = QTTJFL(:lpar,N)
 832        QTTNEW(:lpar) = QTTJFL(:lpar,N)
 833        is_hno3 = n == hno3_ndx
 834        if( is_hno3 ) then
 835          rca_wrk(:lpar) = zero
 836          fca_wrk(:lpar) = zero
 837          rcxa_wrk(:lpar) = zero
 838          fcxa_wrk(:lpar) = zero
 839          rcxb_wrk(:lpar) = zero
 840          fcxb_wrk(:lpar) = zero
 841          rls_wrk(:lpar) = zero
 842          rnew_wrk(:lpar) = zero
 843          cf_trigger(:lpar) = .false.
 844          clwx_wrk(:lpar) = -9999.
 845          deltarime_wrk(:lpar) = -9999.
 846          rax_wrk(:lpar,:) = zero
 847          fax_wrk(:lpar,:) = zero
 848        endif
 849 !-----------------------------------------------------------------------
 850 !  calculate scavenging by large-scale stratiform precipitation
 851 !  check whether mass-limited or henry's law
 852 !-----------------------------------------------------------------------
 853        if( TCKAQB(N) ) then
 854          LWASHTYP = 1
 855        else
 856          LWASHTYP = 2
 857        end if
 858 !-----------------------------------------------------------------------
 859 !  check whether soluble in ice
 860 !-----------------------------------------------------------------------
 861        if( TCNION(N) ) then
 862          LICETYP = 1
 863        else
 864          LICETYP = 2
 865        end if
 866
 867 !-----------------------------------------------------------------------
 868 !  initialization
 869 !-----------------------------------------------------------------------
 870        QTTOPAA = zero
 871        QTTOPCA = zero
 872
 873        RCA  = zero
 874        FCA  = zero
 875        DCA  = zero
 876        RAMA = zero
 877        FAMA = zero
 878        DAMA = zero
 879
 880        AMPCT      = zero
 881        AMCLPCT    = zero
 882        CLNEWPCT   = zero
 883        CLNEWAMPCT = zero
 884        CLOLDPCT   = zero
 885        CLOLDAMPCT = zero
 886 !-----------------------------------------------------------------------
 887 !  Check whether precip in top layer - if so, require CF ge 0.2
 888 !-----------------------------------------------------------------------
 889        if( RLS(LE) > zero ) then
 890          CFXX(LE) = max( CFMIN,CFR(LE) )
 891        else
 892          CFXX(LE) = CFR(LE)
 893        endif
 894
 895        rnew_flag(1:le) = .false.
 896
 897 level_loop : &
 898        do L = LE,1,-1
 899          LM1  = L - 1
 900          FAX  = zero
 901          RAX  = zero
 902          DAX  = zero
 903          FCXA = zero
 904          FCXB = zero
 905          DCXA = zero
 906          DCXB = zero
 907          RCXA = zero
 908          RCXB = zero
 909
 910          QTDISCF   = zero
 911          QTDISRIME = zero
 912          QTDISCXA  = zero
 913
 914          QTEVAPAXP = zero
 915          QTEVAPAXW = zero
 916          QTEVAPAX  = zero
 917          QTWASHAX  = zero
 918
 919          QTEVAPCXAP = zero
 920          QTEVAPCXAW = zero
 921          QTEVAPCXA  = zero
 922          QTRIMECXA  = zero
 923          QTWASHCXA  = zero
 924          QTRAINCXA  = zero
 925          QTRAINCXB  = zero
 926
 927          RAMPCT = zero
 928          RCXPCT = zero
 929
 930          RCXALOC = zero
 931          RCXBLOC = zero
 932          RAXLOC  = zero
 933          RAMALOC = zero
 934          RCALOC  = zero
 935
 936          RPRECIP       = zero
 937          DELTARIMEMASS = zero
 938          DELTARIME     = zero
 939          DOR           = zero
 940          DNEW          = zero
 941
 942          QTTOPAAX = zero
 943          QTTOPCAX = zero
 944
 945 has_rls : &
 946          if( rls_flag(l) ) then
 947 !-----------------------------------------------------------------------
 948 !-----Evaporate ambient precip and decrease area-------------------------
 949 !-----If ice, diam=diam falling from above  If rain, diam=4mm (not used)
 950 !-----Evaporate tracer contained in evaporated precip
 951 !-----Can't evaporate more than we start with-----------------------------
 952 !-----Don't do washout until we adjust ambient precip to match Rbot if needed
 953 !------(after RNEW if statements)
 954 !-----------------------------------------------------------------------
 955            FAX = max( zero,FAMA*(one - evaprate(l)) )
 956            RAX = RAMA                                                                !kg/m2/s
 957            if( FAMA > zero ) then
 958              if( freezing(l) ) then
 959                DAX = DAMA      !mm
 960              else
 961                DAX = four    !mm - not necessary
 962              endif
 963            else
 964              DAX = zero
 965            endif
 966
 967            if( RAMA > zero ) then
 968              QTEVAPAXP = min( QTTOPAA,EVAPRATE(L)*QTTOPAA )
 969            else
 970              QTEVAPAXP = zero
 971            endif
 972            if( is_hno3 ) then
 973              rax_wrk(l,1) = rax
 974              fax_wrk(l,1) = fax
 975            endif
 976
 977
 978 !-----------------------------------------------------------------------
 979 !  Determine how much the in-cloud precip rate has increased------
 980 !-----------------------------------------------------------------------
 981            WRK = RAX*FAX + RCA*FCA
 982            if( WRK > 0. ) then
 983              RNEW_TST = RLS(L)/(GAREA * WRK)
 984            else
 985              RNEW_TST = 10.
 986            endif
 987            RNEW = RLSOG(L) - (RAX*FAX + RCA*FCA)     !GBA*CF
 988            rnew_wrk(l) = rnew_tst
 989 !-----------------------------------------------------------------------
 990 !  if RNEW>0, there is growth and/or new precip formation
 991 !-----------------------------------------------------------------------
 992 has_rnew:  if( rlsog(l) > adj_factor*(rax*fax + rca*fca) ) then
 993 !-----------------------------------------------------------------------
 994 !  Min cloudwater requirement for cloud with new precip
 995 !  Min CF is set at top for LE, at end for other levels
 996 !  CWMIN is only needed for new precip formation - do not need for RNEW<0
 997 !-----------------------------------------------------------------------
 998              if( cfxx(l) == zero ) then
 999                write(*,*) 'offline inputs'
1000                write(*,*) qttjfl(:,n)
1001                write(*,*) qm(:)
1002                write(*,*) pofl(:)
1003                write(*,*) delz(:)
1004                write(*,*) rls(:)
1005                write(*,*) clwc(:)
1006                write(*,*) ciwc(:)
1007                write(*,*) cfr(:)
1008                write(*,*) tem(:)
1009                write(*,*) evaprate(:)
1010                write(*,*) hstar(:,n)
1011                write(message,'('' washout: cloud fraction == 0 @ i,j,l,n = '',4i4)') ii,jj,l,n
1012                call wrf_debug( 15, trim(message) )
1013                QTTJFL(:lpar,N) = QTT(:lpar)
1014                cycle species_loop
1015              endif
1016              rnew_flag(l) = .true.
1017              CLWX = max( CLWC(L)+CIWC(L),CWMIN*CFXX(L) )
1018              if( is_hno3 ) then
1019                clwx_wrk(l) = clwx
1020              endif
1021 !-----------------------------------------------------------------------
1022 !  Area of old cloud and new cloud
1023 !-----------------------------------------------------------------------
1024              FCXA = FCA
1025              FCXB = max( zero,CFXX(L)-FCXA )
1026 !-----------------------------------------------------------------------
1027 !                           ICE
1028 !  For ice and mixed phase, grow precip in old cloud by riming
1029 !  Use only portion of cloudwater in old cloud fraction
1030 !  and rain above old cloud fraction
1031 !  COLEFF from Lohmann and Roeckner (1996), Loss rate from Rotstayn (1997)
1032 !-----------------------------------------------------------------------
1033 is_freezing : &
1034              if( freezing(l) ) then
1035                COLEFFSNOW = exp( 2.5e-2*(TEM(L) - TICE) )
1036                if( TEM(L) <= TFROZ ) then
1037                  RHOSNOW = RHOSNOWFIX
1038                else
1039                  RHOSNOW = 0.303*(TEM(L) - TFROZ)*RHOSNOWFIX
1040                endif
1041                if( FCXA > zero ) then
1042                  if( DCA > zero ) then
1043                    DELTARIMEMASS = CLWX*QM(L)*(FCXA/CFXX(L))* &
1044                      (one - exp( (-COLEFFSNOW/(DCA*1.e-3))*((RCA)/(2.*RHOSNOW))*DTSCAV ))   !uses GBA R
1045                  else
1046                    DELTARIMEMASS = zero
1047                  endif
1048                else
1049                  DELTARIMEMASS = zero
1050                endif
1051 !-----------------------------------------------------------------------
1052 !  Increase in precip rate due to riming (kg/m2/s):
1053 !  Limit to total increase in R in cloud
1054 !-----------------------------------------------------------------------
1055                if( FCXA > zero ) then
1056                  DELTARIME = min( RNEW/FCXA,DELTARIMEMASS/(FCXA*GAREA*DTSCAV) ) !GBA
1057                else
1058                  DELTARIME = zero
1059                endif
1060                if( is_hno3 ) then
1061                  deltarime_wrk(l) = deltarime
1062                endif
1063 !-----------------------------------------------------------------------
1064 !  Find diameter of rimed precip, must be at least .1mm
1065 !-----------------------------------------------------------------------
1066                if( RCA > zero ) then
1067                  DOR = max( DMIN,(((RCA+DELTARIME)/RCA)**VOLPOW)*DCA )
1068                else
1069                  DOR = zero
1070                endif
1071 !-----------------------------------------------------------------------
1072 !  If there is some in-cloud precip left, we have new precip formation
1073 !  Will be spread over whole cloud fraction
1074 !-----------------------------------------------------------------------
1075 !  Calculate precip rate in old and new cloud fractions
1076 !-----------------------------------------------------------------------
1077                RPRECIP = (RNEW-(DELTARIME*FCXA))/CFXX(L) !kg/m2/s    !GBA
1078 !-----------------------------------------------------------------------
1079 !  Calculate precip rate in old and new cloud fractions
1080 !-----------------------------------------------------------------------
1081                RCXA = RCA + DELTARIME + RPRECIP          !kg/m2/s GBA
1082                RCXB = RPRECIP                            !kg/m2/s GBA
1083
1084 !-----------------------------------------------------------------------
1085 !  Find diameter of new precip from empirical relation using Rprecip
1086 !  in given area of box- use density of water, not snow, to convert kg/s
1087 !  to mm/s -> as given in Field and Heymsfield
1088 !  Also calculate diameter of mixed precip,DCXA, from empirical relation
1089 !  using total R in FCXA - this will give larger particles than averaging DOR and
1090 !  DNEW in the next level
1091 !  DNEW and DCXA must be at least .1mm
1092 !-----------------------------------------------------------------------
1093                if( RPRECIP > zero ) then
1094                  WEMP = (CLWX*QM(L))/(GAREA*CFXX(L)*DELZ(L)) !kg/m3
1095                  REMP = RPRECIP/((RHORAIN/1.e3))             !mm/s local
1096                  DNEW = DEMPIRICAL( WEMP, REMP )
1097                  DNEW = max( DMIN,DNEW )
1098                  if( FCXB > zero ) then
1099                    DCXB = DNEW
1100                  else
1101                    DCXB = zero
1102                  endif
1103                else
1104                  DCXB = zero
1105                endif
1106
1107                if( FCXA > zero ) then
1108                  WEMP = (CLWX*QM(L)*(FCXA/CFXX(L)))/(GAREA*FCXA*DELZ(L)) !kg/m3
1109                  REMP = RCXA/((RHORAIN/1.e3))                         !mm/s local
1110                  DEMP = DEMPIRICAL( WEMP, REMP )
1111                  DCXA = ((RCA+DELTARIME)/RCXA)*DOR + (RPRECIP/RCXA)*DNEW
1112                  DCXA = max( DEMP,DCXA )
1113                  DCXA = max( DMIN,DCXA )
1114                else
1115                  DCXA = zero
1116                endif
1117
1118                if( QTT(L) > zero ) then
1119 !-----------------------------------------------------------------------
1120 !                       ICE SCAVENGING
1121 !-----------------------------------------------------------------------
1122 !  For ice, rainout only hno3/aerosols using new precip
1123 !  Tracer dissolved given by Kaercher and Voigt (2006) for T<258K
1124 !  For T>258K, use Henry's Law with Retention coefficient
1125 !  Rain out in whole CF
1126 !-----------------------------------------------------------------------
1127                  if( RPRECIP > zero ) then
1128                    if( LICETYP == 1 ) then
1129                      RRAIN = RPRECIP*GAREA                                  !kg/s local
1130                      call DISGAS( CLWX, CFXX(L), TCMASS(N), HSTAR(L,N), &
1131                                   TEM(L),POFL(L),QM(L),                 &
1132 !++mmb
1133 !                                 QTT(L)*CFXX(L),QTDISCF )
1134                                   QTT(L)*CFXX(L),QTDISCF, is_hno3, RETEFF(N) )
1135 !--mmb
1136                      call RAINGAS( RRAIN, DTSCAV, CLWX, CFXX(L),        &
1137                                    QM(L), QTT(L), QTDISCF, QTRAIN )
1138                      WRK       = QTRAIN/CFXX(L)
1139                      QTRAINCXA = FCXA*WRK
1140                      QTRAINCXB = FCXB*WRK
1141                    elseif( LICETYP == 2 ) then
1142                      QTRAINCXA = zero
1143                      QTRAINCXB = zero
1144                    endif
1145                  endif
1146 !-----------------------------------------------------------------------
1147 !  For ice, accretion removal for hno3 and aerosols is propotional to riming,
1148 !  no accretion removal for gases
1149 !  remove only in mixed portion of cloud
1150 !  Limit DELTARIMEMASS to RNEW*DTSCAV for ice - evaporation of rimed ice to match
1151 !  RNEW precip rate would result in HNO3 escaping from ice (no trapping)
1152 !-----------------------------------------------------------------------
1153                  if( DELTARIME > zero ) then
1154                    if( LICETYP == 1 ) then
1155                      if( TEM(L) <= TFROZ ) then
1156                        RHOSNOW = RHOSNOWFIX
1157                      else
1158                        RHOSNOW = 0.303*(TEM(L) - TFROZ)*RHOSNOWFIX
1159                      endif
1160                      QTCXA = QTT(L)*FCXA
1161                      call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N),   &
1162                                   HSTAR(L,N), TEM(L), POFL(L),            &
1163 !++mmb
1164 !                                 QM(L), QTCXA, QTDISRIME )
1165                                   QM(L), QTCXA, QTDISRIME, is_hno3, RETEFF(N) )
1166 !--mmb
1167                      QTDISSTAR = (QTDISRIME*QTCXA)/(QTDISRIME + QTCXA)
1168                      QTRIMECXA = QTCXA*                             &
1169                         (one - exp((-COLEFFSNOW/(DCA*1.e-3))*       &
1170                         (RCA/(2.*RHOSNOW))*                         &  !uses GBA R
1171                         (QTDISSTAR/QTCXA)*DTSCAV))
1172                      QTRIMECXA = min( QTRIMECXA, &
1173                         ((RNEW*GAREA*DTSCAV)/(CLWX*QM(L)*(FCXA/CFXX(L))))*QTDISSTAR)
1174                    elseif( LICETYP == 2 ) then
1175                      QTRIMECXA = zero
1176                    endif
1177                  endif
1178                else
1179                  QTRAINCXA = zero
1180                  QTRAINCXB = zero
1181                  QTRIMECXA = zero
1182                endif
1183 !-----------------------------------------------------------------------
1184 !  For ice, no washout in interstitial cloud air
1185 !-----------------------------------------------------------------------
1186                QTWASHCXA = zero
1187                QTEVAPCXA = zero
1188
1189 !-----------------------------------------------------------------------
1190 !                      RAIN
1191 !  For rain, accretion increases rain rate but diameter remains constant
1192 !  Diameter is 4mm (not used)
1193 !-----------------------------------------------------------------------
1194              else is_freezing
1195                if( FCXA > zero ) then
1196                  DELTARIMEMASS = (CLWX*QM(L))*(FCXA/CFXX(L))*           &
1197                    (one - exp( -0.24*COLEFFRAIN*((RCA)**0.75)*DTSCAV ))  !local
1198                else
1199                  DELTARIMEMASS = zero
1200                endif
1201 !-----------------------------------------------------------------------
1202 !  Increase in precip rate due to riming (kg/m2/s):
1203 !  Limit to total increase in R in cloud
1204 !-----------------------------------------------------------------------
1205                if( FCXA > zero ) then
1206                  DELTARIME = min( RNEW/FCXA,DELTARIMEMASS/(FCXA*GAREA*DTSCAV) ) !GBA
1207                else
1208                  DELTARIME = zero
1209                endif
1210 !-----------------------------------------------------------------------
1211 !  If there is some in-cloud precip left, we have new precip formation
1212 !-----------------------------------------------------------------------
1213                RPRECIP = (RNEW-(DELTARIME*FCXA))/CFXX(L)       !GBA
1214
1215                RCXA = RCA + DELTARIME + RPRECIP            !kg/m2/s GBA
1216                RCXB = RPRECIP                              !kg/m2/s GBA
1217                DCXA = FOUR
1218                if( FCXB > zero ) then
1219                  DCXB = FOUR
1220                else
1221                  DCXB = zero
1222                endif
1223 !-----------------------------------------------------------------------
1224 !                         RAIN SCAVENGING
1225 !  For rain, rainout both hno3/aerosols and gases using new precip
1226 !-----------------------------------------------------------------------
1227                if( QTT(L) > zero ) then
1228                  if( RPRECIP > zero ) then
1229                    RRAIN = (RPRECIP*GAREA) !kg/s local
1230                    call DISGAS( CLWX, CFXX(L), TCMASS(N), HSTAR(L,N), &
1231                                 TEM(L), POFL(L), QM(L),               &
1232 !++mmb
1233 !                               QTT(L)*CFXX(L), QTDISCF )
1234                                 QTT(L)*CFXX(L), QTDISCF, is_hno3, RETEFF(N) )
1235 !--mmb
1236                    call RAINGAS( RRAIN, DTSCAV, CLWX, CFXX(L),        &
1237                                  QM(L), QTT(L), QTDISCF, QTRAIN )
1238                    WRK       = QTRAIN/CFXX(L)
1239                    QTRAINCXA = FCXA*WRK
1240                    QTRAINCXB = FCXB*WRK
1241                  endif
1242 !-----------------------------------------------------------------------
1243 !  For rain, accretion removal is propotional to riming
1244 !  caclulate for hno3/aerosols and gases
1245 !  Remove only in mixed portion of cloud
1246 !  Limit DELTARIMEMASS to RNEW*DTSCAV
1247 !-----------------------------------------------------------------------
1248                  if( DELTARIME > zero ) then
1249                    QTCXA = QTT(L)*FCXA
1250                    call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N),    &
1251                                 HSTAR(L,N), TEM(L), POFL(L),             &
1252 !++mmb
1253 !                               QM(L), QTCXA, QTDISRIME )
1254                                 QM(L), QTCXA, QTDISRIME, is_hno3, RETEFF(N) )
1255 !--mmb
1256                    QTDISSTAR = (QTDISRIME*QTCXA)/(QTDISRIME + QTCXA)
1257                    QTRIMECXA = QTCXA*                              &
1258                       (one - exp(-0.24*COLEFFRAIN*                 &
1259                       ((RCA)**0.75)*                               & !local
1260                       (QTDISSTAR/QTCXA)*DTSCAV))
1261                    QTRIMECXA = min( QTRIMECXA, &
1262                       ((RNEW*GAREA*DTSCAV)/(CLWX*QM(L)*(FCXA/CFXX(L))))*QTDISSTAR)
1263                  else
1264                    QTRIMECXA = zero
1265                  endif
1266                else
1267                  QTRAINCXA = zero
1268                  QTRAINCXB = zero
1269                  QTRIMECXA = zero
1270                endif
1271 !-----------------------------------------------------------------------
1272 !  For rain, washout gases and HNO3/aerosols using rain from above old cloud
1273 !  Washout for HNO3/aerosols is only on non-dissolved portion, impaction-style
1274 !  Washout for gases is on non-dissolved portion, limited by QTTOP+QTRIME
1275 !-----------------------------------------------------------------------
1276                if( RCA > zero ) then
1277                  QTPRECIP = FCXA*QTT(L) - QTDISRIME
1278                  if( LWASHTYP == 1 ) then
1279                    if( QTPRECIP > zero ) then
1280                      QTWASHCXA = QTPRECIP*(one - exp( -0.24*COLEFFAER*((RCA)**0.75)*DTSCAV ))   !local
1281                    else
1282                      QTWASHCXA = zero
1283                    endif
1284                    QTEVAPCXA = zero
1285                  elseif( LWASHTYP == 2 ) then
1286                    RWASH = RCA*GAREA                                !kg/s local
1287                    if( QTPRECIP > zero ) then
1288                      call WASHGAS( RWASH, FCA, DTSCAV, QTTOPCA+QTRIMECXA, &
1289                                    HSTAR(L,N), TEM(L), POFL(L),           &
1290                                    QM(L), QTPRECIP, QTWASHCXA, QTEVAPCXA )
1291                    else
1292                      QTWASHCXA = zero
1293                      QTEVAPCXA = zero
1294                    endif
1295                  endif
1296                endif
1297              endif is_freezing
1298 !-----------------------------------------------------------------------
1299 !  If RNEW<O, confine precip to area of cloud above
1300 !  FCXA does not require a minimum (could be zero if R(L).le.what
1301 !  evaporated in ambient)
1302 !-----------------------------------------------------------------------
1303            else has_rnew
1304              CLWX = CLWC(L) + CIWC(L)
1305              if( is_hno3 ) then
1306                clwx_wrk(l) = clwx
1307              endif
1308              FCXA = FCA
1309              FCXB = max( zero,CFXX(L)-FCXA )
1310              RCXB = zero
1311              DCXB = zero
1312              QTRAINCXA = zero
1313              QTRAINCXB = zero
1314              QTRIMECXA = zero
1315
1316 !-----------------------------------------------------------------------
1317 !  Put rain into cloud up to RCA so that we evaporate
1318 !  from ambient first
1319 !  Adjust ambient to try to match RLS(L)
1320 !  If no cloud, RAX=R(L)
1321 !-----------------------------------------------------------------------
1322              if( FCXA > zero ) then
1323                RCXA = min( RCA,RLS(L)/(GAREA*FCXA) )     !kg/m2/s  GBA
1324                if( FAX > zero .and. ((RCXA+1.e-12) < RLS(L)/(GAREA*FCXA)) ) then
1325                  RAXADJF = RLS(L)/GAREA - RCXA*FCXA
1326                  RAMPCT = RAXADJF/(RAX*FAX)
1327                  FAXADJ = RAMPCT*FAX
1328                  if( FAXADJ > zero ) then
1329                    RAXADJ = RAXADJF/FAXADJ
1330                  else
1331                    RAXADJ = zero
1332                  endif
1333                else
1334                  RAXADJ = zero
1335                  RAMPCT = zero
1336                  FAXADJ = zero
1337                endif
1338              else
1339                RCXA = zero
1340                if( FAX > zero ) then
1341                  RAXADJF = RLS(L)/GAREA
1342                  RAMPCT = RAXADJF/(RAX*FAX)
1343                  FAXADJ = RAMPCT*FAX
1344                  if( FAXADJ > zero ) then
1345                    RAXADJ = RAXADJF/FAXADJ
1346                  else
1347                    RAXADJ = zero
1348                  endif
1349                else
1350                  RAXADJ = zero
1351                  RAMPCT = zero
1352                  FAXADJ = zero
1353                endif
1354              endif
1355
1356              QTEVAPAXP = min( QTTOPAA,QTTOPAA - (RAMPCT*(QTTOPAA-QTEVAPAXP)) )
1357              FAX = FAXADJ
1358              RAX = RAXADJ
1359
1360 !-----------------------------------------------------------------------
1361 !                IN-CLOUD EVAPORATION/WASHOUT
1362 !  If precip out the bottom of the cloud is 0, evaporate everything
1363 !  If there is no cloud, QTTOPCA=0, so nothing happens
1364 !-----------------------------------------------------------------------
1365              if( RCXA <= zero ) then
1366                QTEVAPCXA = QTTOPCA
1367                RCXA = zero
1368                DCXA = zero
1369              else
1370 !-----------------------------------------------------------------------
1371 !  If rain out the bottom of the cloud is >0 (but .le. RCA):
1372 !  For ice, decrease particle size,
1373 !  no washout
1374 !  no evap for non-ice gases (b/c there is nothing in ice)
1375 !  T<Tmix,release hno3& aerosols
1376 !  release is amount dissolved in ice mass released
1377 !  T>Tmix, hno3&aerosols are incorporated into ice structure:
1378 !  do not release
1379 !  For rain, assume full evaporation of some raindrops
1380 !  proportional evaporation for all species
1381 !  washout for gases using Rbot
1382 !  impact washout for hno3/aerosol portion in gas phase
1383 !-----------------------------------------------------------------------
1384 !              if (TEM(L) < TICE ) then
1385 is_freezing_a : &
1386                if( freezing(l) ) then
1387                  QTWASHCXA = zero
1388                  DCXA = ((RCXA/RCA)**VOLPOW)*DCA
1389                  if( LICETYP == 1 ) then
1390                    if( TEM(L) <= TMIX ) then
1391                      MASSLOSS = (RCA-RCXA)*FCXA*GAREA*DTSCAV
1392 !-----------------------------------------------------------------------
1393 !  note-QTT doesn't matter b/c T<258K
1394 !-----------------------------------------------------------------------
1395                      call DISGAS( (MASSLOSS/QM(L)), FCXA, TCMASS(N),   &
1396                                    HSTAR(L,N), TEM(L), POFL(L),        &
1397 !++mmb
1398 !                                  QM(L), QTT(L), QTEVAPCXA )
1399                                    QM(L), QTT(L), QTEVAPCXA, is_hno3, RETEFF(N) )
1400 !--mmb
1401                      QTEVAPCXA = min( QTTOPCA,QTEVAPCXA )
1402                    else
1403                      QTEVAPCXA = zero
1404                    endif
1405                  elseif( LICETYP == 2 ) then
1406                    QTEVAPCXA = zero
1407                  endif
1408                else is_freezing_a
1409                  QTEVAPCXAP = (RCA - RCXA)/RCA*QTTOPCA
1410                  DCXA = FOUR
1411                  QTCXA = FCXA*QTT(L)
1412                  if( LWASHTYP == 1 ) then
1413                    if( QTT(L) > zero ) then
1414                      call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N),   &
1415                                   HSTAR(L,N), TEM(L), POFL(L),            &
1416 !++mmb
1417 !                                 QM(L), QTCXA, QTDISCXA )
1418                                   QM(L), QTCXA, QTDISCXA, is_hno3, RETEFF(N) )
1419 !--mmb
1420                      if( QTCXA > QTDISCXA ) then
1421                        QTWASHCXA = (QTCXA - QTDISCXA)*(one - exp( -0.24*COLEFFAER*((RCXA)**0.75)*DTSCAV )) !local
1422                      else
1423                        QTWASHCXA = zero
1424                      endif
1425                      QTEVAPCXAW = zero
1426                    else
1427                      QTWASHCXA  = zero
1428                      QTEVAPCXAW = zero
1429                    endif
1430                  elseif (LWASHTYP == 2 ) then
1431                    RWASH = RCXA*GAREA                         !kg/s local
1432                    call WASHGAS( RWASH, FCXA, DTSCAV, QTTOPCA, HSTAR(L,N), &
1433                                  TEM(L), POFL(L), QM(L),                   &
1434                                  QTCXA-QTDISCXA, QTWASHCXA, QTEVAPCXAW )
1435                  endif
1436                  QTEVAPCXA = QTEVAPCXAP + QTEVAPCXAW
1437                endif is_freezing_a
1438              endif
1439            endif has_rnew
1440
1441 !-----------------------------------------------------------------------
1442 !                 AMBIENT WASHOUT
1443 !  Ambient precip is finalized - if it is rain, washout
1444 !  no ambient washout for ice, since gases are in vapor phase
1445 !-----------------------------------------------------------------------
1446            if( RAX > zero ) then
1447 !            if( TEM(L) >= TICE ) then
1448              if( .not. freezing(l) ) then
1449                QTAX = FAX*QTT(L)
1450                if( LWASHTYP == 1 ) then
1451                  QTWASHAX = QTAX*                        &
1452                     (one - exp(-0.24*COLEFFAER*       &
1453                    ((RAX)**0.75)*DTSCAV))  !local
1454                  QTEVAPAXW = zero
1455                elseif( LWASHTYP == 2 ) then
1456                  RWASH = RAX*GAREA   !kg/s local
1457                  call WASHGAS( RWASH, FAX, DTSCAV, QTTOPAA, HSTAR(L,N), &
1458                                TEM(L), POFL(L), QM(L), QTAX,            &
1459                                QTWASHAX, QTEVAPAXW )
1460                endif
1461              else
1462                QTEVAPAXW = zero
1463                QTWASHAX  = zero
1464              endif
1465            else
1466              QTEVAPAXW = zero
1467              QTWASHAX  = zero
1468            endif
1469            QTEVAPAX = QTEVAPAXP + QTEVAPAXW
1470
1471 !-----------------------------------------------------------------------
1472 !                  END SCAVENGING
1473 !  Require CF if our ambient evaporation rate would give less
1474 !  precip than R from model.
1475 !-----------------------------------------------------------------------
1476            if( is_hno3 ) then
1477              rls_wrk(l) = rls(l)/garea
1478              rca_wrk(l) = rca
1479              fca_wrk(l) = fca
1480              rcxa_wrk(l) = rcxa
1481              fcxa_wrk(l) = fcxa
1482              rcxb_wrk(l) = rcxb
1483              fcxb_wrk(l) = fcxb
1484              rax_wrk(l,2) = rax
1485              fax_wrk(l,2) = fax
1486            endif
1487 upper_level : &
1488            if( L > 1 ) then
1489              FAMA = max( FCXA + FCXB + FAX - CFR(LM1),zero )
1490              if( FAX > zero ) then
1491                RAXLOC = RAX/FAX
1492              else
1493                RAXLOC = zero
1494              endif
1495              if( FCXA > zero ) then
1496                RCXALOC = RCXA/FCXA
1497              else
1498                RCXALOC = zero
1499              endif
1500              if( FCXB > zero ) then
1501                RCXBLOC = RCXB/FCXB
1502              else
1503                RCXBLOC = zero
1504              endif
1505
1506              if( CFR(LM1) >= CFMIN ) then
1507                CFXX(LM1) = CFR(LM1)
1508              else
1509                if( adj_factor*RLSOG(LM1) >= (RCXA*FCXA + RCXB*FCXB + RAX*FAX)*(one - EVAPRATE(LM1)) ) then
1510                  CFXX(LM1) = CFMIN
1511                  cf_trigger(lm1) = .true.
1512                else
1513                  CFXX(LM1) = CFR(LM1)
1514                endif
1515              endif
1516 !-----------------------------------------------------------------------
1517 !  Figure out what will go into ambient and cloud below
1518 !  Don't do for lowest level
1519 !-----------------------------------------------------------------------
1520              if( FAX > zero ) then
1521                RAXLOC = RAX/FAX
1522                AMPCT = max( zero,min( one,(CFXX(L) + FAX - CFXX(LM1))/FAX ) )
1523                AMCLPCT = one - AMPCT
1524              else
1525                RAXLOC  = zero
1526                AMPCT   = zero
1527                AMCLPCT = zero
1528              endif
1529              if( FCXB > zero ) then
1530                RCXBLOC = RCXB/FCXB
1531                CLNEWPCT = max( zero,min( (CFXX(LM1) - FCXA)/FCXB,one ) )
1532                CLNEWAMPCT = one - CLNEWPCT
1533              else
1534                RCXBLOC    = zero
1535                CLNEWPCT   = zero
1536                CLNEWAMPCT = zero
1537              endif
1538              if( FCXA > zero ) then
1539                RCXALOC = RCXA/FCXA
1540                CLOLDPCT = max( zero,min( CFXX(LM1)/FCXA,one ) )
1541                CLOLDAMPCT = one - CLOLDPCT
1542              else
1543                RCXALOC    = zero
1544                CLOLDPCT   = zero
1545                CLOLDAMPCT = zero
1546              endif
1547 !-----------------------------------------------------------------------
1548 !  Remix everything for the next level
1549 !-----------------------------------------------------------------------
1550              FCA = min( CFXX(LM1),FCXA*CLOLDPCT + CLNEWPCT*FCXB + AMCLPCT*FAX )
1551              if( FCA > zero ) then
1552 !-----------------------------------------------------------------------
1553 !  Maintain cloud core by reducing NC and AM area going into cloud below
1554 !-----------------------------------------------------------------------
1555                RCA = (RCXA*FCXA*CLOLDPCT + RCXB*FCXB*CLNEWPCT + RAX*FAX*AMCLPCT)/FCA
1556                if( RCA > zero ) then
1557                  DCA = (RCXA*FCXA*CLOLDPCT)/(RCA*FCA)*DCXA + &
1558                        (RCXB*FCXB*CLNEWPCT)/(RCA*FCA)*DCXB + &
1559                        (RAX*FAX*AMCLPCT)/(RCA*FCA)*DAX
1560                else
1561                  DCA = zero
1562                endif
1563              else
1564                FCA = zero
1565                DCA = zero
1566                RCA = zero
1567              endif
1568
1569              FAMA = FCXA + FCXB + FAX - CFXX(LM1)
1570              if( FAMA > zero ) then
1571                RAMA = (RCXA*FCXA*CLOLDAMPCT + RCXB*FCXB*CLNEWAMPCT + RAX*FAX*AMPCT)/FAMA
1572                if( RAMA > zero ) then
1573                  DAMA = (RCXA*FCXA*CLOLDAMPCT)/(RAMA*FAMA)*DCXA +  &
1574                         (RCXB*FCXB*CLNEWAMPCT)/(RAMA*FAMA)*DCXB +  &
1575                         (RAX*FAX*AMPCT)/(RAMA*FAMA)*DAX
1576                else
1577                   FAMA = zero
1578                   DAMA = zero
1579                endif
1580              else
1581                FAMA = zero
1582                DAMA = zero
1583                RAMA = zero
1584              endif
1585            else upper_level
1586              AMPCT      = zero
1587              AMCLPCT    = zero
1588              CLNEWPCT   = zero
1589              CLNEWAMPCT = zero
1590              CLOLDPCT   = zero
1591              CLOLDAMPCT = zero
1592            endif upper_level
1593          else has_rls
1594            RNEW = zero
1595            QTEVAPCXA = QTTOPCA
1596            QTEVAPAX = QTTOPAA
1597            if( L > 1 ) then
1598              if( RLS(LM1) > zero ) then
1599                CFXX(LM1) = max( CFMIN,CFR(LM1) )
1600 !              if( CFR(LM1) >= CFMIN ) then
1601 !                CFXX(LM1) = CFR(LM1)
1602 !              else
1603 !                CFXX(LM1) = CFMIN
1604 !              endif
1605              else
1606                CFXX(LM1) = CFR(LM1)
1607              endif
1608            endif
1609            AMPCT      = zero
1610            AMCLPCT    = zero
1611            CLNEWPCT   = zero
1612            CLNEWAMPCT = zero
1613            CLOLDPCT   = zero
1614            CLOLDAMPCT = zero
1615            RCA        = zero
1616            RAMA       = zero
1617            FCA        = zero
1618            FAMA       = zero
1619            DCA        = zero
1620            DAMA       = zero
1621          endif has_rls
1622
1623          if( is_hno3 ) then
1624            fama_wrk(l) = fama
1625            rama_wrk(l) = rama
1626          endif
1627 !-----------------------------------------------------------------------
1628 !  Net loss can not exceed QTT in each region
1629 !-----------------------------------------------------------------------
1630          QTNETLCXA = QTRAINCXA + QTRIMECXA + QTWASHCXA - QTEVAPCXA
1631          QTNETLCXA = min( QTT(L)*FCXA,QTNETLCXA )
1632
1633          QTNETLCXB =QTRAINCXB
1634          QTNETLCXB = min( QTT(L)*FCXB,QTNETLCXB )
1635
1636          QTNETLAX = QTWASHAX - QTEVAPAX
1637          QTNETLAX = min( QTT(L)*FAX,QTNETLAX )
1638
1639          QTTNEW(L) = QTT(L) - (QTNETLCXA + QTNETLCXB + QTNETLAX)
1640
1641
1642          QTTOPCAX = (QTTOPCA + QTNETLCXA)*CLOLDPCT + QTNETLCXB*CLNEWPCT + (QTTOPAA + QTNETLAX)*AMCLPCT
1643          QTTOPAAX = (QTTOPCA + QTNETLCXA)*CLOLDAMPCT + QTNETLCXB*CLNEWAMPCT + (QTTOPAA + QTNETLAX)*AMPCT
1644          QTTOPCA = QTTOPCAX
1645          QTTOPAA = QTTOPAAX
1646        end do level_loop
1647
1648 !-----------------------------------------------------------------------
1649 !  reload new tracer mass and rescale moments: check upper limits (LE)
1650 !-----------------------------------------------------------------------
1651        QTTJFL(:le,N) = QTTNEW(:le)
1652      end do species_loop
1653
1654 end subroutine WASHOUT
1655
1656 subroutine DISGAS( CLWX, CFX, MOLMASS, HSTAR, &
1657                    TM, PR, QM, &
1658 !++mmb
1659 !                  QT, QTDIS )
1660                    QT, QTDIS, is_hno3, RETEFF )
1661 !--mmb
1662
1663 !-----------------------------------------------------------------------
1664 !  dummy arguments
1665 !-----------------------------------------------------------------------
1666       real, intent(in) :: CLWX,CFX    !cloud water,cloud fraction
1667       real, intent(in) :: MOLMASS     !molecular mass of tracer
1668       real, intent(in) :: HSTAR       !Henry's Law coeffs A*exp(-B/T)
1669       real, intent(in) :: TM          !temperature of box (K)
1670       real, intent(in) :: PR          !pressure of box (hPa)
1671       real, intent(in) :: QM          !air mass in box (kg)
1672       real, intent(in) :: QT          !tracer in box (kg)
1673       real, intent(out) :: QTDIS      !tracer dissolved in aqueous phase
1674
1675 !++mmb
1676       logical, intent(in) :: is_hno3
1677       real, intent(in) :: RETEFF      !Ice retention parameter
1678 !--mmb
1679
1680 !-----------------------------------------------------------------------
1681 !  local variables
1682 !-----------------------------------------------------------------------
1683
1684 !++mmb
1685 !     real, parameter  :: RETEFF = 0.5
1686 !     real, parameter  :: RETEFF = 1.0
1687 !--mmb
1688
1689       real  :: MUEMP
1690
1691 !-----------------------------------------------------------------------
1692 !  calculate rate of uptake of tracer
1693 !-----------------------------------------------------------------------
1694
1695 !-----------------------------------------------------------------------
1696 !  effective Henry's Law constant: H* = moles-T / liter-precip / press(atm-T)
1697 !  p(atm of tracer-T) = (QT/QM) * (.029/MolWt-T) * pressr(hPa)/1000
1698 !  limit temperature effects to T above freezing
1699 !  MU from fit to Kaercher and Voigt (2006)
1700 !-----------------------------------------------------------------------
1701       if( TM >= TICE ) then
1702         QTDIS = (HSTAR*(QT/(QM*CFX))*0.029*(PR/1.0e3))*(CLWX*QM)
1703 !++mmb
1704 !     elseif( TM <= TMIX ) then
1705       elseif( TM <= TMIX .and. is_hno3 ) then
1706 !--mmb
1707         MUEMP = exp( -14.2252 + TM*(1.55704e-1 - 7.1929e-4*TM) )
1708         QTDIS = MUEMP*(MOLMASS/18.)*(CLWX*QM)
1709       else
1710         QTDIS = RETEFF*((HSTAR*(QT/(QM*CFX))*0.029*(PR/1.0e3))*(CLWX*QM))
1711       endif
1712
1713 end subroutine DISGAS
1714
1715 subroutine RAINGAS( RRAIN, DTSCAV, CLWX, CFX, QM, QT, QTDIS, QTRAIN )
1716 !-----------------------------------------------------------------------
1717 !---New trace-gas rainout from large-scale precip with two time scales,
1718 !---one based on precip formation from cloud water and one based on
1719 !---Henry's Law solubility: correct limit for delta-t
1720 !---
1721 !---NB this code does not consider the aqueous dissociation (eg, C-q)
1722 !---   that makes uptake of HNO3 and H2SO4 so complete.  To do so would
1723 !---   require that we keep track of the pH of the falling rain.
1724 !---THUS the Henry's Law coefficient KHA needs to be enhanced to incldue this!
1725 !---ALSO the possible formation of other soluble species from, eg, CH2O, H2O2
1726 !---   can be considered with enhanced values of KHA.
1727 !---
1728 !---Does NOT now use RMC (moist conv rain) but could, assuming 30% coverage
1729 !-----------------------------------------------------------------------
1730
1731 !-----------------------------------------------------------------------
1732 !  dummy arguments
1733 !-----------------------------------------------------------------------
1734       real, intent(in) :: RRAIN       !new rain formation in box (kg/s)
1735       real, intent(in) :: DTSCAV      !time step (s)
1736       real, intent(in) :: CLWX,CFX    !cloud water and cloud fraction
1737       real, intent(in) :: QM          !air mass in box (kg)
1738       real, intent(in) :: QT          !tracer in box (kg)
1739       real, intent(in) :: QTDIS       !tracer in aqueous phase (kg)
1740       real, intent(out) :: QTRAIN     !tracer picked up by new rain
1741
1742 !-----------------------------------------------------------------------
1743 !  local variables
1744 !-----------------------------------------------------------------------
1745       real ::  QTLF, QTDISSTAR
1746
1747       QTDISSTAR = (QTDIS*(QT*CFX))/(QTDIS+(QT*CFX))
1748
1749 !-----------------------------------------------------------------------
1750 !  Tracer Loss frequency (1/s) within cloud fraction:
1751 !-----------------------------------------------------------------------
1752       QTLF = (RRAIN*QTDISSTAR)/(CLWX*QM*QT*CFX)
1753
1754 !-----------------------------------------------------------------------
1755 !  in time = DTSCAV, the amount of QTT scavenged is calculated
1756 !  from CF*AMOUNT OF UPTAKE
1757 !-----------------------------------------------------------------------
1758       QTRAIN = QT*CFX*(one - exp( -DTSCAV*QTLF ))
1759
1760 end subroutine RAINGAS
1761
1762 subroutine WASHGAS( RWASH, BOXF, DTSCAV, QTRTOP, HSTAR, TM, PR, QM, &
1763                     QT, QTWASH, QTEVAP )
1764 !-----------------------------------------------------------------------
1765 !---for most gases below-cloud washout assume Henry-Law equilib with precip
1766 !---assumes that precip is liquid, if frozen, do not call this sub
1767 !---since solubility is moderate, fraction of box with rain does not matter
1768 !---NB this code does not consider the aqueous dissociation (eg, C-q)
1769 !---   that makes uptake of HNO3 and H2SO4 so complete.  To do so would
1770 !---   require that we keep track of the pH of the falling rain.
1771 !---THUS the Henry's Law coefficient KHA needs to be enhanced to incldue this!
1772 !---ALSO the possible formation of other soluble species from, eg, CH2O, H2O2
1773 !---   can be considered with enhanced values of KHA.
1774 !-----------------------------------------------------------------------
1775
1776 !-----------------------------------------------------------------------
1777 !  dummy arguments
1778 !-----------------------------------------------------------------------
1779       real, intent(in)  :: RWASH   ! precip leaving bottom of box (kg/s)
1780       real, intent(in)  :: BOXF    ! fraction of box with washout
1781       real, intent(in)  :: DTSCAV  ! time step (s)
1782       real, intent(in)  :: QTRTOP  ! tracer-T in rain entering top of box over time step (kg)
1783       real, intent(in)  :: HSTAR   ! Henry's Law coeffs A*exp(-B/T)
1784       real, intent(in)  :: TM      ! temperature of box (K)
1785       real, intent(in)  :: PR      ! pressure of box (hPa)
1786       real, intent(in)  :: QT      ! tracer in box (kg)
1787       real, intent(in)  :: QM      ! air mass in box (kg)
1788       real, intent(out) :: QTWASH  ! tracer picked up by precip (kg)
1789       real, intent(out) :: QTEVAP  ! tracer evaporated from precip (kg)
1790
1791 !-----------------------------------------------------------------------
1792 !  local variables
1793 !-----------------------------------------------------------------------
1794       real            :: FWASH, QTMAX, QTDIF
1795
1796 !-----------------------------------------------------------------------
1797 !  effective Henry's Law constant: H* = moles-T / liter-precip / press(atm-T)
1798 !  p(atm of tracer-T) = (QT/QM) * (.029/MolWt-T) * pressr(hPa)/1000
1799 !  limit temperature effects to T above freezing
1800 !-----------------------------------------------------------------------
1801
1802 !-----------------------------------------------------------------------
1803 !  effective washout frequency (1/s):
1804 !-----------------------------------------------------------------------
1805       FWASH = (RWASH*HSTAR*29.e-6*PR)/(QM*BOXF)
1806 !-----------------------------------------------------------------------
1807 !  equilib amount of T (kg) in rain thru bottom of box over time step
1808 !-----------------------------------------------------------------------
1809       QTMAX = QT*FWASH*DTSCAV
1810       if( QTMAX > QTRTOP ) then
1811 !-----------------------------------------------------------------------
1812 !  more of tracer T can go into rain
1813 !-----------------------------------------------------------------------
1814          QTDIF = min( QT,QTMAX-QTRTOP )
1815          QTWASH = QTDIF * (one - exp( -DTSCAV*FWASH ))
1816          QTEVAP = zero
1817       else
1818 !-----------------------------------------------------------------------
1819 !  too much of T in rain, must degas/evap T
1820 !-----------------------------------------------------------------------
1821          QTWASH = zero
1822          QTEVAP = QTRTOP - QTMAX
1823       endif
1824
1825 end subroutine WASHGAS
1826
1827 function DEMPIRICAL( CWATER, RRATE )
1828
1829 !-----------------------------------------------------------------------
1830 !  dummy arguments
1831 !-----------------------------------------------------------------------
1832       real, intent(in)  :: CWATER
1833       real, intent(in)  :: RRATE
1834
1835 !-----------------------------------------------------------------------
1836 !  local variables
1837 !-----------------------------------------------------------------------
1838       real(8), parameter :: big_diameter = 100._8            ! mm
1839       real(8), parameter :: const0       = .638_8
1840       real(8), parameter :: const1       = oner8 + const0
1841
1842       real(8) :: RRATEX, WX, THETA, PHI, ETA, BETA, ALPHA, BEE
1843       real(8) :: GAMTHETA, GAMBETA
1844       real(8) :: numer, denom
1845       real(8) :: diameter                      ! mm
1846
1847       real :: DEMPIRICAL
1848
1849       if( cwater > 0._8 ) then
1850         RRATEX = real(RRATE,kind=8)*3600._8       !mm/hr
1851         WX     = real(CWATER,kind=8)*1.0e3_8      !g/m3
1852
1853         if( RRATEX > 0.04_8 ) then
1854            THETA = exp( -1.43_8*log10( 7._8*RRATEX ) ) + 2.8_8
1855         else
1856            THETA = 5._8
1857         endif
1858
1859         PHI = RRATEX/(3600._8*10._8) !cgs units
1860         ETA = exp( 3.01_8*THETA - 10.5_8 )
1861         BETA  = THETA/const1
1862         ALPHA = exp( FOURR8*(BETA - 3.5_8) )
1863         BEE   = const0*BETA - ONER8
1864         GAMTHETA = real( GAMMA( real(THETA,kind=4) ),kind=8 )
1865         GAMBETA  = real( GAMMA( real(BETA + ONER8,kind=4) ),kind=8 )
1866
1867         numer = WX*ETA*GAMTHETA
1868         denom = 1.0e6_8*ALPHA*PHI*GAMBETA
1869         diameter = ((numer/denom)**(-oner8/BEE))*10._8
1870         diameter = min( big_diameter,diameter )
1871         DEMPIRICAL = real( diameter )
1872 !     DEMPIRICAL = (((WX*ETA*GAMTHETA)/(1.0e6*ALPHA*PHI*GAMBETA))** (-one/BEE))*10.  ! in mm (wx/1e6 for cgs)
1873       else
1874         DEMPIRICAL = 0.
1875       endif
1876
1877 end function DEMPIRICAL
1878
1879 function GAMMA( X )
1880 !-----------------------------------------------------------------------
1881 !       Purpose: Compute the gamma function Ahat(x)
1882 !       Input :  x  --- Argument of Ahat(x)
1883 !                       ( x is not equal to 0,-1,-2, ... )
1884 !       Output:  GA --- Ahat(x)
1885 !-----------------------------------------------------------------------
1886
1887 !-----------------------------------------------------------------------
1888 !  dummy arguments
1889 !-----------------------------------------------------------------------
1890         real, intent(in)  :: X
1891
1892 !-----------------------------------------------------------------------
1893 !  local variables
1894 !-----------------------------------------------------------------------
1895         real, parameter :: PI = 3.141592653589793e0
1896
1897         integer :: k, M, M1
1898         real    :: GR, R, Z
1899         real    :: G(26)
1900
1901 !-----------------------------------------------------------------------
1902 !  function definition
1903 !-----------------------------------------------------------------------
1904         real :: GAMMA
1905
1906         DATA G/1.0e0,0.5772156649015329,                     &
1907             -0.6558780715202538e0, -0.420026350340952e-1,    &
1908             0.1665386113822915e0,-.421977345555443e-1,       &
1909             -.96219715278770e-2, .72189432466630e-2,         &
1910             -.11651675918591e-2, -.2152416741149e-3,         &
1911             .1280502823882e-3, -.201348547807e-4,            &
1912             -.12504934821e-5, .11330272320e-5,               &
1913             -.2056338417e-6, .61160950e-8,                   &
1914             .50020075e-8, -.11812746e-8,                     &
1915             .1043427e-9, .77823e-11,                         &
1916             -.36968e-11, .51e-12,                            &
1917             -.206e-13, -.54e-14, .14e-14, .1e-15/
1918
1919 is_integer : &
1920         IF( x == real( int(x) ) ) then
1921           IF( X > zero ) THEN
1922             GAMMA = ONE
1923             M1 = INT(X) - 1
1924             DO K = 2,M1
1925               GAMMA = GAMMA*real(K)
1926             END DO
1927           ELSE
1928             GAMMA = 1.0e36
1929           ENDIF
1930         ELSE is_integer
1931           IF( ABS(X) > ONE ) THEN
1932             Z = ABS(X)
1933             M = INT(Z)
1934             R = ONE
1935             DO K = 1,M
1936               R = R*(Z - real(k))
1937             END DO
1938             Z = Z - real(M)
1939           ELSE
1940             Z = X
1941           ENDIF
1942           GR = G(26)
1943           DO K = 25,1,-1
1944             GR = GR*Z + G(K)
1945           end DO
1946           GAMMA = ONE/(GR*Z)
1947           IF( ABS(X) > ONE ) THEN
1948             GAMMA = GAMMA*R
1949             IF( X < zero ) then
1950               GAMMA = -PI/(X*GAMMA*SIN( PI*X ))
1951             ENDIF
1952           ENDIF
1953         ENDIF is_integer
1954
1955 END function GAMMA
1956
1957 END MODULE module_mozcart_wetscav