LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_neu_wetdep.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 534 | 793 | 67.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 7 | 7 | 100.0 %

              Line data    Source code
           1             : !
           2             : ! code written by J.-F. Lamarque, S. Walters and F. Vitt
           3             : ! based on the original code from J. Neu developed for UC Irvine
           4             : ! model
           5             : !
           6             : ! LKE 2/23/2018 - correct setting flag for mass-limited (HNO3,etc.) vs Henry's Law washout
           7             : !
           8             : module mo_neu_wetdep
           9             : !
          10             :   use shr_kind_mod,     only : r8 => shr_kind_r8
          11             :   use cam_logfile,      only : iulog
          12             :   use constituents,     only : pcnst
          13             :   use spmd_utils,       only : masterproc
          14             :   use cam_abortutils,   only : endrun
          15             :   use shr_drydep_mod,   only : n_species_table, species_name_table, dheff
          16             :   use gas_wetdep_opts,  only : gas_wetdep_method, gas_wetdep_list, gas_wetdep_cnt
          17             : !
          18             :   implicit none
          19             : !
          20             :   private
          21             :   public :: neu_wetdep_init
          22             :   public :: neu_wetdep_tend
          23             : !
          24             :   save
          25             : !
          26             :   integer, allocatable, dimension(:) :: mapping_to_heff,mapping_to_mmr
          27             :   real(r8),allocatable, dimension(:) :: mol_weight
          28             :   logical ,allocatable, dimension(:) :: ice_uptake
          29             :   integer                     :: index_cldice,index_cldliq,nh3_ndx,co2_ndx,so2_ndx
          30             :   integer                     :: so4_ndx,so4s_ndx ! geos-chem
          31             :   logical                     :: debug   = .false.
          32             :   integer                     :: hno3_ndx = 0
          33             : !
          34             : ! diagnostics
          35             : !
          36             :   logical                     :: do_diag = .false.
          37             :   integer, parameter          :: kdiag = 18
          38             : !
          39             :   real(r8), parameter :: zero = 0._r8
          40             :   real(r8), parameter :: one  = 1._r8
          41             : !
          42             :   logical :: do_neu_wetdep
          43             : !
          44             :   real(r8), parameter  :: TICE=263._r8
          45             :
          46             : contains
          47             :
          48             : !-----------------------------------------------------------------------
          49             : !-----------------------------------------------------------------------
          50             : !
          51        1536 : subroutine neu_wetdep_init
          52             : !
          53             :   use constituents, only : cnst_get_ind,cnst_mw
          54             :   use cam_history,  only : addfld, add_default, horiz_only
          55             :   use phys_control, only : phys_getopts, cam_chempkg_is
          56             : !
          57             :   integer :: m,l
          58             :   character*20 :: test_name
          59             :
          60             :   logical :: history_chemistry
          61             :
          62        1536 :   call phys_getopts(history_chemistry_out=history_chemistry)
          63             :
          64        1536 :   do_neu_wetdep = gas_wetdep_method == 'NEU' .and. gas_wetdep_cnt>0
          65             :
          66        1536 :   if (.not.do_neu_wetdep) return
          67             :
          68        4608 :   allocate( mapping_to_heff(gas_wetdep_cnt) )
          69        3072 :   allocate( mapping_to_mmr(gas_wetdep_cnt) )
          70        3072 :   allocate( ice_uptake(gas_wetdep_cnt) )
          71        4608 :   allocate( mol_weight(gas_wetdep_cnt) )
          72             :
          73             : !
          74             : ! find mapping to heff table
          75             : !
          76        1536 :   if ( debug ) then
          77           0 :     print '(a,i4)','gas_wetdep_cnt=',gas_wetdep_cnt
          78           0 :     print '(a,i4)','n_species_table=',n_species_table
          79             :   end if
          80        9216 :   mapping_to_heff = -99
          81        9216 :   do m=1,gas_wetdep_cnt
          82             : !
          83        7680 :     test_name = gas_wetdep_list(m)
          84        7680 :     if ( debug ) print '(i4,a)',m,trim(test_name)
          85             : !
          86             : ! mapping based on the MOZART4 wet removal subroutine;
          87             : ! this might need to be redone (JFL: Sep 2010)
          88             : !
          89             : ! Skip mapping if using GEOS-Chem; all GEOS-Chem species are in dep_data_file
          90             : ! (heff table) specified in namelist drv_flds_in (EWL: Dec 2022)
          91        7680 :   if ( .not. cam_chempkg_is('geoschem_mam4') ) then
          92       15360 :     select case( trim(test_name) )
          93             : !
          94             : ! CCMI: added SO2t and NH_50W
          95             : !
          96             :       case ( 'SOGB','SOGI','SOGM','SOGT','SOGX' )
          97           0 :          test_name = 'H2O2'
          98             :       case ( 'SO2t' )
          99           0 :          test_name = 'SO2'
         100             :       case ( 'CLONO2','BRONO2','HCL','HOCL','HOBR','HBR', 'Pb', 'HF', 'COF2', 'COFCL')
         101           0 :          test_name = 'HNO3'
         102             :       case ( 'NH_50W', 'NDEP', 'NHDEP', 'NH4', 'NH4NO3' )
         103           0 :          test_name = 'HNO3'
         104             :       case(  'SOAGbb0' )  ! Henry's Law coeff. added for VBS SOA's, biomass burning is the same as fossil fuels
         105           0 :          test_name = 'SOAGff0'
         106             :       case(  'SOAGbb1' )
         107           0 :          test_name = 'SOAGff1'
         108             :       case(  'SOAGbb2' )
         109           0 :          test_name = 'SOAGff2'
         110             :       case(  'SOAGbb3' )
         111           0 :          test_name = 'SOAGff3'
         112             :       case(  'SOAGbb4' )
         113       15360 :          test_name = 'SOAGff4'
         114             :     end select
         115             :   endif
         116             : !
         117      789504 :     do l = 1,n_species_table
         118             : !
         119             : !      if ( debug ) print '(i4,a)',l,trim(species_name_table(l))
         120             : !
         121      789504 :        if( trim(test_name) == trim( species_name_table(l) ) ) then
         122        7680 :           mapping_to_heff(m)  = l
         123        7680 :           if ( debug ) print '(a,a,i4)','mapping to heff of ',trim(species_name_table(l)),l
         124        7680 :           exit
         125             :        end if
         126             :     end do
         127        7680 :     if ( mapping_to_heff(m) == -99 ) then
         128           0 :       if (masterproc) print *,'problem with mapping_to_heff of ',trim(test_name)
         129             : !      call endrun()
         130             :     end if
         131             : !
         132             : ! special cases for NH3 and CO2
         133             : !
         134        7680 :     if ( trim(test_name) == 'NH3' ) then
         135           0 :       nh3_ndx = m
         136             :     end if
         137        7680 :     if ( trim(test_name) == 'CO2' ) then
         138           0 :       co2_ndx = m
         139             :     end if
         140        7680 :     if ( trim(gas_wetdep_list(m)) == 'HNO3' ) then
         141           0 :       hno3_ndx = m
         142             :     end if
         143        7680 :     if ( trim(test_name) == 'SO2' ) then
         144        1536 :       so2_ndx = m
         145             :     end if
         146        7680 :     if ( trim(test_name) == 'SO4' ) then ! GEOS-Chem bulk sulfate
         147           0 :       so4_ndx = m
         148             :     end if
         149        9216 :     if ( trim(test_name) == 'SO4S' ) then ! GEOS-Chem bulk sulfate on surface seasalt
         150           0 :       so4s_ndx = m
         151             :     end if
         152             : !
         153             :   end do
         154             :
         155        9216 :    if (any ( mapping_to_heff(:) == -99 ))  call endrun('mo_neu_wet->depwetdep_init: unmapped species error' )
         156             : !
         157        1536 :   if ( debug .and. masterproc ) then
         158           0 :     write(iulog, '(a,i4)') 'co2_ndx',co2_ndx
         159           0 :     write(iulog, '(a,i4)') 'nh3_ndx',nh3_ndx
         160           0 :     write(iulog, '(a,i4)') 'so2_ndx',so2_ndx
         161             :   end if
         162             : !
         163             : ! find mapping to species
         164             : !
         165        9216 :   mapping_to_mmr = -99
         166        9216 :   do m=1,gas_wetdep_cnt
         167        7680 :     if ( debug .and. masterproc ) write(iulog, '(i4,a)') m,trim(gas_wetdep_list(m))
         168        7680 :     call cnst_get_ind(gas_wetdep_list(m), mapping_to_mmr(m), abort=.false. )
         169        7680 :     if ( debug .and. masterproc ) write(iulog, '(a,i4)') 'mapping_to_mmr ',mapping_to_mmr(m)
         170        9216 :     if ( mapping_to_mmr(m) <= 0 ) then
         171           0 :       if (masterproc) write(iulog,*) 'problem with mapping_to_mmr of ',gas_wetdep_list(m)
         172           0 :       call endrun('neu_wetdep_init: problem with mapping_to_mmr of '//trim(gas_wetdep_list(m)))
         173             :     end if
         174             :   end do
         175             : !
         176             : ! define species-dependent arrays
         177             : !
         178        9216 :   do m=1,gas_wetdep_cnt
         179             : !
         180        7680 :     mol_weight(m) = cnst_mw(mapping_to_mmr(m))
         181        7680 :     if ( debug .and. masterproc ) write(iulog, '(i4,a,f8.4)') m,' mol_weight ',mol_weight(m)
         182        7680 :     ice_uptake(m) = .false.
         183        9216 :     if ( trim(gas_wetdep_list(m)) == 'HNO3' ) then
         184           0 :       ice_uptake(m) = .true.
         185             :     end if
         186             : !
         187             : !
         188             :   end do
         189             : !
         190             : ! indices for cloud quantities
         191             : !
         192        1536 :   call cnst_get_ind( 'CLDICE', index_cldice )
         193        1536 :   call cnst_get_ind( 'CLDLIQ', index_cldliq )
         194             : !
         195             : ! define output
         196             : !
         197        9216 :   do m=1,gas_wetdep_cnt
         198       15360 :     call addfld     ('DTWR_'//trim(gas_wetdep_list(m)),(/ 'lev' /), 'A','kg/kg/s','wet removal Neu scheme tendency')
         199        7680 :     call addfld     ('WD_'//trim(gas_wetdep_list(m)),horiz_only, 'A','kg/m2/s','vertical integrated wet deposition flux')
         200       15360 :     call addfld     ('HEFF_'//trim(gas_wetdep_list(m)),(/ 'lev' /), 'A','M/atm','Effective Henrys Law coeff.')
         201        9216 :     if (history_chemistry) then
         202        7680 :        call add_default('WD_'//trim(gas_wetdep_list(m)), 1, ' ')
         203             :     end if
         204             :   end do
         205             : !
         206        1536 :   if ( do_diag ) then
         207           0 :     call addfld     ('QT_RAIN_HNO3',(/ 'lev' /), 'A','mol/mol/s','wet removal Neu scheme rain tendency')
         208           0 :     call addfld     ('QT_RIME_HNO3',(/ 'lev' /), 'A','mol/mol/s','wet removal Neu scheme rain tendency')
         209           0 :     call addfld     ('QT_WASH_HNO3',(/ 'lev' /), 'A','mol/mol/s','wet removal Neu scheme rain tendency')
         210           0 :     call addfld     ('QT_EVAP_HNO3',(/ 'lev' /), 'A','mol/mol/s','wet removal Neu scheme rain tendency')
         211           0 :     if (history_chemistry) then
         212           0 :        call add_default('QT_RAIN_HNO3',1,' ')
         213           0 :        call add_default('QT_RIME_HNO3',1,' ')
         214           0 :        call add_default('QT_WASH_HNO3',1,' ')
         215           0 :        call add_default('QT_EVAP_HNO3',1,' ')
         216             :     end if
         217             :   end if
         218             : !
         219             :   return
         220             : !
         221        3072 : end subroutine neu_wetdep_init
         222             : !
         223     1489176 : subroutine neu_wetdep_tend(lchnk,ncol,mmr,pmid,pdel,zint,tfld,delt, &
         224     1489176 :      prain, nevapr, cld, cmfdqr, wd_tend, wd_tend_int)
         225             : !
         226        1536 :   use ppgrid,           only : pcols, pver
         227             :   use phys_grid,        only : get_area_all_p, get_rlat_all_p
         228             :   use shr_const_mod,    only : SHR_CONST_REARTH,SHR_CONST_G
         229             :   use cam_history,      only : outfld
         230             :   use shr_const_mod,    only : pi => shr_const_pi
         231             : !
         232             :   implicit none
         233             : !
         234             :   integer,        intent(in)    :: lchnk,ncol
         235             :   real(r8),       intent(in)    :: mmr(pcols,pver,pcnst)    ! mass mixing ratio (kg/kg)
         236             :   real(r8),       intent(in)    :: pmid(pcols,pver)         ! midpoint pressures (Pa)
         237             :   real(r8),       intent(in)    :: pdel(pcols,pver)         ! pressure delta about midpoints (Pa)
         238             :   real(r8),       intent(in)    :: zint(pcols,pver+1)       ! interface geopotential height above the surface (m)
         239             :   real(r8),       intent(in)    :: tfld(pcols,pver)         ! midpoint temperature (K)
         240             :   real(r8),       intent(in)    :: delt                     ! timestep (s)
         241             : !
         242             :   real(r8),       intent(in)    :: prain(ncol, pver)
         243             :   real(r8),       intent(in)    :: nevapr(ncol, pver)
         244             :   real(r8),       intent(in)    :: cld(ncol, pver)
         245             :   real(r8),       intent(in)    :: cmfdqr(ncol, pver)
         246             :   real(r8),       intent(inout) :: wd_tend(pcols,pver,pcnst)
         247             :   real(r8),       intent(inout) :: wd_tend_int(pcols,pcnst)
         248             : !
         249             : ! local arrays and variables
         250             : !
         251             :   integer :: i,k,l,kk,m
         252             :   real(r8), parameter                       :: rearth = SHR_CONST_REARTH    ! radius earth (m)
         253             :   real(r8), parameter                       :: gravit = SHR_CONST_G         ! m/s^2
         254     2978352 :   real(r8), dimension(ncol)                 :: area, wk_out
         255     2978352 :   real(r8), dimension(ncol,pver)            :: cldice,cldliq,cldfrc,totprec,totevap,delz,p
         256     2978352 :   real(r8), dimension(ncol,pver)            :: rls,evaprate,mass_in_layer,temp
         257     2978352 :   real(r8), dimension(ncol,pver,gas_wetdep_cnt) :: trc_mass,heff,dtwr
         258     2978352 :   real(r8), dimension(ncol,pver,gas_wetdep_cnt) :: wd_mmr
         259     2978352 :   logical , dimension(gas_wetdep_cnt)           :: tckaqb
         260     2978352 :   integer , dimension(ncol)                 :: test_flag
         261             : !
         262             : ! arrays for HNO3 diagnostics
         263             : !
         264     2978352 :   real(r8), dimension(ncol,pver)            :: qt_rain,qt_rime,qt_wash,qt_evap
         265             : !
         266             : ! for Henry's law calculations
         267             : !
         268             :   real(r8), parameter       :: t0     = 298._r8
         269             :   real(r8), parameter       :: ph     = 1.e-5_r8
         270             :   real(r8), parameter       :: ph_inv = 1._r8/ph
         271             :   real(r8)                  :: e298, dhr
         272     2978352 :   real(r8), dimension(ncol) :: dk1s,dk2s,wrk
         273             :   real(r8) :: lats(pcols)
         274             :
         275             :   real(r8), parameter :: rad2deg = 180._r8/pi
         276             :
         277             : !
         278             : ! from cam/src/physics/cam/stratiform.F90
         279             : !
         280             :
         281     1489176 :   if (.not.do_neu_wetdep) return
         282             : !
         283             : ! don't do anything if there are no species to be removed
         284             : !
         285     1489176 :   if ( gas_wetdep_cnt == 0 ) return
         286             : !
         287             : ! reset output variables
         288             : !
         289     1489176 :    wd_tend_int = 0._r8
         290             : !
         291             : ! get area (in radians square)
         292             : !
         293     1489176 :   call get_area_all_p(lchnk, ncol, area)
         294    24865776 :   area = area * rearth**2                     ! in m^2
         295             : !
         296             : ! reverse order along the vertical before calling
         297             : ! J. Neu's wet removal subroutine
         298             : !
         299   139982544 :   do k=1,pver
         300   138493368 :     kk = pver - k + 1
         301  2314006344 :     do i=1,ncol
         302             : !
         303  2174023800 :       mass_in_layer(i,k) = area(i) * pdel(i,kk)/gravit          ! kg
         304             : !
         305  2174023800 :       cldice (i,k) = mmr(i,kk,index_cldice)                     ! kg/kg
         306  2174023800 :       cldliq (i,k) = mmr(i,kk,index_cldliq)                     ! kg/kg
         307  2174023800 :       cldfrc (i,k) = cld(i,kk)                                  ! unitless
         308             : !
         309             :       totprec(i,k) = (prain(i,kk)+cmfdqr(i,kk)) &
         310  2174023800 :                                   * mass_in_layer(i,k)          ! kg/s
         311  2174023800 :       totevap(i,k) = nevapr(i,kk) * mass_in_layer(i,k)          ! kg/s
         312             : !
         313  2174023800 :       delz(i,k) = zint(i,kk) - zint(i,kk+1)                     ! in m
         314             : !
         315  2174023800 :       temp(i,k) = tfld(i,kk)
         316             : !
         317             : ! convert tracer mass to kg to kg/kg
         318             : !
         319 13044142800 :       trc_mass(i,k,:) = mmr(i,kk,mapping_to_mmr(:)) * mass_in_layer(i,k)
         320             : !
         321  2312517168 :       p   (i,k) = pmid(i,kk) * 0.01_r8          ! in hPa
         322             : !
         323             :     end do
         324             :   end do
         325             : !
         326             : ! define array for tendency calculation (on model grid)
         327             : !
         328 11571520896 :   dtwr(1:ncol,:,:) = mmr(1:ncol,:,mapping_to_mmr(:))
         329             : !
         330             : ! compute 1) integrated precipitation flux across the interfaces (rls)
         331             : !         2) evaporation rate
         332             : !
         333    24865776 :   rls      (:,pver) = 0._r8
         334    24865776 :   evaprate (:,pver) = 0._r8
         335   138493368 :   do k=pver-1,1,-1
         336  2287651392 :     rls     (:,k) = max(0._r8,totprec(:,k)-totevap(:,k)+rls(:,k+1))
         337             :     !evaprate(:,k) = min(1._r8,totevap(:,k)/(rls(:,k+1)+totprec(:,k)+1.e-36_r8))
         338  2289140568 :     evaprate(:,k) = min(1._r8,totevap(:,k)/(rls(:,k+1)+1.e-36_r8))
         339             :   end do
         340             : !
         341             : ! compute effective Henry's law coefficients
         342             : !
         343 11571520896 :   heff = 0._r8
         344   139982544 :   do k=1,pver
         345             : !
         346   138493368 :     kk = pver - k + 1
         347             : !
         348  2312517168 :     wrk(:) = (t0-tfld(1:ncol,kk))/(t0*tfld(1:ncol,kk))
         349             : !
         350   832449384 :     do m=1,gas_wetdep_cnt
         351             : !
         352   692466840 :       l    = mapping_to_heff(m)
         353   692466840 :       e298 = dheff(1,l)
         354   692466840 :       dhr  = dheff(2,l)
         355 11562585840 :       heff(:,k,m) = e298*exp( dhr*wrk(:) )
         356 11562585840 :       test_flag = -99
         357   692466840 :       if( dheff(3,l) /= 0._r8 .and. dheff(5,l) == 0._r8 ) then
         358   138493368 :         e298 = dheff(3,l)
         359   138493368 :         dhr  = dheff(4,l)
         360  2312517168 :         dk1s(:) = e298*exp( dhr*wrk(:) )
         361 11285599104 :         where( heff(:,k,m) /= 0._r8 )
         362   138493368 :           heff(:,k,m) = heff(:,k,m)*(1._r8 + dk1s(:)*ph_inv)
         363             :         elsewhere
         364             :           test_flag = 1
         365   138493368 :           heff(:,k,m) = dk1s(:)*ph_inv
         366             :         endwhere
         367             :       end if
         368             : !
         369 11562585840 :       if (k.eq.1 .and. maxval(test_flag) > 0 .and. debug .and. masterproc ) then
         370           0 :          write(iulog, '(a,i4)') 'heff for m=',m
         371             :       endif
         372             : !
         373   830960208 :       if( dheff(5,l) /= 0._r8 ) then
         374   138493368 :         if( nh3_ndx > 0 .or. co2_ndx > 0 .or. so2_ndx > 0 .or. so4_ndx > 0 .or. so4s_ndx > 0 ) then
         375   138493368 :           e298 = dheff(3,l)
         376   138493368 :           dhr  = dheff(4,l)
         377  2312517168 :           dk1s(:) = e298*exp( dhr*wrk(:) )
         378   138493368 :           e298 = dheff(5,l)
         379   138493368 :           dhr  = dheff(6,l)
         380  2312517168 :           dk2s(:) = e298*exp( dhr*wrk(:) )
         381   138493368 :           if( m == co2_ndx .or. m == so2_ndx .or. m == so4_ndx .or. m == so4s_ndx ) then
         382  2312517168 :              heff(:,k,m) = heff(:,k,m)*(1._r8 + dk1s(:)*ph_inv*(1._r8 + dk2s(:)*ph_inv))
         383           0 :           else if( m == nh3_ndx ) then
         384           0 :              heff(:,k,m) = heff(:,k,m)*(1._r8 + dk1s(:)*ph/dk2s(:))
         385             :           else
         386           0 :              if ( masterproc ) write(iulog,*) 'error in assigning henrys law coefficients'
         387           0 :              if ( masterproc ) write(iulog,*) 'species ',m
         388             :           end if
         389             :         end if
         390             :       end if
         391             : !
         392             :     end do
         393             :   end do
         394             : !
         395     1489176 :   if ( debug .and. masterproc ) then
         396           0 :     write(iulog,'(a,50f8.2)')  'tckaqb     ',tckaqb
         397           0 :     write(iulog,'(a,50e12.4)') 'heff       ',heff(1,1,:)
         398           0 :     write(iulog,'(a,50i4)')    'ice_uptake ',ice_uptake
         399           0 :     write(iulog,'(a,50f8.2)')  'mol_weight ',mol_weight(:)
         400           0 :     write(iulog,'(a,50f8.2)')  'temp       ',temp(1,:)
         401           0 :     write(iulog,'(a,50f8.2)')  'p          ',p   (1,:)
         402             :   end if
         403             : !
         404             : ! call J. Neu's subroutine
         405             : !
         406    24865776 :   do i=1,ncol
         407             : !
         408    93506400 :     call washo(pver,gas_wetdep_cnt,delt,trc_mass(i,:,:),mass_in_layer(i,:),p(i,:),delz(i,:) &
         409   140259600 :               ,rls(i,:),cldliq(i,:),cldice(i,:),cldfrc(i,:),temp(i,:),evaprate(i,:) &
         410    23376600 :               ,area(i),heff(i,:,:),mol_weight(:),tckaqb(:),ice_uptake(:) &
         411 70201418976 :               ,qt_rain(i,:),qt_rime(i,:),qt_wash(i,:),qt_evap(i,:) )
         412             : !
         413             :   end do
         414             : !
         415             : ! compute tendencies and convert back to mmr
         416             : ! on original vertical grid
         417             : !
         418   139982544 :   do k=1,pver
         419   138493368 :     kk = pver - k + 1
         420  2314006344 :     do i=1,ncol
         421             : !
         422             : ! convert tracer mass from kg
         423             : !
         424 13182636168 :       wd_mmr(i,kk,:) = trc_mass(i,k,:) / mass_in_layer(i,k)
         425             : !
         426             :     end do
         427             :   end do
         428             : !
         429             : ! tendency calculation (on model grid)
         430             : !
         431 11571520896 :   dtwr(1:ncol,:,:) = wd_mmr(1:ncol,:,:) - dtwr(1:ncol,:,:)
         432 11571520896 :   dtwr(1:ncol,:,:) = dtwr(1:ncol,:,:) / delt
         433             :
         434             : ! polarward of 60S, 60N and <200hPa set to zero!
         435     1489176 :   call get_rlat_all_p(lchnk, pcols, lats )
         436   139982544 :   do k = 1, pver
         437  2314006344 :     do i= 1, ncol
         438  2312517168 :       if ( abs( lats(i)*rad2deg ) > 60._r8 ) then
         439   265535088 :         if ( pmid(i,k) < 20000._r8) then
         440   967462686 :            dtwr(i,k,:) = 0._r8
         441             :         endif
         442             :       endif
         443             :     end do
         444             :   end do
         445             : !
         446             : ! output tendencies
         447             : !
         448     8935056 :   do m=1,gas_wetdep_cnt
         449 11570031720 :     wd_tend(1:ncol,:,mapping_to_mmr(m)) = wd_tend(1:ncol,:,mapping_to_mmr(m)) + dtwr(1:ncol,:,m)
         450     7445880 :     call outfld( 'DTWR_'//trim(gas_wetdep_list(m)),dtwr(:,:,m),ncol,lchnk )
         451             :
         452 11570031720 :     call outfld( 'HEFF_'//trim(gas_wetdep_list(m)),heff(:,pver:1:-1,m),ncol,lchnk )
         453             : !
         454             : ! vertical integrated wet deposition rate [kg/m2/s]
         455             : !
         456   124328880 :     wk_out = 0._r8
         457   699912720 :     do k=1,pver
         458   692466840 :       kk = pver - k + 1
         459 11570031720 :       wk_out(1:ncol) = wk_out(1:ncol) + (dtwr(1:ncol,k,m) * mass_in_layer(1:ncol,kk)/area(1:ncol))
         460             :     end do
         461     7445880 :     call outfld( 'WD_'//trim(gas_wetdep_list(m)),wk_out,ncol,lchnk )
         462             : !
         463             : ! to be used in mo_chm_diags to compute wet_deposition_NOy_as_N and wet_deposition_NHx_as_N (units: kg/m2/s)
         464             : !
         465     7445880 :     if ( debug .and. masterproc ) then
         466           0 :        write(iulog,*) 'mo_neu ',mapping_to_mmr(m),(wk_out(1:ncol))
         467             :     endif
         468   125818056 :     wd_tend_int(1:ncol,mapping_to_mmr(m)) = wk_out(1:ncol)
         469             : !
         470             :   end do
         471             : !
         472     1489176 :   if ( do_diag ) then
         473           0 :     call outfld('QT_RAIN_HNO3', qt_rain, ncol, lchnk )
         474           0 :     call outfld('QT_RIME_HNO3', qt_rime, ncol, lchnk )
         475           0 :     call outfld('QT_WASH_HNO3', qt_wash, ncol, lchnk )
         476           0 :     call outfld('QT_EVAP_HNO3', qt_evap, ncol, lchnk )
         477             :   end if
         478             : !
         479             :   return
         480     1489176 : end subroutine neu_wetdep_tend
         481             :
         482             : !-----------------------------------------------------------------------
         483             : !
         484             : ! Original code from Jessica Neu
         485             : ! Updated by S. Walters and J.-F. Lamarque (March-April 2011)
         486             : !
         487             : !-----------------------------------------------------------------------
         488             :
         489    23376600 :       subroutine WASHO(LPAR,NTRACE,DTSCAV,QTTJFL,QM,POFL,DELZ,  &
         490    23376600 :       RLS,CLWC,CIWC,CFR,TEM,EVAPRATE,GAREA,HSTAR,TCMASS,TCKAQB, &
         491    23376600 :       TCNION, qt_rain, qt_rime, qt_wash, qt_evap)
         492             : !
         493             :       implicit none
         494             :
         495             : !-----------------------------------------------------------------------
         496             : !---p-conde 5.4 (2007)   -----called from main-----
         497             : !---called from pmain to calculate rainout and washout of tracers
         498             : !---revised by JNEU 8/2007
         499             : !---
         500             : !-LAER has been removed - no scavenging for aerosols
         501             : !-LAER could be used as LWASHTYP
         502             : !---WILL THIS WORK FOR T42->T21???????????
         503             : !-----------------------------------------------------------------------
         504             :
         505             :       integer LPAR, NTRACE
         506             :       real(r8),  intent(inout) ::  QTTJFL(LPAR,NTRACE)
         507             :       real(r8),  intent(in) :: DTSCAV, QM(LPAR),POFL(LPAR),DELZ(LPAR),GAREA
         508             :       real(r8),  intent(in) :: RLS(LPAR),CLWC(LPAR),CIWC(LPAR),CFR(LPAR),TEM(LPAR),      &
         509             :                                EVAPRATE(LPAR)
         510             :       real(r8),  intent(in) :: HSTAR(LPAR,NTRACE),TCMASS(NTRACE)
         511             :       logical ,  intent(in) :: TCKAQB(NTRACE),TCNION(NTRACE)
         512             : !
         513             :       real(r8),  intent(inout) :: qt_rain(lpar)
         514             :       real(r8),  intent(inout) :: qt_rime(lpar)
         515             :       real(r8),  intent(inout) :: qt_wash(lpar)
         516             :       real(r8),  intent(inout) :: qt_evap(lpar)
         517             : !
         518             :       integer L,N,LE, LM1
         519    46753200 :       real(r8), dimension(LPAR) :: CFXX
         520    46753200 :       real(r8), dimension(LPAR) :: QTT, QTTNEW
         521             :
         522             :       real(r8) WRK, RNEW_TST
         523             :       real(r8) CLWX
         524             :       real(r8) RNEW,RPRECIP,DELTARIMEMASS,DELTARIME,RAMPCT
         525             :       real(r8) MASSLOSS
         526             :       real(r8) DOR,DNEW,DEMP,COLEFFSNOW,RHOSNOW
         527             :       real(r8) WEMP,REMP,RRAIN,RWASH
         528             :       real(r8) QTPRECIP,QTRAIN,QTCXA,QTAX
         529             :
         530             :       real(r8) FAMA,RAMA,DAMA,FCA,RCA,DCA
         531             :       real(r8) FAX,RAX,DAX,FCXA,RCXA,DCXA,FCXB,RCXB,DCXB
         532             :       real(r8) RAXADJ,FAXADJ,RAXADJF
         533             :       real(r8) QTDISCF,QTDISRIME,QTDISCXA
         534             :       real(r8) QTEVAPAXP,QTEVAPAXW,QTEVAPAX
         535             :       real(r8) QTWASHAX
         536             :       real(r8) QTEVAPCXAP,QTEVAPCXAW,QTEVAPCXA
         537             :       real(r8) QTWASHCXA,QTRIMECXA
         538             :       real(r8) QTRAINCXA,QTRAINCXB
         539             :       real(r8) QTTOPCA,QTTOPAA,QTTOPCAX,QTTOPAAX
         540             :
         541             :       real(r8) AMPCT,AMCLPCT,CLNEWPCT,CLNEWAMPCT,CLOLDPCT,CLOLDAMPCT
         542             :
         543             :       real(r8) QTNETLCXA,QTNETLCXB,QTNETLAX
         544             :       real(r8) QTDISSTAR
         545             :
         546             :
         547             :       real(r8), parameter  :: CFMIN=0.1_r8
         548             :       real(r8), parameter  :: CWMIN=1.0e-5_r8
         549             :       real(r8), parameter  :: DMIN=1.0e-1_r8       !mm
         550             :       real(r8), parameter  :: VOLPOW=1._r8/3._r8
         551             :       real(r8), parameter  :: RHORAIN=1.0e3_r8     !kg/m3
         552             :       real(r8), parameter  :: RHOSNOWFIX=1.0e2_r8     !kg/m3
         553             :       real(r8), parameter  :: COLEFFRAIN=0.7_r8
         554             :       real(r8), parameter  :: TMIX=258._r8
         555             :       real(r8), parameter  :: TFROZ=240._r8
         556             :       real(r8), parameter  :: COLEFFAER=0.05_r8
         557             : !
         558             : ! additional work arrays and diagnostics
         559             : !
         560    46753200 :       real(r8) :: rls_wrk(lpar)
         561    46753200 :       real(r8) :: rnew_wrk(lpar)
         562    46753200 :       real(r8) :: rca_wrk(lpar)
         563    46753200 :       real(r8) :: fca_wrk(lpar)
         564    46753200 :       real(r8) :: rcxa_wrk(lpar)
         565    46753200 :       real(r8) :: fcxa_wrk(lpar)
         566    46753200 :       real(r8) :: rcxb_wrk(lpar)
         567    46753200 :       real(r8) :: fcxb_wrk(lpar)
         568    46753200 :       real(r8) :: rax_wrk(lpar,2)
         569    46753200 :       real(r8) :: fax_wrk(lpar,2)
         570    46753200 :       real(r8) :: rama_wrk(lpar)
         571    46753200 :       real(r8) :: fama_wrk(lpar)
         572    46753200 :       real(r8) :: deltarime_wrk(lpar)
         573    46753200 :       real(r8) :: clwx_wrk(lpar)
         574    46753200 :       real(r8) :: frc(lpar,3)
         575    46753200 :       real(r8) :: rlsog(lpar)
         576             : !
         577             :       logical :: is_hno3
         578    46753200 :       logical :: rls_flag(lpar)
         579    46753200 :       logical :: rnew_flag(lpar)
         580    46753200 :       logical :: cf_trigger(lpar)
         581    46753200 :       logical :: freezing(lpar)
         582             : !
         583             :       real(r8), parameter :: four = 4._r8
         584             :       real(r8), parameter :: adj_factor = one + 10._r8*epsilon( one )
         585             : !
         586             :       integer :: LICETYP
         587             : !
         588    23376600 :       if ( debug .and. masterproc ) then
         589           0 :         write(iulog,'(a,50f8.2)')  'tckaqb     ',tckaqb
         590           0 :         write(iulog,'(a,50e12.4)') 'hstar      ',hstar(1,:)
         591           0 :         write(iulog,'(a,50i4)')    'ice_uptake ',TCNION
         592           0 :         write(iulog,'(a,50f8.2)')  'mol_weight ',TCMASS(:)
         593           0 :         write(iulog,'(a,50f8.2)')  'temp       ',tem(:)
         594           0 :         write(iulog,'(a,50f8.2)')  'p          ',pofl(:)
         595             :       end if
         596             :
         597             : !-----------------------------------------------------------------------
         598    23376600 :       LE = LPAR-1
         599             : !
         600  2174023800 :       rls_flag(1:le) = rls(1:le) > zero
         601  2174023800 :       freezing(1:le) = tem(1:le) < tice
         602  2174023800 :       rlsog(1:le) = rls(1:le)/garea
         603             : !
         604             : species_loop : &
         605   140259600 :      do N = 1,NTRACE
         606 10987002000 :        QTT(:lpar)    = QTTJFL(:lpar,N)
         607 10987002000 :        QTTNEW(:lpar) = QTTJFL(:lpar,N)
         608   116883000 :        is_hno3 = n == hno3_ndx
         609   116883000 :        if( is_hno3 ) then
         610           0 :          qt_rain(:lpar) = zero
         611           0 :          qt_rime(:lpar) = zero
         612           0 :          qt_wash(:lpar) = zero
         613           0 :          qt_evap(:lpar) = zero
         614           0 :          rca_wrk(:lpar) = zero
         615           0 :          fca_wrk(:lpar) = zero
         616           0 :          rcxa_wrk(:lpar) = zero
         617           0 :          fcxa_wrk(:lpar) = zero
         618           0 :          rcxb_wrk(:lpar) = zero
         619           0 :          fcxb_wrk(:lpar) = zero
         620           0 :          rls_wrk(:lpar) = zero
         621           0 :          rnew_wrk(:lpar) = zero
         622           0 :          cf_trigger(:lpar) = .false.
         623           0 :          clwx_wrk(:lpar) = -9999._r8
         624           0 :          deltarime_wrk(:lpar) = -9999._r8
         625           0 :          rax_wrk(:lpar,:) = zero
         626           0 :          fax_wrk(:lpar,:) = zero
         627             :        endif
         628             :
         629             : !-----------------------------------------------------------------------
         630             : !  check whether soluble in ice
         631             : !-----------------------------------------------------------------------
         632   116883000 :        if( TCNION(N) ) then
         633             :          LICETYP = 1
         634             :        else
         635   116883000 :          LICETYP = 2
         636             :        end if
         637             :
         638             : !-----------------------------------------------------------------------
         639             : !  initialization
         640             : !-----------------------------------------------------------------------
         641   116883000 :        QTTOPAA = zero
         642   116883000 :        QTTOPCA = zero
         643             :
         644   116883000 :        RCA  = zero
         645   116883000 :        FCA  = zero
         646   116883000 :        DCA  = zero
         647   116883000 :        RAMA = zero
         648   116883000 :        FAMA = zero
         649   116883000 :        DAMA = zero
         650             :
         651   116883000 :        AMPCT      = zero
         652   116883000 :        AMCLPCT    = zero
         653   116883000 :        CLNEWPCT   = zero
         654   116883000 :        CLNEWAMPCT = zero
         655   116883000 :        CLOLDPCT   = zero
         656   116883000 :        CLOLDAMPCT = zero
         657             : !-----------------------------------------------------------------------
         658             : !  Check whether precip in top layer - if so, require CF ge 0.2
         659             : !-----------------------------------------------------------------------
         660   116883000 :        if( RLS(LE) > zero ) then
         661           0 :          CFXX(LE) = max( CFMIN,CFR(LE) )
         662             :        else
         663   116883000 :          CFXX(LE) = CFR(LE)
         664             :        endif
         665             :
         666 10870119000 :        rnew_flag(1:le) = .false.
         667             :
         668             : level_loop : &
         669 10870119000 :        do L = LE,1,-1
         670 10753236000 :          LM1  = L - 1
         671 10753236000 :          FAX  = zero
         672 10753236000 :          RAX  = zero
         673 10753236000 :          DAX  = zero
         674 10753236000 :          FCXA = zero
         675 10753236000 :          FCXB = zero
         676 10753236000 :          DCXA = zero
         677 10753236000 :          DCXB = zero
         678 10753236000 :          RCXA = zero
         679 10753236000 :          RCXB = zero
         680             :
         681 10753236000 :          QTDISCF   = zero
         682 10753236000 :          QTDISRIME = zero
         683 10753236000 :          QTDISCXA  = zero
         684             :
         685 10753236000 :          QTEVAPAXP = zero
         686 10753236000 :          QTEVAPAXW = zero
         687 10753236000 :          QTEVAPAX  = zero
         688 10753236000 :          QTWASHAX  = zero
         689             :
         690 10753236000 :          QTEVAPCXAP = zero
         691 10753236000 :          QTEVAPCXAW = zero
         692 10753236000 :          QTEVAPCXA  = zero
         693 10753236000 :          QTRIMECXA  = zero
         694 10753236000 :          QTWASHCXA  = zero
         695 10753236000 :          QTRAINCXA  = zero
         696 10753236000 :          QTRAINCXB  = zero
         697             :
         698 10753236000 :          RAMPCT = zero
         699             :
         700 10753236000 :          RPRECIP       = zero
         701 10753236000 :          DELTARIMEMASS = zero
         702 10753236000 :          DELTARIME     = zero
         703 10753236000 :          DOR           = zero
         704 10753236000 :          DNEW          = zero
         705             :
         706 10753236000 :          QTTOPAAX = zero
         707 10753236000 :          QTTOPCAX = zero
         708             :
         709             : has_rls : &
         710 10753236000 :          if( rls_flag(l) ) then
         711             : !-----------------------------------------------------------------------
         712             : !-----Evaporate ambient precip and decrease area-------------------------
         713             : !-----If ice, diam=diam falling from above  If rain, diam=4mm (not used)
         714             : !-----Evaporate tracer contained in evaporated precip
         715             : !-----Can't evaporate more than we start with-----------------------------
         716             : !-----Don't do washout until we adjust ambient precip to match Rbot if needed
         717             : !------(after RNEW if statements)
         718             : !-----------------------------------------------------------------------
         719  2964352485 :            FAX = max( zero,FAMA*(one - evaprate(l)) )
         720  2964352485 :            RAX = RAMA                                                                !kg/m2/s
         721  2964352485 :            if ( debug ) then
         722           0 :              if( (l == 3 .or. l == 2) ) then
         723           0 :                write(*,*) 'washout: l,rls,fax = ',l,rls(l),fax
         724             :              endif
         725             :            endif
         726  2964352485 :            if( FAMA > zero ) then
         727  1463780260 :              if( freezing(l) ) then
         728   469266215 :                DAX = DAMA      !mm
         729             :              else
         730   994514045 :                DAX = four    !mm - not necessary
         731             :              endif
         732             :            else
         733  1500572225 :              DAX = zero
         734             :            endif
         735             :
         736  2964352485 :            if( RAMA > zero ) then
         737  1463780260 :              QTEVAPAXP = min( QTTOPAA,EVAPRATE(L)*QTTOPAA )
         738             :            else
         739             :              QTEVAPAXP = zero
         740             :            endif
         741  2964352485 :            if( is_hno3 ) then
         742           0 :              rax_wrk(l,1) = rax
         743           0 :              fax_wrk(l,1) = fax
         744             :            endif
         745             :
         746             :
         747             : !-----------------------------------------------------------------------
         748             : !  Determine how much the in-cloud precip rate has increased------
         749             : !-----------------------------------------------------------------------
         750  2964352485 :            WRK = RAX*FAX + RCA*FCA
         751  2964352485 :            if( WRK > 0._r8 ) then
         752  2795998695 :              RNEW_TST = RLS(L)/(GAREA * WRK)
         753             :            else
         754             :              RNEW_TST = 10._r8
         755             :            endif
         756  2964352485 :            RNEW = RLSOG(L) - (RAX*FAX + RCA*FCA)     !GBA*CF
         757  2964352485 :            rnew_wrk(l) = rnew_tst
         758  2964352485 :            if ( debug ) then
         759           0 :              if( is_hno3 .and. l == kdiag-1 ) then
         760           0 :                write(*,*) ' '
         761           0 :                write(*,*) 'washout: rls,rax,fax,rca,fca'
         762           0 :                write(*,'(1p,5g15.7)') rls(l),rax,fax,rca,fca
         763           0 :                write(*,*) ' '
         764             :              endif
         765             :            endif
         766             : !-----------------------------------------------------------------------
         767             : !  if RNEW>0, there is growth and/or new precip formation
         768             : !-----------------------------------------------------------------------
         769  2964352485 : has_rnew:  if( rlsog(l) > adj_factor*(rax*fax + rca*fca) ) then
         770             : !-----------------------------------------------------------------------
         771             : !  Min cloudwater requirement for cloud with new precip
         772             : !  Min CF is set at top for LE, at end for other levels
         773             : !  CWMIN is only needed for new precip formation - do not need for RNEW<0
         774             : !-----------------------------------------------------------------------
         775  2154041850 :              if( cfxx(l) == zero ) then
         776           0 :                if ( do_diag ) then
         777           0 :                  write(*,*) 'cfxx(l) == zero',l
         778           0 :                  write(*,*) qttjfl(:,n)
         779           0 :                  write(*,*) qm(:)
         780           0 :                  write(*,*) pofl(:)
         781           0 :                  write(*,*) delz(:)
         782           0 :                  write(*,*) rls(:)
         783           0 :                  write(*,*) clwc(:)
         784           0 :                  write(*,*) ciwc(:)
         785           0 :                  write(*,*) cfr(:)
         786           0 :                  write(*,*) tem(:)
         787           0 :                  write(*,*) evaprate(:)
         788           0 :                  write(*,*) hstar(:,n)
         789             :                end if
         790             : !
         791             : ! if we are here,, that means that there is
         792             : ! a inconsistency and this will lead to a division
         793             : ! by 0 later on! This column should then be skipped
         794             : !
         795           0 :                QTTJFL(:lpar,n) = QTT(:lpar)
         796             :                cycle species_loop
         797             : !
         798             : !              call endrun()
         799             : !
         800             :              endif
         801  2154041850 :              rnew_flag(l) = .true.
         802  2154041850 :              CLWX = max( CLWC(L)+CIWC(L),CWMIN*CFXX(L) )
         803  2154041850 :              if( is_hno3 ) then
         804           0 :                clwx_wrk(l) = clwx
         805             :              endif
         806             : !-----------------------------------------------------------------------
         807             : !  Area of old cloud and new cloud
         808             : !-----------------------------------------------------------------------
         809  2154041850 :              FCXA = FCA
         810  2154041850 :              FCXB = max( zero,CFXX(L)-FCXA )
         811             : !-----------------------------------------------------------------------
         812             : !                           ICE
         813             : !  For ice and mixed phase, grow precip in old cloud by riming
         814             : !  Use only portion of cloudwater in old cloud fraction
         815             : !  and rain above old cloud fraction
         816             : !  COLEFF from Lohmann and Roeckner (1996), Loss rate from Rotstayn (1997)
         817             : !-----------------------------------------------------------------------
         818             : is_freezing : &
         819  2154041850 :              if( freezing(l) ) then
         820  1160372665 :                COLEFFSNOW = exp( 2.5e-2_r8*(TEM(L) - TICE) )
         821  1160372665 :                if( TEM(L) <= TFROZ ) then
         822             :                  RHOSNOW = RHOSNOWFIX
         823             :                else
         824   504179785 :                  RHOSNOW = 0.303_r8*(TEM(L) - TFROZ)*RHOSNOWFIX
         825             :                endif
         826  1160372665 :                if( FCXA > zero ) then
         827  1041872415 :                  if( DCA > zero ) then
         828             :                    DELTARIMEMASS = CLWX*QM(L)*(FCXA/CFXX(L))* &
         829  1041872415 :                      (one - exp( (-COLEFFSNOW/(DCA*1.e-3_r8))*((RCA)/(2._r8*RHOSNOW))*DTSCAV ))   !uses GBA R
         830             :                  else
         831             :                    DELTARIMEMASS = zero
         832             :                  endif
         833             :                else
         834             :                  DELTARIMEMASS = zero
         835             :                endif
         836             : !-----------------------------------------------------------------------
         837             : !  Increase in precip rate due to riming (kg/m2/s):
         838             : !  Limit to total increase in R in cloud
         839             : !-----------------------------------------------------------------------
         840  1160372665 :                if( FCXA > zero ) then
         841  1041872415 :                  DELTARIME = min( RNEW/FCXA,DELTARIMEMASS/(FCXA*GAREA*DTSCAV) ) !GBA
         842             :                else
         843   118500250 :                  DELTARIME = zero
         844             :                endif
         845  1160372665 :                if( is_hno3 ) then
         846           0 :                  deltarime_wrk(l) = deltarime
         847             :                endif
         848             : !-----------------------------------------------------------------------
         849             : !  Find diameter of rimed precip, must be at least .1mm
         850             : !-----------------------------------------------------------------------
         851  1160372665 :                if( RCA > zero ) then
         852  1041872415 :                  DOR = max( DMIN,(((RCA+DELTARIME)/RCA)**VOLPOW)*DCA )
         853             :                else
         854   118500250 :                  DOR = zero
         855             :                endif
         856             : !-----------------------------------------------------------------------
         857             : !  If there is some in-cloud precip left, we have new precip formation
         858             : !  Will be spread over whole cloud fraction
         859             : !-----------------------------------------------------------------------
         860             : !  Calculate precip rate in old and new cloud fractions
         861             : !-----------------------------------------------------------------------
         862  1160372665 :                RPRECIP = (RNEW-(DELTARIME*FCXA))/CFXX(L) !kg/m2/s    !GBA
         863             : !-----------------------------------------------------------------------
         864             : !  Calculate precip rate in old and new cloud fractions
         865             : !-----------------------------------------------------------------------
         866  1160372665 :                RCXA = RCA + DELTARIME + RPRECIP          !kg/m2/s GBA
         867  1160372665 :                RCXB = RPRECIP                            !kg/m2/s GBA
         868             :
         869             : !-----------------------------------------------------------------------
         870             : !  Find diameter of new precip from empirical relation using Rprecip
         871             : !  in given area of box- use density of water, not snow, to convert kg/s
         872             : !  to mm/s -> as given in Field and Heymsfield
         873             : !  Also calculate diameter of mixed precip,DCXA, from empirical relation
         874             : !  using total R in FCXA - this will give larger particles than averaging DOR and
         875             : !  DNEW in the next level
         876             : !  DNEW and DCXA must be at least .1mm
         877             : !-----------------------------------------------------------------------
         878  1160372665 :                if( RPRECIP > zero ) then
         879  1024096410 :                  WEMP = (CLWX*QM(L))/(GAREA*CFXX(L)*DELZ(L)) !kg/m3
         880  1024096410 :                  REMP = RPRECIP/((RHORAIN/1.e3_r8))             !mm/s local
         881  1024096410 :                  DNEW = DEMPIRICAL( WEMP, REMP )
         882  1024096410 :                  if ( debug ) then
         883           0 :                    if( is_hno3 .and. l >= 15 ) then
         884           0 :                      write(*,*) ' '
         885           0 :                      write(*,*) 'washout: wemp,remp.dnew @ l = ',l
         886           0 :                      write(*,'(1p,3g15.7)') wemp,remp,dnew
         887           0 :                      write(*,*) ' '
         888             :                    endif
         889             :                  endif
         890  1024096410 :                  DNEW = max( DMIN,DNEW )
         891  1024096410 :                  if( FCXB > zero ) then
         892   255081490 :                    DCXB = DNEW
         893             :                  else
         894   769014920 :                    DCXB = zero
         895             :                  endif
         896             :                else
         897   136276255 :                  DCXB = zero
         898             :                endif
         899             :
         900  1160372665 :                if( FCXA > zero ) then
         901  1041872415 :                  WEMP = (CLWX*QM(L)*(FCXA/CFXX(L)))/(GAREA*FCXA*DELZ(L)) !kg/m3
         902  1041872415 :                  REMP = RCXA/((RHORAIN/1.e3_r8))                         !mm/s local
         903  1041872415 :                  DEMP = DEMPIRICAL( WEMP, REMP )
         904  1041872415 :                  DCXA = ((RCA+DELTARIME)/RCXA)*DOR + (RPRECIP/RCXA)*DNEW
         905  1041872415 :                  DCXA = max( DEMP,DCXA )
         906  1041872415 :                  DCXA = max( DMIN,DCXA )
         907             :                else
         908   118500250 :                  DCXA = zero
         909             :                endif
         910  1160372665 :                if ( debug ) then
         911           0 :                  if( is_hno3 .and. l >= 15 ) then
         912           0 :                    write(*,*) ' '
         913           0 :                    write(*,*) 'washout: rca,rcxa,deltarime,dor,rprecip,dnew @ l = ',l
         914           0 :                    write(*,'(1p,6g15.7)') rca,rcxa,deltarime,dor,rprecip,dnew
         915           0 :                    write(*,*) 'washout: dcxa,dcxb,wemp,remp,demp'
         916           0 :                    write(*,'(1p,5g15.7)') dcxa,dcxb,wemp,remp,demp
         917           0 :                    write(*,*) ' '
         918             :                  end if
         919             :                endif
         920             :
         921  1160372665 :                if( QTT(L) > zero ) then
         922             : !-----------------------------------------------------------------------
         923             : !                       ICE SCAVENGING
         924             : !-----------------------------------------------------------------------
         925             : !  For ice, rainout only hno3/aerosols using new precip
         926             : !  Tracer dissolved given by Kaercher and Voigt (2006) for T<258K
         927             : !  For T>258K, use Henry's Law with Retention coefficient
         928             : !  Rain out in whole CF
         929             : !-----------------------------------------------------------------------
         930  1160372665 :                  if( RPRECIP > zero ) then
         931  1024096410 :                    if( LICETYP == 1 ) then
         932           0 :                      RRAIN = RPRECIP*GAREA                                  !kg/s local
         933             :                      call DISGAS( CLWX, CFXX(L), TCMASS(N), HSTAR(L,N), &
         934             :                                   TEM(L),POFL(L),QM(L),                 &
         935           0 :                                   QTT(L)*CFXX(L),QTDISCF )
         936           0 :                      call RAINGAS( RRAIN, DTSCAV, CLWX, CFXX(L),        &
         937           0 :                                    QM(L), QTT(L), QTDISCF, QTRAIN )
         938           0 :                      WRK       = QTRAIN/CFXX(L)
         939           0 :                      QTRAINCXA = FCXA*WRK
         940           0 :                      QTRAINCXB = FCXB*WRK
         941  1024096410 :                    elseif( LICETYP == 2 ) then
         942  1024096410 :                      QTRAINCXA = zero
         943  1024096410 :                      QTRAINCXB = zero
         944             :                    endif
         945  1024096410 :                    if( debug .and. is_hno3 .and. l == kdiag ) then
         946           0 :                      write(*,*) ' '
         947           0 :                      write(*,*) 'washout: Ice Scavenging'
         948           0 :                      write(*,*) 'washout: qtraincxa, qtraincxb, fcxa, fcxb, qt_rain, cfxx(l), wrk @ level = ',l
         949           0 :                      write(*,'(1p,7g15.7)') qtraincxa, qtraincxb, fcxa, fcxb, qt_rain(l), cfxx(l), wrk
         950           0 :                      write(*,*) ' '
         951             :                    endif
         952             :                  endif
         953             : !-----------------------------------------------------------------------
         954             : !  For ice, accretion removal for hno3 and aerosols is propotional to riming,
         955             : !  no accretion removal for gases
         956             : !  remove only in mixed portion of cloud
         957             : !  Limit DELTARIMEMASS to RNEW*DTSCAV for ice - evaporation of rimed ice to match
         958             : !  RNEW precip rate would result in HNO3 escaping from ice (no trapping)
         959             : !-----------------------------------------------------------------------
         960  1160372665 :                  if( DELTARIME > zero ) then
         961   963070675 :                    if( LICETYP == 1 ) then
         962           0 :                      if( TEM(L) <= TFROZ ) then
         963             :                        RHOSNOW = RHOSNOWFIX
         964             :                      else
         965           0 :                        RHOSNOW = 0.303_r8*(TEM(L) - TFROZ)*RHOSNOWFIX
         966             :                      endif
         967           0 :                      QTCXA = QTT(L)*FCXA
         968             :                      call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N),   &
         969             :                                   HSTAR(L,N), TEM(L), POFL(L),            &
         970           0 :                                   QM(L), QTCXA, QTDISRIME )
         971           0 :                      QTDISSTAR = (QTDISRIME*QTCXA)/(QTDISRIME + QTCXA)
         972           0 :                      if ( debug ) then
         973           0 :                        if( is_hno3 .and. l >= 15 ) then
         974           0 :                          write(*,*) ' '
         975           0 :                          write(*,*) 'washout: fcxa,dca,rca,qtdisstar @ l = ',l
         976           0 :                          write(*,'(1p,4g15.7)') fcxa,dca,rca,qtdisstar
         977           0 :                          write(*,*) ' '
         978             :                        endif
         979             :                      endif
         980             :                      QTRIMECXA = QTCXA*                             &
         981             :                         (one - exp((-COLEFFSNOW/(DCA*1.e-3_r8))*       &
         982             :                         (RCA/(2._r8*RHOSNOW))*                         &  !uses GBA R
         983           0 :                         (QTDISSTAR/QTCXA)*DTSCAV))
         984             :                      QTRIMECXA = min( QTRIMECXA, &
         985           0 :                         ((RNEW*GAREA*DTSCAV)/(CLWX*QM(L)*(FCXA/CFXX(L))))*QTDISSTAR)
         986   963070675 :                    elseif( LICETYP == 2 ) then
         987   963070675 :                      QTRIMECXA = zero
         988             :                    endif
         989             :                  endif
         990             :                else
         991           0 :                  QTRAINCXA = zero
         992           0 :                  QTRAINCXB = zero
         993           0 :                  QTRIMECXA = zero
         994             :                endif
         995             : !-----------------------------------------------------------------------
         996             : !  For ice, no washout in interstitial cloud air
         997             : !-----------------------------------------------------------------------
         998  1160372665 :                QTWASHCXA = zero
         999  1160372665 :                QTEVAPCXA = zero
        1000             :
        1001             : !-----------------------------------------------------------------------
        1002             : !                      RAIN
        1003             : !  For rain, accretion increases rain rate but diameter remains constant
        1004             : !  Diameter is 4mm (not used)
        1005             : !-----------------------------------------------------------------------
        1006             :              else is_freezing
        1007   993669185 :                if( FCXA > zero ) then
        1008             :                  DELTARIMEMASS = (CLWX*QM(L))*(FCXA/CFXX(L))*           &
        1009   943815645 :                    (one - exp( -0.24_r8*COLEFFRAIN*((RCA)**0.75_r8)*DTSCAV ))  !local
        1010             :                else
        1011             :                  DELTARIMEMASS = zero
        1012             :                endif
        1013             : !-----------------------------------------------------------------------
        1014             : !  Increase in precip rate due to riming (kg/m2/s):
        1015             : !  Limit to total increase in R in cloud
        1016             : !-----------------------------------------------------------------------
        1017   993669185 :                if( FCXA > zero ) then
        1018   943815645 :                  DELTARIME = min( RNEW/FCXA,DELTARIMEMASS/(FCXA*GAREA*DTSCAV) ) !GBA
        1019             :                else
        1020    49853540 :                  DELTARIME = zero
        1021             :                endif
        1022             : !-----------------------------------------------------------------------
        1023             : !  If there is some in-cloud precip left, we have new precip formation
        1024             : !-----------------------------------------------------------------------
        1025   993669185 :                RPRECIP = (RNEW-(DELTARIME*FCXA))/CFXX(L)       !GBA
        1026             :
        1027   993669185 :                RCXA = RCA + DELTARIME + RPRECIP            !kg/m2/s GBA
        1028   993669185 :                RCXB = RPRECIP                              !kg/m2/s GBA
        1029   993669185 :                DCXA = FOUR
        1030   993669185 :                if( FCXB > zero ) then
        1031   153130420 :                  DCXB = FOUR
        1032             :                else
        1033   840538765 :                  DCXB = zero
        1034             :                endif
        1035             : !-----------------------------------------------------------------------
        1036             : !                         RAIN SCAVENGING
        1037             : !  For rain, rainout both hno3/aerosols and gases using new precip
        1038             : !-----------------------------------------------------------------------
        1039   993669185 :                if( QTT(L) > zero ) then
        1040   993669185 :                  if( RPRECIP > zero ) then
        1041   804624160 :                    RRAIN = (RPRECIP*GAREA) !kg/s local
        1042             :                    call DISGAS( CLWX, CFXX(L), TCMASS(N), HSTAR(L,N), &
        1043             :                                 TEM(L), POFL(L), QM(L),               &
        1044   804624160 :                                 QTT(L)*CFXX(L), QTDISCF )
        1045   804624160 :                    call RAINGAS( RRAIN, DTSCAV, CLWX, CFXX(L),        &
        1046  1609248320 :                                  QM(L), QTT(L), QTDISCF, QTRAIN )
        1047   804624160 :                    WRK       = QTRAIN/CFXX(L)
        1048   804624160 :                    QTRAINCXA = FCXA*WRK
        1049   804624160 :                    QTRAINCXB = FCXB*WRK
        1050   804624160 :                    if( debug .and. is_hno3 .and. l == kdiag ) then
        1051           0 :                      write(*,*) ' '
        1052           0 :                      write(*,*) 'washout: Rain Scavenging'
        1053           0 :                      write(*,*) 'washout: qtraincxa, qtraincxb, fcxa, fcxb, qt_rain, cfxx(l), wrk @ level = ',l
        1054           0 :                      write(*,'(1p,7g15.7)') qtraincxa, qtraincxb, fcxa, fcxb, qt_rain(l), cfxx(l), wrk
        1055           0 :                      write(*,*) ' '
        1056             :                    endif
        1057             :                  endif
        1058             : !-----------------------------------------------------------------------
        1059             : !  For rain, accretion removal is propotional to riming
        1060             : !  caclulate for hno3/aerosols and gases
        1061             : !  Remove only in mixed portion of cloud
        1062             : !  Limit DELTARIMEMASS to RNEW*DTSCAV
        1063             : !-----------------------------------------------------------------------
        1064   993669185 :                  if( DELTARIME > zero ) then
        1065   943815350 :                    QTCXA = QTT(L)*FCXA
        1066             :                    call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N),    &
        1067             :                                 HSTAR(L,N), TEM(L), POFL(L),             &
        1068   943815350 :                                 QM(L), QTCXA, QTDISRIME )
        1069   943815350 :                    QTDISSTAR = (QTDISRIME*QTCXA)/(QTDISRIME + QTCXA)
        1070             :                    QTRIMECXA = QTCXA*                              &
        1071             :                       (one - exp(-0.24_r8*COLEFFRAIN*                 &
        1072             :                       ((RCA)**0.75_r8)*                               & !local
        1073   943815350 :                       (QTDISSTAR/QTCXA)*DTSCAV))
        1074             :                    QTRIMECXA = min( QTRIMECXA, &
        1075   943815350 :                       ((RNEW*GAREA*DTSCAV)/(CLWX*QM(L)*(FCXA/CFXX(L))))*QTDISSTAR)
        1076             :                  else
        1077    49853835 :                    QTRIMECXA = zero
        1078             :                  endif
        1079             :                else
        1080           0 :                  QTRAINCXA = zero
        1081           0 :                  QTRAINCXB = zero
        1082           0 :                  QTRIMECXA = zero
        1083             :                endif
        1084             : !-----------------------------------------------------------------------
        1085             : !  For rain, washout gases and HNO3/aerosols using rain from above old cloud
        1086             : !  Washout for HNO3/aerosols is only on non-dissolved portion, impaction-style
        1087             : !  Washout for gases is on non-dissolved portion, limited by QTTOP+QTRIME
        1088             : !-----------------------------------------------------------------------
        1089   993669185 :                if( RCA > zero ) then
        1090   943815645 :                  QTPRECIP = FCXA*QTT(L) - QTDISRIME
        1091   943815645 :                  if( HSTAR(L,N) > 1.e4_r8 ) then
        1092   612669258 :                    if( QTPRECIP > zero ) then
        1093   210573067 :                      QTWASHCXA = QTPRECIP*(one - exp( -0.24_r8*COLEFFAER*((RCA)**0.75_r8)*DTSCAV ))   !local
        1094             :                    else
        1095   402096191 :                      QTWASHCXA = zero
        1096             :                    endif
        1097   612669258 :                    QTEVAPCXA = zero
        1098             :                  else
        1099   331146387 :                    RWASH = RCA*GAREA                                !kg/s local
        1100   331146387 :                    if( QTPRECIP > zero ) then
        1101             :                      call WASHGAS( RWASH, FCA, DTSCAV, QTTOPCA+QTRIMECXA, &
        1102             :                                    HSTAR(L,N), TEM(L), POFL(L),           &
        1103   331146387 :                                    QM(L), QTPRECIP, QTWASHCXA, QTEVAPCXA )
        1104             :                    else
        1105           0 :                      QTWASHCXA = zero
        1106           0 :                      QTEVAPCXA = zero
        1107             :                    endif
        1108             :                  endif
        1109             :                endif
        1110             :              endif is_freezing
        1111             : !-----------------------------------------------------------------------
        1112             : !  If RNEW<O, confine precip to area of cloud above
        1113             : !  FCXA does not require a minimum (could be zero if R(L).le.what
        1114             : !  evaporated in ambient)
        1115             : !-----------------------------------------------------------------------
        1116             :            else has_rnew
        1117   810310635 :              CLWX = CLWC(L) + CIWC(L)
        1118   810310635 :              if( is_hno3 ) then
        1119           0 :                clwx_wrk(l) = clwx
        1120             :              endif
        1121   810310635 :              FCXA = FCA
        1122   810310635 :              FCXB = max( zero,CFXX(L)-FCXA )
        1123   810310635 :              RCXB = zero
        1124   810310635 :              DCXB = zero
        1125   810310635 :              QTRAINCXA = zero
        1126   810310635 :              QTRAINCXB = zero
        1127   810310635 :              QTRIMECXA = zero
        1128             :
        1129             : !-----------------------------------------------------------------------
        1130             : !  Put rain into cloud up to RCA so that we evaporate
        1131             : !  from ambient first
        1132             : !  Adjust ambient to try to match RLS(L)
        1133             : !  If no cloud, RAX=R(L)
        1134             : !-----------------------------------------------------------------------
        1135   810310635 :              if( FCXA > zero ) then
        1136   795669630 :                RCXA = min( RCA,RLS(L)/(GAREA*FCXA) )     !kg/m2/s  GBA
        1137   795669630 :                if( FAX > zero .and. ((RCXA+1.e-12_r8) < RLS(L)/(GAREA*FCXA)) ) then
        1138   397757880 :                  RAXADJF = RLS(L)/GAREA - RCXA*FCXA
        1139   397757880 :                  RAMPCT = RAXADJF/(RAX*FAX)
        1140   397757880 :                  FAXADJ = RAMPCT*FAX
        1141   397757880 :                  if( FAXADJ > zero ) then
        1142   397757880 :                    RAXADJ = RAXADJF/FAXADJ
        1143             :                  else
        1144             :                    RAXADJ = zero
        1145             :                  endif
        1146             :                else
        1147             :                  RAXADJ = zero
        1148             :                  RAMPCT = zero
        1149             :                  FAXADJ = zero
        1150             :                endif
        1151             :              else
        1152    14641005 :                RCXA = zero
        1153    14641005 :                if( FAX > zero ) then
        1154    14641005 :                  RAXADJF = RLS(L)/GAREA
        1155    14641005 :                  RAMPCT = RAXADJF/(RAX*FAX)
        1156    14641005 :                  FAXADJ = RAMPCT*FAX
        1157    14641005 :                  if( FAXADJ > zero ) then
        1158    14641005 :                    RAXADJ = RAXADJF/FAXADJ
        1159             :                  else
        1160             :                    RAXADJ = zero
        1161             :                  endif
        1162             :                else
        1163             :                  RAXADJ = zero
        1164             :                  RAMPCT = zero
        1165             :                  FAXADJ = zero
        1166             :                endif
        1167             :              endif
        1168             :
        1169   810310635 :              QTEVAPAXP = min( QTTOPAA,QTTOPAA - (RAMPCT*(QTTOPAA-QTEVAPAXP)) )
        1170   810310635 :              FAX = FAXADJ
        1171   810310635 :              RAX = RAXADJ
        1172   810310635 :              if ( debug ) then
        1173           0 :                if( (l == 3 .or. l == 2) ) then
        1174           0 :                  write(*,*) 'washout: l,fcxa,fax = ',l,fcxa,fax
        1175             :                endif
        1176             :              endif
        1177             :
        1178             : !-----------------------------------------------------------------------
        1179             : !                IN-CLOUD EVAPORATION/WASHOUT
        1180             : !  If precip out the bottom of the cloud is 0, evaporate everything
        1181             : !  If there is no cloud, QTTOPCA=0, so nothing happens
        1182             : !-----------------------------------------------------------------------
        1183   810310635 :              if( RCXA <= zero ) then
        1184    14641005 :                QTEVAPCXA = QTTOPCA
        1185    14641005 :                RCXA = zero
        1186    14641005 :                DCXA = zero
        1187             :              else
        1188             : !-----------------------------------------------------------------------
        1189             : !  If rain out the bottom of the cloud is >0 (but .le. RCA):
        1190             : !  For ice, decrease particle size,
        1191             : !  no washout
        1192             : !  no evap for non-ice gases (b/c there is nothing in ice)
        1193             : !  T<Tmix,release hno3& aerosols
        1194             : !  release is amount dissolved in ice mass released
        1195             : !  T>Tmix, hno3&aerosols are incorporated into ice structure:
        1196             : !  do not release
        1197             : !  For rain, assume full evaporation of some raindrops
        1198             : !  proportional evaporation for all species
        1199             : !  washout for gases using Rbot
        1200             : !  impact washout for hno3/aerosol portion in gas phase
        1201             : !-----------------------------------------------------------------------
        1202             : !              if (TEM(L) < TICE ) then
        1203             : is_freezing_a : &
        1204   795669630 :                if( freezing(l) ) then
        1205   157963950 :                  QTWASHCXA = zero
        1206   157963950 :                  DCXA = ((RCXA/RCA)**VOLPOW)*DCA
        1207   157963950 :                  if( LICETYP == 1 ) then
        1208           0 :                    if( TEM(L) <= TMIX ) then
        1209           0 :                      MASSLOSS = (RCA-RCXA)*FCXA*GAREA*DTSCAV
        1210             : !-----------------------------------------------------------------------
        1211             : !  note-QTT doesn't matter b/c T<258K
        1212             : !-----------------------------------------------------------------------
        1213             :                      call DISGAS( (MASSLOSS/QM(L)), FCXA, TCMASS(N),   &
        1214             :                                    HSTAR(L,N), TEM(L), POFL(L),        &
        1215           0 :                                    QM(L), QTT(L), QTEVAPCXA )
        1216           0 :                      QTEVAPCXA = min( QTTOPCA,QTEVAPCXA )
        1217             :                    else
        1218           0 :                      QTEVAPCXA = zero
        1219             :                    endif
        1220   157963950 :                  elseif( LICETYP == 2 ) then
        1221   157963950 :                    QTEVAPCXA = zero
        1222             :                  endif
        1223             :                else is_freezing_a
        1224   637705680 :                  QTEVAPCXAP = (RCA - RCXA)/RCA*QTTOPCA
        1225   637705680 :                  DCXA = FOUR
        1226   637705680 :                  QTCXA = FCXA*QTT(L)
        1227   637705680 :                  if( HSTAR(L,N) > 1.e4_r8 ) then
        1228   388125220 :                    if( QTT(L) > zero ) then
        1229             :                      call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N),   &
        1230             :                                   HSTAR(L,N), TEM(L), POFL(L),            &
        1231   388125220 :                                   QM(L), QTCXA, QTDISCXA )
        1232   388125220 :                      if( QTCXA > QTDISCXA ) then
        1233   340926283 :                        QTWASHCXA = (QTCXA - QTDISCXA)*(one - exp( -0.24_r8*COLEFFAER*((RCXA)**0.75_r8)*DTSCAV )) !local
        1234             :                      else
        1235    47198937 :                        QTWASHCXA = zero
        1236             :                      endif
        1237   388125220 :                      QTEVAPCXAW = zero
        1238             :                    else
        1239           0 :                      QTWASHCXA  = zero
        1240           0 :                      QTEVAPCXAW = zero
        1241             :                    endif
        1242             :                  else
        1243   249580460 :                    RWASH = RCXA*GAREA                         !kg/s local
        1244             :                    call WASHGAS( RWASH, FCXA, DTSCAV, QTTOPCA, HSTAR(L,N), &
        1245             :                                  TEM(L), POFL(L), QM(L),                   &
        1246   249580460 :                                  QTCXA-QTDISCXA, QTWASHCXA, QTEVAPCXAW )
        1247             :                  endif
        1248   637705680 :                  QTEVAPCXA = QTEVAPCXAP + QTEVAPCXAW
        1249             :                endif is_freezing_a
        1250             :              endif
        1251             :            endif has_rnew
        1252             :
        1253             : !-----------------------------------------------------------------------
        1254             : !                 AMBIENT WASHOUT
        1255             : !  Ambient precip is finalized - if it is rain, washout
        1256             : !  no ambient washout for ice, since gases are in vapor phase
        1257             : !-----------------------------------------------------------------------
        1258  2964352485 :            if( RAX > zero ) then
        1259  1420544355 :              if( .not. freezing(l) ) then
        1260   973697745 :                QTAX = FAX*QTT(L)
        1261   973697745 :                if( HSTAR(L,N) > 1.e4_r8 ) then
        1262             :                  QTWASHAX = QTAX*                        &
        1263             :                     (one - exp(-0.24_r8*COLEFFAER*       &
        1264   618947865 :                    ((RAX)**0.75_r8)*DTSCAV))  !local
        1265   618947865 :                  QTEVAPAXW = zero
        1266             :                else
        1267   354749880 :                  RWASH = RAX*GAREA   !kg/s local
        1268             :                  call WASHGAS( RWASH, FAX, DTSCAV, QTTOPAA, HSTAR(L,N), &
        1269             :                                TEM(L), POFL(L), QM(L), QTAX,            &
        1270   354749880 :                                QTWASHAX, QTEVAPAXW )
        1271             :                endif
        1272             :              else
        1273   446846610 :                QTEVAPAXW = zero
        1274   446846610 :                QTWASHAX  = zero
        1275             :              endif
        1276             :            else
        1277  1543808130 :              QTEVAPAXW = zero
        1278  1543808130 :              QTWASHAX  = zero
        1279             :            endif
        1280  2964352485 :            QTEVAPAX = QTEVAPAXP + QTEVAPAXW
        1281             :
        1282             : !-----------------------------------------------------------------------
        1283             : !                  END SCAVENGING
        1284             : !  Require CF if our ambient evaporation rate would give less
        1285             : !  precip than R from model.
        1286             : !-----------------------------------------------------------------------
        1287  2964352485 :            if( do_diag .and. is_hno3 ) then
        1288           0 :              rls_wrk(l) = rls(l)/garea
        1289           0 :              rca_wrk(l) = rca
        1290           0 :              fca_wrk(l) = fca
        1291           0 :              rcxa_wrk(l) = rcxa
        1292           0 :              fcxa_wrk(l) = fcxa
        1293           0 :              rcxb_wrk(l) = rcxb
        1294           0 :              fcxb_wrk(l) = fcxb
        1295           0 :              rax_wrk(l,2) = rax
        1296           0 :              fax_wrk(l,2) = fax
        1297             :            endif
        1298             : upper_level : &
        1299  2964352485 :            if( L > 1 ) then
        1300  2870933190 :              if( CFR(LM1) >= CFMIN ) then
        1301  1177133540 :                CFXX(LM1) = CFR(LM1)
        1302             :              else
        1303  1693799650 :                if( adj_factor*RLSOG(LM1) >= (RCXA*FCXA + RCXB*FCXB + RAX*FAX)*(one - EVAPRATE(LM1)) ) then
        1304  1664488605 :                  CFXX(LM1) = CFMIN
        1305  1664488605 :                  cf_trigger(lm1) = .true.
        1306             :                else
        1307    29311045 :                  CFXX(LM1) = CFR(LM1)
        1308             :                endif
        1309  1693799650 :                if( is_hno3 .and. lm1 == kdiag .and. debug ) then
        1310           0 :                  write(*,*) ' '
        1311           0 :                  write(*,*) 'washout: rls,garea,rcxa,fcxa,rcxb,fcxb,rax,fax'
        1312           0 :                  write(*,'(1p,8g15.7)') rls(lm1),garea,rcxa,fcxa,rcxb,fcxb,rax,fax
        1313           0 :                  write(*,*) ' '
        1314             :                endif
        1315             :              endif
        1316             : !-----------------------------------------------------------------------
        1317             : !  Figure out what will go into ambient and cloud below
        1318             : !  Don't do for lowest level
        1319             : !-----------------------------------------------------------------------
        1320  2870933190 :              if( FAX > zero ) then
        1321  1361945000 :                AMPCT = max( zero,min( one,(CFXX(L) + FAX - CFXX(LM1))/FAX ) )
        1322  1361945000 :                AMCLPCT = one - AMPCT
        1323             :              else
        1324  1508988190 :                AMPCT   = zero
        1325  1508988190 :                AMCLPCT = zero
        1326             :              endif
        1327  2870933190 :              if( FCXB > zero ) then
        1328   419822435 :                CLNEWPCT = max( zero,min( (CFXX(LM1) - FCXA)/FCXB,one ) )
        1329   419822435 :                CLNEWAMPCT = one - CLNEWPCT
        1330             :              else
        1331  2451110755 :                CLNEWPCT   = zero
        1332  2451110755 :                CLNEWAMPCT = zero
        1333             :              endif
        1334  2870933190 :              if( FCXA > zero ) then
        1335  2688263890 :                CLOLDPCT = max( zero,min( CFXX(LM1)/FCXA,one ) )
        1336  2688263890 :                CLOLDAMPCT = one - CLOLDPCT
        1337             :              else
        1338   182669300 :                CLOLDPCT   = zero
        1339   182669300 :                CLOLDAMPCT = zero
        1340             :              endif
        1341             : !-----------------------------------------------------------------------
        1342             : !  Remix everything for the next level
        1343             : !-----------------------------------------------------------------------
        1344  2870933190 :              FCA = min( CFXX(LM1),FCXA*CLOLDPCT + CLNEWPCT*FCXB + AMCLPCT*FAX )
        1345  2870933190 :              if( FCA > zero ) then
        1346             : !-----------------------------------------------------------------------
        1347             : !  Maintain cloud core by reducing NC and AM area going into cloud below
        1348             : !-----------------------------------------------------------------------
        1349  2856291500 :                RCA = (RCXA*FCXA*CLOLDPCT + RCXB*FCXB*CLNEWPCT + RAX*FAX*AMCLPCT)/FCA
        1350  2856291500 :                if ( debug ) then
        1351           0 :                  if( is_hno3 ) then
        1352           0 :                    write(*,*) ' '
        1353           0 :                    write(*,*) 'washout: rcxa,fcxa,cloldpctrca,rca,fca,dcxa @ l = ',l
        1354           0 :                    write(*,'(1p,6g15.7)') rcxa,fcxa,cloldpct,rca,fca,dcxa
        1355           0 :                    write(*,*) 'washout: rcxb,fcxb,clnewpct,dcxb'
        1356           0 :                    write(*,'(1p,4g15.7)') rcxb,fcxb,clnewpct,dcxb
        1357           0 :                    write(*,*) 'washout: rax,fax,amclpct,dax'
        1358           0 :                    write(*,'(1p,4g15.7)') rax,fax,amclpct,dax
        1359           0 :                    write(*,*) ' '
        1360             :                  endif
        1361             :                endif
        1362             :
        1363  2856291500 :                if (RCA > zero) then
        1364             :                  DCA = (RCXA*FCXA*CLOLDPCT)/(RCA*FCA)*DCXA + &
        1365             :                        (RCXB*FCXB*CLNEWPCT)/(RCA*FCA)*DCXB + &
        1366  2856291500 :                        (RAX*FAX*AMCLPCT)/(RCA*FCA)*DAX
        1367             :                else
        1368           0 :                  DCA = zero
        1369           0 :                  FCA = zero
        1370             :                endif
        1371             :
        1372             :              else
        1373    14641690 :                FCA = zero
        1374    14641690 :                DCA = zero
        1375    14641690 :                RCA = zero
        1376             :              endif
        1377             :
        1378  2870933190 :              FAMA = FCXA + FCXB + FAX - CFXX(LM1)
        1379  2870933190 :              if( FAMA > zero ) then
        1380  1478504940 :                RAMA = (RCXA*FCXA*CLOLDAMPCT + RCXB*FCXB*CLNEWAMPCT + RAX*FAX*AMPCT)/FAMA
        1381  1478504940 :                if( RAMA > zero ) then
        1382             :                  DAMA = (RCXA*FCXA*CLOLDAMPCT)/(RAMA*FAMA)*DCXA +  &
        1383             :                         (RCXB*FCXB*CLNEWAMPCT)/(RAMA*FAMA)*DCXB +  &
        1384  1474153215 :                         (RAX*FAX*AMPCT)/(RAMA*FAMA)*DAX
        1385             :                else
        1386     4351725 :                   FAMA = zero
        1387     4351725 :                   DAMA = zero
        1388             :                endif
        1389             :              else
        1390  1392428250 :                FAMA = zero
        1391  1392428250 :                DAMA = zero
        1392  1392428250 :                RAMA = zero
        1393             :              endif
        1394             :            else upper_level
        1395    93419295 :              AMPCT      = zero
        1396    93419295 :              AMCLPCT    = zero
        1397    93419295 :              CLNEWPCT   = zero
        1398    93419295 :              CLNEWAMPCT = zero
        1399    93419295 :              CLOLDPCT   = zero
        1400    93419295 :              CLOLDAMPCT = zero
        1401             :            endif upper_level
        1402             :          else has_rls
        1403  7788883515 :            RNEW = zero
        1404  7788883515 :            QTEVAPCXA = QTTOPCA
        1405  7788883515 :            QTEVAPAX = QTTOPAA
        1406  7788883515 :            if( L > 1 ) then
        1407  7765419810 :              if( RLS(LM1) > zero ) then
        1408   168353790 :                CFXX(LM1) = max( CFMIN,CFR(LM1) )
        1409             : !              if( CFR(LM1) >= CFMIN ) then
        1410             : !                CFXX(LM1) = CFR(LM1)
        1411             : !              else
        1412             : !                CFXX(LM1) = CFMIN
        1413             : !              endif
        1414             :              else
        1415  7597066020 :                CFXX(LM1) = CFR(LM1)
        1416             :              endif
        1417             :            endif
        1418  7788883515 :            AMPCT      = zero
        1419  7788883515 :            AMCLPCT    = zero
        1420  7788883515 :            CLNEWPCT   = zero
        1421  7788883515 :            CLNEWAMPCT = zero
        1422  7788883515 :            CLOLDPCT   = zero
        1423  7788883515 :            CLOLDAMPCT = zero
        1424  7788883515 :            RCA        = zero
        1425  7788883515 :            RAMA       = zero
        1426  7788883515 :            FCA        = zero
        1427  7788883515 :            FAMA       = zero
        1428  7788883515 :            DCA        = zero
        1429  7788883515 :            DAMA       = zero
        1430             :          endif has_rls
        1431             :
        1432 10753236000 :          if( do_diag .and. is_hno3 ) then
        1433           0 :            fama_wrk(l) = fama
        1434           0 :            rama_wrk(l) = rama
        1435             :          endif
        1436             : !-----------------------------------------------------------------------
        1437             : !  Net loss can not exceed QTT in each region
        1438             : !-----------------------------------------------------------------------
        1439 10753236000 :          QTNETLCXA = QTRAINCXA + QTRIMECXA + QTWASHCXA - QTEVAPCXA
        1440 10753236000 :          QTNETLCXA = min( QTT(L)*FCXA,QTNETLCXA )
        1441             :
        1442 10753236000 :          QTNETLCXB =QTRAINCXB
        1443 10753236000 :          QTNETLCXB = min( QTT(L)*FCXB,QTNETLCXB )
        1444             :
        1445 10753236000 :          QTNETLAX = QTWASHAX - QTEVAPAX
        1446 10753236000 :          QTNETLAX = min( QTT(L)*FAX,QTNETLAX )
        1447             :
        1448 10753236000 :          QTTNEW(L) = QTT(L) - (QTNETLCXA + QTNETLCXB + QTNETLAX)
        1449             :
        1450 10753236000 :          if( do_diag .and. is_hno3 ) then
        1451           0 :            qt_rain(l) = qtraincxa + qtraincxb
        1452           0 :            qt_rime(l) = qtrimecxa
        1453           0 :            qt_wash(l) = qtwashcxa + qtwashax
        1454           0 :            qt_evap(l) = qtevapcxa + qtevapax
        1455           0 :            frc(l,1) = qtnetlcxa
        1456           0 :            frc(l,2) = qtnetlcxb
        1457           0 :            frc(l,3) = qtnetlax
        1458             :          endif
        1459 10753236000 :          if( debug .and. is_hno3 .and. l == kdiag ) then
        1460           0 :            write(*,*) ' '
        1461           0 :            write(*,*) 'washout: qtraincxa, qtraincxb, qtrimecxa @ level = ',l
        1462           0 :            write(*,'(1p,3g15.7)') qtraincxa, qtraincxb, qtrimecxa
        1463           0 :            write(*,*) ' '
        1464             :          endif
        1465 10753236000 :          if ( debug ) then
        1466           0 :            if( (l == 3 .or. l == 2) ) then
        1467           0 :              write(*,*) 'washout: hno3, hno3, qtnetlca,b, qtnetlax @ level = ',l
        1468           0 :              write(*,'(1p,5g15.7)') qttnew(l), qtt(l), qtnetlcxa, qtnetlcxb, qtnetlax
        1469           0 :              write(*,*) 'washout: qtwashax, qtevapax,fax,fama'
        1470           0 :              write(*,'(1p,5g15.7)') qtwashax, qtevapax, fax, fama
        1471             :            endif
        1472             :          endif
        1473             :
        1474 10753236000 :          QTTOPCAX = (QTTOPCA + QTNETLCXA)*CLOLDPCT + QTNETLCXB*CLNEWPCT + (QTTOPAA + QTNETLAX)*AMCLPCT
        1475 10753236000 :          QTTOPAAX = (QTTOPCA + QTNETLCXA)*CLOLDAMPCT + QTNETLCXB*CLNEWAMPCT + (QTTOPAA + QTNETLAX)*AMPCT
        1476 10753236000 :          QTTOPCA = QTTOPCAX
        1477 10870119000 :          QTTOPAA = QTTOPAAX
        1478             :        end do level_loop
        1479             :
        1480   116883000 :        if ( debug ) then
        1481           0 :          if( is_hno3 ) then
        1482           0 :            write(*,*) ' '
        1483           0 :            write(*,*) 'washout: clwx_wrk'
        1484           0 :            write(*,'(1p,5g15.7)') clwx_wrk(1:le)
        1485           0 :            write(*,*) 'washout: cfr'
        1486           0 :            write(*,'(1p,5g15.7)') cfr(1:le)
        1487           0 :            write(*,*) 'washout: cfxx'
        1488           0 :            write(*,'(1p,5g15.7)') cfxx(1:le)
        1489           0 :            write(*,*) 'washout: cf trigger'
        1490           0 :            write(*,'(10l4)') cf_trigger(1:le)
        1491           0 :            write(*,*) 'washout: evaprate'
        1492           0 :            write(*,'(1p,5g15.7)') evaprate(1:le)
        1493           0 :            write(*,*) 'washout: rls'
        1494           0 :            write(*,'(1p,5g15.7)') rls(1:le)
        1495           0 :            write(*,*) 'washout: rls/garea'
        1496           0 :            write(*,'(1p,5g15.7)') rls_wrk(1:le)
        1497           0 :            write(*,*) 'washout: rnew_wrk'
        1498           0 :            write(*,'(1p,5g15.7)') rnew_wrk(1:le)
        1499           0 :            write(*,*) 'washout: rnew_flag'
        1500           0 :            write(*,'(10l4)') rnew_flag(1:le)
        1501           0 :            write(*,*) 'washout: deltarime_wrk'
        1502           0 :            write(*,'(1p,5g15.7)') deltarime_wrk(1:le)
        1503           0 :            write(*,*) 'washout: rama_wrk'
        1504           0 :            write(*,'(1p,5g15.7)') rama_wrk(1:le)
        1505           0 :            write(*,*) 'washout: fama_wrk'
        1506           0 :            write(*,'(1p,5g15.7)') fama_wrk(1:le)
        1507           0 :            write(*,*) 'washout: rca_wrk'
        1508           0 :            write(*,'(1p,5g15.7)') rca_wrk(1:le)
        1509           0 :            write(*,*) 'washout: fca_wrk'
        1510           0 :            write(*,'(1p,5g15.7)') fca_wrk(1:le)
        1511           0 :            write(*,*) 'washout: rcxa_wrk'
        1512           0 :            write(*,'(1p,5g15.7)') rcxa_wrk(1:le)
        1513           0 :            write(*,*) 'washout: fcxa_wrk'
        1514           0 :            write(*,'(1p,5g15.7)') fcxa_wrk(1:le)
        1515           0 :            write(*,*) 'washout: rcxb_wrk'
        1516           0 :            write(*,'(1p,5g15.7)') rcxb_wrk(1:le)
        1517           0 :            write(*,*) 'washout: fcxb_wrk'
        1518           0 :            write(*,'(1p,5g15.7)') fcxb_wrk(1:le)
        1519           0 :            write(*,*) 'washout: rax1_wrk'
        1520           0 :            write(*,'(1p,5g15.7)') rax_wrk(1:le,1)
        1521           0 :            write(*,*) 'washout: fax1_wrk'
        1522           0 :            write(*,'(1p,5g15.7)') fax_wrk(1:le,1)
        1523           0 :            write(*,*) 'washout: rax2_wrk'
        1524           0 :            write(*,'(1p,5g15.7)') rax_wrk(1:le,2)
        1525           0 :            write(*,*) 'washout: fax2_wrk'
        1526           0 :            write(*,'(1p,5g15.7)') fax_wrk(1:le,2)
        1527           0 :            write(*,*) 'washout: rls_flag'
        1528           0 :            write(*,'(1p,10l4)') rls_flag(1:le)
        1529           0 :            write(*,*) 'washout: freezing'
        1530           0 :            write(*,'(1p,10l4)') freezing(1:le)
        1531           0 :            write(*,*) 'washout: qtnetlcxa'
        1532           0 :            write(*,'(1p,5g15.7)') frc(1:le,1)
        1533           0 :            write(*,*) 'washout: qtnetlcxb'
        1534           0 :            write(*,'(1p,5g15.7)') frc(1:le,2)
        1535           0 :            write(*,*) 'washout: qtnetlax'
        1536           0 :            write(*,'(1p,5g15.7)') frc(1:le,3)
        1537           0 :            write(*,*) ' '
        1538             :          endif
        1539             :        endif
        1540             : !-----------------------------------------------------------------------
        1541             : !  reload new tracer mass and rescale moments: check upper limits (LE)
        1542             : !-----------------------------------------------------------------------
        1543 10893495600 :        QTTJFL(:le,N) = QTTNEW(:le)
        1544             :
        1545             :      end do species_loop
        1546             : !
        1547    23376600 :      return
        1548     1489176 :    end subroutine washo
        1549             : !---------------------------------------------------------------------
        1550  2136564730 :       subroutine DISGAS (CLWX,CFX,MOLMASS,HSTAR,TM,PR,QM,QT,QTDIS)
        1551             : !---------------------------------------------------------------------
        1552             :       implicit none
        1553             :       real(r8), intent(in) :: CLWX,CFX    !cloud water,cloud fraction
        1554             :       real(r8), intent(in) :: MOLMASS     !molecular mass of tracer
        1555             :       real(r8), intent(in) :: HSTAR       !Henry's Law coeffs A*exp(-B/T)
        1556             :       real(r8), intent(in) :: TM          !temperature of box (K)
        1557             :       real(r8), intent(in) :: PR          !pressure of box (hPa)
        1558             :       real(r8), intent(in) :: QM          !air mass in box (kg)
        1559             :       real(r8), intent(in) :: QT          !tracer in box (kg)
        1560             :       real(r8), intent(out) :: QTDIS      !tracer dissolved in aqueous phase
        1561             :
        1562             :       real(r8)  MUEMP
        1563             :       real(r8), parameter :: INV298 = 1._r8/298._r8
        1564             :       real(r8), parameter  :: TMIX=258._r8
        1565             :       real(r8), parameter  :: RETEFF=0.5_r8
        1566             : !---Next calculate rate of uptake of tracer
        1567             :
        1568             : !---effective Henry's Law constant: H* = moles-T / liter-precip / press(atm-T)
        1569             : !---p(atm of tracer-T) = (QT/QM) * (.029/MolWt-T) * pressr(hPa)/1000
        1570             : !---limit temperature effects to T above freezing
        1571             : !----MU from fit to Kaercher and Voigt (2006)
        1572             :
        1573  2136564730 :       if(TM .ge. TICE) then
        1574  2136564730 :          QTDIS=(HSTAR*(QT/(QM*CFX))*0.029_r8*(PR/1.0e3_r8))*(CLWX*QM)
        1575           0 :       elseif (TM .le. TMIX) then
        1576           0 :          MUEMP=exp(-14.2252_r8+(1.55704e-1_r8*TM)-(7.1929e-4_r8*(TM**2.0_r8)))
        1577           0 :          QTDIS=MUEMP*(MOLMASS/18._r8)*(CLWX*QM)
        1578             :       else
        1579           0 :        QTDIS=RETEFF*((HSTAR*(QT/(QM*CFX))*0.029_r8*(PR/1.0e3_r8))*(CLWX*QM))
        1580             :       endif
        1581             :
        1582  2136564730 :       return
        1583             :       end subroutine DISGAS
        1584             :
        1585             : !-----------------------------------------------------------------------
        1586   804624160 :       subroutine RAINGAS (RRAIN,DTSCAV,CLWX,CFX,QM,QT,QTDIS,QTRAIN)
        1587             : !-----------------------------------------------------------------------
        1588             : !---New trace-gas rainout from large-scale precip with two time scales,
        1589             : !---one based on precip formation from cloud water and one based on
        1590             : !---Henry's Law solubility: correct limit for delta-t
        1591             : !---
        1592             : !---NB this code does not consider the aqueous dissociation (eg, C-q)
        1593             : !---   that makes uptake of HNO3 and H2SO4 so complete.  To do so would
        1594             : !---   require that we keep track of the pH of the falling rain.
        1595             : !---THUS the Henry's Law coefficient KHA needs to be enhanced to incldue this!
        1596             : !---ALSO the possible formation of other soluble species from, eg, CH2O, H2O2
        1597             : !---   can be considered with enhanced values of KHA.
        1598             : !---
        1599             : !---Does NOT now use RMC (moist conv rain) but could, assuming 30% coverage
        1600             : !-----------------------------------------------------------------------
        1601             :       implicit none
        1602             :       real(r8), intent(in) :: RRAIN       !new rain formation in box (kg/s)
        1603             :       real(r8), intent(in) :: DTSCAV      !time step (s)
        1604             :       real(r8), intent(in) :: CLWX,CFX !cloud water and cloud fraction
        1605             :       real(r8), intent(in) :: QM          !air mass in box (kg)
        1606             :       real(r8), intent(in) :: QT          !tracer in box (kg)
        1607             :       real(r8), intent(in) :: QTDIS          !tracer in aqueous phase (kg)
        1608             :       real(r8), intent(out) :: QTRAIN      !tracer picked up by new rain
        1609             :
        1610             :       real(r8)   QTLF,QTDISSTAR
        1611             :
        1612             :
        1613             :
        1614             :
        1615             :
        1616   804624160 :       QTDISSTAR=(QTDIS*(QT*CFX))/(QTDIS+(QT*CFX))
        1617             :
        1618             : !---Tracer Loss frequency (1/s) within cloud fraction:
        1619   804624160 :       QTLF = (RRAIN*QTDISSTAR)/(CLWX*QM*QT*CFX)
        1620             :
        1621             : !---in time = DTSCAV, the amount of QTT scavenged is calculated
        1622             : !---from CF*AMOUNT OF UPTAKE
        1623   804624160 :       QTRAIN = QT*CFX*(1._r8 - exp(-DTSCAV*QTLF))
        1624             :
        1625   804624160 :       return
        1626             :       end subroutine RAINGAS
        1627             :
        1628             :
        1629             : !-----------------------------------------------------------------------
        1630   935476727 :       subroutine WASHGAS (RWASH,BOXF,DTSCAV,QTRTOP,HSTAR,TM,PR,QM, &
        1631             :                             QT,QTWASH,QTEVAP)
        1632             : !-----------------------------------------------------------------------
        1633             : !---for most gases below-cloud washout assume Henry-Law equilib with precip
        1634             : !---assumes that precip is liquid, if frozen, do not call this sub
        1635             : !---since solubility is moderate, fraction of box with rain does not matter
        1636             : !---NB this code does not consider the aqueous dissociation (eg, C-q)
        1637             : !---   that makes uptake of HNO3 and H2SO4 so complete.  To do so would
        1638             : !---   require that we keep track of the pH of the falling rain.
        1639             : !---THUS the Henry's Law coefficient KHA needs to be enhanced to incldue this!
        1640             : !---ALSO the possible formation of other soluble species from, eg, CH2O, H2O2
        1641             : !---   can be considered with enhanced values of KHA.
        1642             : !-----------------------------------------------------------------------
        1643             :       implicit none
        1644             :       real(r8), intent(in)  :: RWASH   ! precip leaving bottom of box (kg/s)
        1645             :       real(r8), intent(in)  :: BOXF   ! fraction of box with washout
        1646             :       real(r8), intent(in)  :: DTSCAV  ! time step (s)
        1647             :       real(r8), intent(in)  :: QTRTOP  ! tracer-T in rain entering top of box
        1648             : !                                              over time step (kg)
        1649             :       real(r8), intent(in)  :: HSTAR ! Henry's Law coeffs A*exp(-B/T)
        1650             :       real(r8), intent(in)  :: TM      ! temperature of box (K)
        1651             :       real(r8), intent(in)  :: PR      ! pressure of box (hPa)
        1652             :       real(r8), intent(in)  :: QT      ! tracer in box (kg)
        1653             :       real(r8), intent(in)  :: QM      ! air mass in box (kg)
        1654             :       real(r8), intent(out) :: QTWASH  ! tracer picked up by precip (kg)
        1655             :       real(r8), intent(out) :: QTEVAP  ! tracer evaporated from precip (kg)
        1656             :
        1657             :       real(r8), parameter :: INV298 = 1._r8/298._r8
        1658             :       real(r8)            :: FWASH, QTMAX, QTDIF
        1659             :
        1660             : !---effective Henry's Law constant: H* = moles-T / liter-precip / press(atm-T)
        1661             : !---p(atm of tracer-T) = (QT/QM) * (.029/MolWt-T) * pressr(hPa)/1000
        1662             : !---limit temperature effects to T above freezing
        1663             :
        1664             : !
        1665             : ! jfl
        1666             : !
        1667             : ! added test for BOXF = 0.
        1668             : !
        1669   935476727 :       if ( BOXF == 0._r8 ) then
        1670       89402 :         QTWASH = 0._r8
        1671       89402 :         QTEVAP = 0._r8
        1672       89402 :         return
        1673             :       end if
        1674             :
        1675             : !---effective washout frequency (1/s):
        1676   935387325 :         FWASH = (RWASH*HSTAR*29.e-6_r8*PR)/(QM*BOXF)
        1677             : !---equilib amount of T (kg) in rain thru bottom of box over time step
        1678   935387325 :         QTMAX = QT*FWASH*DTSCAV
        1679   935387325 :       if (QTMAX .gt. QTRTOP) then
        1680             : !---more of tracer T can go into rain
        1681   758406523 :          QTDIF = min (QT, QTMAX-QTRTOP)
        1682   758406523 :          QTWASH = QTDIF * (1._r8 - exp(-DTSCAV*FWASH))
        1683   758406523 :          QTEVAP=0._r8
        1684             :       else
        1685             : !--too much of T in rain, must degas/evap T
        1686   176980802 :          QTWASH = 0._r8
        1687   176980802 :          QTEVAP = QTRTOP - QTMAX
        1688             :       endif
        1689             :
        1690             :       return
        1691             :       end subroutine WASHGAS
        1692             :
        1693             : !-----------------------------------------------------------------------
        1694  2065968825 :       function DEMPIRICAL (CWATER,RRATE)
        1695             : !-----------------------------------------------------------------------
        1696             :       use shr_spfn_mod, only: shr_spfn_gamma
        1697             :
        1698             :       implicit none
        1699             :       real(r8), intent(in)  :: CWATER
        1700             :       real(r8), intent(in)  :: RRATE
        1701             :
        1702             :       real(r8) :: DEMPIRICAL
        1703             :
        1704             :       real(r8) RRATEX,WX,THETA,PHI,ETA,BETA,ALPHA,BEE
        1705             :       real(r8) GAMTHETA,GAMBETA
        1706             :
        1707             :
        1708             :
        1709  2065968825 :       RRATEX=RRATE*3600._r8       !mm/hr
        1710  2065968825 :       WX=CWATER*1.0e3_r8  !g/m3
        1711             :
        1712  2065968825 :       if(RRATEX .gt. 0.04_r8) then
        1713   201216160 :          THETA=exp(-1.43_r8*dlog10(7._r8*RRATEX))+2.8_r8
        1714             :       else
        1715  1864752665 :          THETA=5._r8
        1716             :       endif
        1717  2065968825 :       PHI=RRATEX/(3600._r8*10._r8) !cgs units
        1718  2065968825 :       ETA=exp((3.01_r8*THETA)-10.5_r8)
        1719  2065968825 :       BETA=THETA/(1._r8+0.638_r8)
        1720  2065968825 :       ALPHA=exp(4._r8*(BETA-3.5_r8))
        1721  2065968825 :       BEE=(.638_r8*THETA/(1._r8+.638_r8))-1.0_r8
        1722  2065968825 :       GAMTHETA = shr_spfn_gamma(THETA)
        1723  2065968825 :       GAMBETA  = shr_spfn_gamma(BETA+1._r8)
        1724             :       DEMPIRICAL=(((WX*ETA*GAMTHETA)/(1.0e6_r8*ALPHA*PHI*GAMBETA))** &
        1725  2065968825 :                  (-1._r8/BEE))*10._r8      ! in mm (wx/1e6 for cgs)
        1726             :
        1727             :
        1728             :       return
        1729             :       end function DEMPIRICAL
        1730             : !
        1731             : end module mo_neu_wetdep

Generated by: LCOV version 1.14