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