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
|
