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 188928 : mapping_to_heff = -99
81 188928 : do m=1,gas_wetdep_cnt
82 : !
83 187392 : test_name = gas_wetdep_list(m)
84 187392 : 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 187392 : if ( .not. cam_chempkg_is('geoschem_mam4') ) then
92 374784 : 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 13824 : test_name = 'HNO3'
102 : case ( 'NH_50W', 'NDEP', 'NHDEP', 'NH4', 'NH4NO3' )
103 4608 : test_name = 'HNO3'
104 : case( 'SOAGbb0' ) ! Henry's Law coeff. added for VBS SOA's, biomass burning is the same as fossil fuels
105 1536 : test_name = 'SOAGff0'
106 : case( 'SOAGbb1' )
107 1536 : test_name = 'SOAGff1'
108 : case( 'SOAGbb2' )
109 1536 : test_name = 'SOAGff2'
110 : case( 'SOAGbb3' )
111 1536 : test_name = 'SOAGff3'
112 : case( 'SOAGbb4' )
113 374784 : test_name = 'SOAGff4'
114 : end select
115 : endif
116 : !
117 12678144 : do l = 1,n_species_table
118 : !
119 : ! if ( debug ) print '(i4,a)',l,trim(species_name_table(l))
120 : !
121 12678144 : if( trim(test_name) == trim( species_name_table(l) ) ) then
122 187392 : mapping_to_heff(m) = l
123 187392 : if ( debug ) print '(a,a,i4)','mapping to heff of ',trim(species_name_table(l)),l
124 187392 : exit
125 : end if
126 : end do
127 187392 : 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 187392 : if ( trim(test_name) == 'NH3' ) then
135 1536 : nh3_ndx = m
136 : end if
137 187392 : if ( trim(test_name) == 'CO2' ) then
138 0 : co2_ndx = m
139 : end if
140 187392 : if ( trim(gas_wetdep_list(m)) == 'HNO3' ) then
141 1536 : hno3_ndx = m
142 : end if
143 187392 : if ( trim(test_name) == 'SO2' ) then
144 1536 : so2_ndx = m
145 : end if
146 187392 : if ( trim(test_name) == 'SO4' ) then ! GEOS-Chem bulk sulfate
147 0 : so4_ndx = m
148 : end if
149 188928 : 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 188928 : 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 188928 : mapping_to_mmr = -99
166 188928 : do m=1,gas_wetdep_cnt
167 187392 : if ( debug .and. masterproc ) write(iulog, '(i4,a)') m,trim(gas_wetdep_list(m))
168 187392 : call cnst_get_ind(gas_wetdep_list(m), mapping_to_mmr(m), abort=.false. )
169 187392 : if ( debug .and. masterproc ) write(iulog, '(a,i4)') 'mapping_to_mmr ',mapping_to_mmr(m)
170 188928 : 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 188928 : do m=1,gas_wetdep_cnt
179 : !
180 187392 : mol_weight(m) = cnst_mw(mapping_to_mmr(m))
181 187392 : if ( debug .and. masterproc ) write(iulog, '(i4,a,f8.4)') m,' mol_weight ',mol_weight(m)
182 187392 : ice_uptake(m) = .false.
183 188928 : if ( trim(gas_wetdep_list(m)) == 'HNO3' ) then
184 1536 : 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 188928 : do m=1,gas_wetdep_cnt
198 374784 : call addfld ('DTWR_'//trim(gas_wetdep_list(m)),(/ 'lev' /), 'A','kg/kg/s','wet removal Neu scheme tendency')
199 187392 : call addfld ('WD_'//trim(gas_wetdep_list(m)),horiz_only, 'A','kg/m2/s','vertical integrated wet deposition flux')
200 374784 : call addfld ('HEFF_'//trim(gas_wetdep_list(m)),(/ 'lev' /), 'A','M/atm','Effective Henrys Law coeff.')
201 188928 : if (history_chemistry) then
202 0 : 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 72960 : subroutine neu_wetdep_tend(lchnk,ncol,mmr,pmid,pdel,zint,tfld,delt, &
224 72960 : 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 145920 : real(r8), dimension(ncol) :: area, wk_out
255 145920 : real(r8), dimension(ncol,pver) :: cldice,cldliq,cldfrc,totprec,totevap,delz,p
256 145920 : real(r8), dimension(ncol,pver) :: rls,evaprate,mass_in_layer,temp
257 145920 : real(r8), dimension(ncol,pver,gas_wetdep_cnt) :: trc_mass,heff,dtwr
258 145920 : real(r8), dimension(ncol,pver,gas_wetdep_cnt) :: wd_mmr
259 145920 : logical , dimension(gas_wetdep_cnt) :: tckaqb
260 145920 : integer , dimension(ncol) :: test_flag
261 : !
262 : ! arrays for HNO3 diagnostics
263 : !
264 145920 : 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 145920 : 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 72960 : if (.not.do_neu_wetdep) return
282 : !
283 : ! don't do anything if there are no species to be removed
284 : !
285 72960 : if ( gas_wetdep_cnt == 0 ) return
286 : !
287 : ! reset output variables
288 : !
289 72960 : wd_tend_int = 0._r8
290 : !
291 : ! get area (in radians square)
292 : !
293 72960 : call get_area_all_p(lchnk, ncol, area)
294 1123584 : 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 2407680 : do k=1,pver
300 2334720 : kk = pver - k + 1
301 36027648 : do i=1,ncol
302 : !
303 33619968 : mass_in_layer(i,k) = area(i) * pdel(i,kk)/gravit ! kg
304 : !
305 33619968 : cldice (i,k) = mmr(i,kk,index_cldice) ! kg/kg
306 33619968 : cldliq (i,k) = mmr(i,kk,index_cldliq) ! kg/kg
307 33619968 : cldfrc (i,k) = cld(i,kk) ! unitless
308 : !
309 : totprec(i,k) = (prain(i,kk)+cmfdqr(i,kk)) &
310 33619968 : * mass_in_layer(i,k) ! kg/s
311 33619968 : totevap(i,k) = nevapr(i,kk) * mass_in_layer(i,k) ! kg/s
312 : !
313 33619968 : delz(i,k) = zint(i,kk) - zint(i,kk+1) ! in m
314 : !
315 33619968 : temp(i,k) = tfld(i,kk)
316 : !
317 : ! convert tracer mass to kg to kg/kg
318 : !
319 4135256064 : trc_mass(i,k,:) = mmr(i,kk,mapping_to_mmr(:)) * mass_in_layer(i,k)
320 : !
321 35954688 : 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 4395446016 : 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 1123584 : rls (:,pver) = 0._r8
334 1123584 : evaprate (:,pver) = 0._r8
335 2334720 : do k=pver-1,1,-1
336 34831104 : 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 34904064 : 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 4395446016 : heff = 0._r8
344 2407680 : do k=1,pver
345 : !
346 2334720 : kk = pver - k + 1
347 : !
348 35954688 : wrk(:) = (t0-tfld(1:ncol,kk))/(t0*tfld(1:ncol,kk))
349 : !
350 287243520 : do m=1,gas_wetdep_cnt
351 : !
352 284835840 : l = mapping_to_heff(m)
353 284835840 : e298 = dheff(1,l)
354 284835840 : dhr = dheff(2,l)
355 4386471936 : heff(:,k,m) = e298*exp( dhr*wrk(:) )
356 4386471936 : test_flag = -99
357 284835840 : if( dheff(3,l) /= 0._r8 .and. dheff(5,l) == 0._r8 ) then
358 39690240 : e298 = dheff(3,l)
359 39690240 : dhr = dheff(4,l)
360 611229696 : dk1s(:) = e298*exp( dhr*wrk(:) )
361 2976768000 : where( heff(:,k,m) /= 0._r8 )
362 39690240 : heff(:,k,m) = heff(:,k,m)*(1._r8 + dk1s(:)*ph_inv)
363 : elsewhere
364 : test_flag = 1
365 39690240 : heff(:,k,m) = dk1s(:)*ph_inv
366 : endwhere
367 : end if
368 : !
369 4386471936 : 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 287170560 : if( dheff(5,l) /= 0._r8 ) then
374 4669440 : if( nh3_ndx > 0 .or. co2_ndx > 0 .or. so2_ndx > 0 .or. so4_ndx > 0 .or. so4s_ndx > 0 ) then
375 4669440 : e298 = dheff(3,l)
376 4669440 : dhr = dheff(4,l)
377 71909376 : dk1s(:) = e298*exp( dhr*wrk(:) )
378 4669440 : e298 = dheff(5,l)
379 4669440 : dhr = dheff(6,l)
380 71909376 : dk2s(:) = e298*exp( dhr*wrk(:) )
381 4669440 : if( m == co2_ndx .or. m == so2_ndx .or. m == so4_ndx .or. m == so4s_ndx ) then
382 35954688 : heff(:,k,m) = heff(:,k,m)*(1._r8 + dk1s(:)*ph_inv*(1._r8 + dk2s(:)*ph_inv))
383 2334720 : else if( m == nh3_ndx ) then
384 35954688 : 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 72960 : 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 1123584 : do i=1,ncol
407 : !
408 4202496 : call washo(pver,gas_wetdep_cnt,delt,trc_mass(i,:,:),mass_in_layer(i,:),p(i,:),delz(i,:) &
409 6303744 : ,rls(i,:),cldliq(i,:),cldice(i,:),cldfrc(i,:),temp(i,:),evaprate(i,:) &
410 1050624 : ,area(i),heff(i,:,:),mol_weight(:),tckaqb(:),ice_uptake(:) &
411 13273656576 : ,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 2407680 : do k=1,pver
419 2334720 : kk = pver - k + 1
420 36027648 : do i=1,ncol
421 : !
422 : ! convert tracer mass from kg
423 : !
424 4137590784 : 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 4395446016 : dtwr(1:ncol,:,:) = wd_mmr(1:ncol,:,:) - dtwr(1:ncol,:,:)
432 4395446016 : dtwr(1:ncol,:,:) = dtwr(1:ncol,:,:) / delt
433 :
434 : ! polarward of 60S, 60N and <200hPa set to zero!
435 72960 : call get_rlat_all_p(lchnk, pcols, lats )
436 2407680 : do k = 1, pver
437 36027648 : do i= 1, ncol
438 35954688 : if ( abs( lats(i)*rad2deg ) > 60._r8 ) then
439 11206656 : if ( pmid(i,k) < 20000._r8) then
440 646133760 : dtwr(i,k,:) = 0._r8
441 : endif
442 : endif
443 : end do
444 : end do
445 : !
446 : ! output tendencies
447 : !
448 8974080 : do m=1,gas_wetdep_cnt
449 4395373056 : wd_tend(1:ncol,:,mapping_to_mmr(m)) = wd_tend(1:ncol,:,mapping_to_mmr(m)) + dtwr(1:ncol,:,m)
450 8901120 : call outfld( 'DTWR_'//trim(gas_wetdep_list(m)),dtwr(:,:,m),ncol,lchnk )
451 :
452 4395373056 : 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 137077248 : wk_out = 0._r8
457 293736960 : do k=1,pver
458 284835840 : kk = pver - k + 1
459 4395373056 : 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 8901120 : 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 8901120 : if ( debug .and. masterproc ) then
466 0 : write(iulog,*) 'mo_neu ',mapping_to_mmr(m),(wk_out(1:ncol))
467 : endif
468 137150208 : wd_tend_int(1:ncol,mapping_to_mmr(m)) = wk_out(1:ncol)
469 : !
470 : end do
471 : !
472 72960 : 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 72960 : 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 1050624 : subroutine WASHO(LPAR,NTRACE,DTSCAV,QTTJFL,QM,POFL,DELZ, &
490 1050624 : RLS,CLWC,CIWC,CFR,TEM,EVAPRATE,GAREA,HSTAR,TCMASS,TCKAQB, &
491 1050624 : 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 2101248 : real(r8), dimension(LPAR) :: CFXX
520 2101248 : 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 2101248 : real(r8) :: rls_wrk(lpar)
561 2101248 : real(r8) :: rnew_wrk(lpar)
562 2101248 : real(r8) :: rca_wrk(lpar)
563 2101248 : real(r8) :: fca_wrk(lpar)
564 2101248 : real(r8) :: rcxa_wrk(lpar)
565 2101248 : real(r8) :: fcxa_wrk(lpar)
566 2101248 : real(r8) :: rcxb_wrk(lpar)
567 2101248 : real(r8) :: fcxb_wrk(lpar)
568 2101248 : real(r8) :: rax_wrk(lpar,2)
569 2101248 : real(r8) :: fax_wrk(lpar,2)
570 2101248 : real(r8) :: rama_wrk(lpar)
571 2101248 : real(r8) :: fama_wrk(lpar)
572 2101248 : real(r8) :: deltarime_wrk(lpar)
573 2101248 : real(r8) :: clwx_wrk(lpar)
574 2101248 : real(r8) :: frc(lpar,3)
575 2101248 : real(r8) :: rlsog(lpar)
576 : !
577 : logical :: is_hno3
578 2101248 : logical :: rls_flag(lpar)
579 2101248 : logical :: rnew_flag(lpar)
580 2101248 : logical :: cf_trigger(lpar)
581 2101248 : 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 1050624 : 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 1050624 : LE = LPAR-1
599 : !
600 33619968 : rls_flag(1:le) = rls(1:le) > zero
601 33619968 : freezing(1:le) = tem(1:le) < tice
602 33619968 : rlsog(1:le) = rls(1:le)/garea
603 : !
604 : species_loop : &
605 129226752 : do N = 1,NTRACE
606 4229812224 : QTT(:lpar) = QTTJFL(:lpar,N)
607 4229812224 : QTTNEW(:lpar) = QTTJFL(:lpar,N)
608 128176128 : is_hno3 = n == hno3_ndx
609 128176128 : if( is_hno3 ) then
610 34670592 : qt_rain(:lpar) = zero
611 34670592 : qt_rime(:lpar) = zero
612 34670592 : qt_wash(:lpar) = zero
613 34670592 : qt_evap(:lpar) = zero
614 34670592 : rca_wrk(:lpar) = zero
615 34670592 : fca_wrk(:lpar) = zero
616 34670592 : rcxa_wrk(:lpar) = zero
617 34670592 : fcxa_wrk(:lpar) = zero
618 34670592 : rcxb_wrk(:lpar) = zero
619 34670592 : fcxb_wrk(:lpar) = zero
620 34670592 : rls_wrk(:lpar) = zero
621 34670592 : rnew_wrk(:lpar) = zero
622 34670592 : cf_trigger(:lpar) = .false.
623 34670592 : clwx_wrk(:lpar) = -9999._r8
624 34670592 : deltarime_wrk(:lpar) = -9999._r8
625 70391808 : rax_wrk(:lpar,:) = zero
626 70391808 : fax_wrk(:lpar,:) = zero
627 : endif
628 :
629 : !-----------------------------------------------------------------------
630 : ! check whether soluble in ice
631 : !-----------------------------------------------------------------------
632 128176128 : if( TCNION(N) ) then
633 : LICETYP = 1
634 : else
635 127125504 : LICETYP = 2
636 : end if
637 :
638 : !-----------------------------------------------------------------------
639 : ! initialization
640 : !-----------------------------------------------------------------------
641 128176128 : QTTOPAA = zero
642 128176128 : QTTOPCA = zero
643 :
644 128176128 : RCA = zero
645 128176128 : FCA = zero
646 128176128 : DCA = zero
647 128176128 : RAMA = zero
648 128176128 : FAMA = zero
649 128176128 : DAMA = zero
650 :
651 128176128 : AMPCT = zero
652 128176128 : AMCLPCT = zero
653 128176128 : CLNEWPCT = zero
654 128176128 : CLNEWAMPCT = zero
655 128176128 : CLOLDPCT = zero
656 128176128 : CLOLDAMPCT = zero
657 : !-----------------------------------------------------------------------
658 : ! Check whether precip in top layer - if so, require CF ge 0.2
659 : !-----------------------------------------------------------------------
660 128176128 : if( RLS(LE) > zero ) then
661 2936418 : CFXX(LE) = max( CFMIN,CFR(LE) )
662 : else
663 125239710 : CFXX(LE) = CFR(LE)
664 : endif
665 :
666 4101636096 : rnew_flag(1:le) = .false.
667 :
668 : level_loop : &
669 4101636096 : do L = LE,1,-1
670 3973459968 : LM1 = L - 1
671 3973459968 : FAX = zero
672 3973459968 : RAX = zero
673 3973459968 : DAX = zero
674 3973459968 : FCXA = zero
675 3973459968 : FCXB = zero
676 3973459968 : DCXA = zero
677 3973459968 : DCXB = zero
678 3973459968 : RCXA = zero
679 3973459968 : RCXB = zero
680 :
681 3973459968 : QTDISCF = zero
682 3973459968 : QTDISRIME = zero
683 3973459968 : QTDISCXA = zero
684 :
685 3973459968 : QTEVAPAXP = zero
686 3973459968 : QTEVAPAXW = zero
687 3973459968 : QTEVAPAX = zero
688 3973459968 : QTWASHAX = zero
689 :
690 3973459968 : QTEVAPCXAP = zero
691 3973459968 : QTEVAPCXAW = zero
692 3973459968 : QTEVAPCXA = zero
693 3973459968 : QTRIMECXA = zero
694 3973459968 : QTWASHCXA = zero
695 3973459968 : QTRAINCXA = zero
696 3973459968 : QTRAINCXB = zero
697 :
698 3973459968 : RAMPCT = zero
699 :
700 3973459968 : RPRECIP = zero
701 3973459968 : DELTARIMEMASS = zero
702 3973459968 : DELTARIME = zero
703 3973459968 : DOR = zero
704 3973459968 : DNEW = zero
705 :
706 3973459968 : QTTOPAAX = zero
707 3973459968 : QTTOPCAX = zero
708 :
709 : has_rls : &
710 3973459968 : 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 1603382682 : FAX = max( zero,FAMA*(one - evaprate(l)) )
720 1603382682 : RAX = RAMA !kg/m2/s
721 1603382682 : 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 1603382682 : if( FAMA > zero ) then
727 741612868 : if( freezing(l) ) then
728 289291036 : DAX = DAMA !mm
729 : else
730 452321832 : DAX = four !mm - not necessary
731 : endif
732 : else
733 861769814 : DAX = zero
734 : endif
735 :
736 1603382682 : if( RAMA > zero ) then
737 741612868 : QTEVAPAXP = min( QTTOPAA,EVAPRATE(L)*QTTOPAA )
738 : else
739 : QTEVAPAXP = zero
740 : endif
741 1603382682 : if( is_hno3 ) then
742 13142481 : rax_wrk(l,1) = rax
743 13142481 : fax_wrk(l,1) = fax
744 : endif
745 :
746 :
747 : !-----------------------------------------------------------------------
748 : ! Determine how much the in-cloud precip rate has increased------
749 : !-----------------------------------------------------------------------
750 1603382682 : WRK = RAX*FAX + RCA*FCA
751 1603382682 : if( WRK > 0._r8 ) then
752 1431329254 : RNEW_TST = RLS(L)/(GAREA * WRK)
753 : else
754 : RNEW_TST = 10._r8
755 : endif
756 1603382682 : RNEW = RLSOG(L) - (RAX*FAX + RCA*FCA) !GBA*CF
757 1603382682 : rnew_wrk(l) = rnew_tst
758 1603382682 : 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 1603382682 : 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 1172617884 : 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 1172617884 : rnew_flag(l) = .true.
802 1172617884 : CLWX = max( CLWC(L)+CIWC(L),CWMIN*CFXX(L) )
803 1172617884 : if( is_hno3 ) then
804 9611622 : clwx_wrk(l) = clwx
805 : endif
806 : !-----------------------------------------------------------------------
807 : ! Area of old cloud and new cloud
808 : !-----------------------------------------------------------------------
809 1172617884 : FCXA = FCA
810 1172617884 : 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 1172617884 : if( freezing(l) ) then
820 717829456 : COLEFFSNOW = exp( 2.5e-2_r8*(TEM(L) - TICE) )
821 717829456 : if( TEM(L) <= TFROZ ) then
822 : RHOSNOW = RHOSNOWFIX
823 : else
824 330624880 : RHOSNOW = 0.303_r8*(TEM(L) - TFROZ)*RHOSNOWFIX
825 : endif
826 717829456 : if( FCXA > zero ) then
827 586001624 : if( DCA > zero ) then
828 : DELTARIMEMASS = CLWX*QM(L)*(FCXA/CFXX(L))* &
829 586001624 : (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 717829456 : if( FCXA > zero ) then
841 586001624 : DELTARIME = min( RNEW/FCXA,DELTARIMEMASS/(FCXA*GAREA*DTSCAV) ) !GBA
842 : else
843 131827832 : DELTARIME = zero
844 : endif
845 717829456 : if( is_hno3 ) then
846 5883848 : deltarime_wrk(l) = deltarime
847 : endif
848 : !-----------------------------------------------------------------------
849 : ! Find diameter of rimed precip, must be at least .1mm
850 : !-----------------------------------------------------------------------
851 717829456 : if( RCA > zero ) then
852 586001624 : DOR = max( DMIN,(((RCA+DELTARIME)/RCA)**VOLPOW)*DCA )
853 : else
854 131827832 : 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 717829456 : RPRECIP = (RNEW-(DELTARIME*FCXA))/CFXX(L) !kg/m2/s !GBA
863 : !-----------------------------------------------------------------------
864 : ! Calculate precip rate in old and new cloud fractions
865 : !-----------------------------------------------------------------------
866 717829456 : RCXA = RCA + DELTARIME + RPRECIP !kg/m2/s GBA
867 717829456 : 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 717829456 : if( RPRECIP > zero ) then
879 632965524 : WEMP = (CLWX*QM(L))/(GAREA*CFXX(L)*DELZ(L)) !kg/m3
880 632965524 : REMP = RPRECIP/((RHORAIN/1.e3_r8)) !mm/s local
881 632965524 : DNEW = DEMPIRICAL( WEMP, REMP )
882 632965524 : 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 632965524 : DNEW = max( DMIN,DNEW )
891 632965524 : if( FCXB > zero ) then
892 249436442 : DCXB = DNEW
893 : else
894 383529082 : DCXB = zero
895 : endif
896 : else
897 84863932 : DCXB = zero
898 : endif
899 :
900 717829456 : if( FCXA > zero ) then
901 586001624 : WEMP = (CLWX*QM(L)*(FCXA/CFXX(L)))/(GAREA*FCXA*DELZ(L)) !kg/m3
902 586001624 : REMP = RCXA/((RHORAIN/1.e3_r8)) !mm/s local
903 586001624 : DEMP = DEMPIRICAL( WEMP, REMP )
904 586001624 : DCXA = ((RCA+DELTARIME)/RCXA)*DOR + (RPRECIP/RCXA)*DNEW
905 586001624 : DCXA = max( DEMP,DCXA )
906 586001624 : DCXA = max( DMIN,DCXA )
907 : else
908 131827832 : DCXA = zero
909 : endif
910 717829456 : 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 717829456 : 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 717829456 : if( RPRECIP > zero ) then
931 632965524 : if( LICETYP == 1 ) then
932 5188242 : 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 5188242 : QTT(L)*CFXX(L),QTDISCF )
936 5188242 : call RAINGAS( RRAIN, DTSCAV, CLWX, CFXX(L), &
937 10376484 : QM(L), QTT(L), QTDISCF, QTRAIN )
938 5188242 : WRK = QTRAIN/CFXX(L)
939 5188242 : QTRAINCXA = FCXA*WRK
940 5188242 : QTRAINCXB = FCXB*WRK
941 627777282 : elseif( LICETYP == 2 ) then
942 627777282 : QTRAINCXA = zero
943 627777282 : QTRAINCXB = zero
944 : endif
945 632965524 : 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 717829456 : if( DELTARIME > zero ) then
961 550307962 : if( LICETYP == 1 ) then
962 4510721 : if( TEM(L) <= TFROZ ) then
963 : RHOSNOW = RHOSNOWFIX
964 : else
965 2487297 : RHOSNOW = 0.303_r8*(TEM(L) - TFROZ)*RHOSNOWFIX
966 : endif
967 4510721 : QTCXA = QTT(L)*FCXA
968 : call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N), &
969 : HSTAR(L,N), TEM(L), POFL(L), &
970 4510721 : QM(L), QTCXA, QTDISRIME )
971 4510721 : QTDISSTAR = (QTDISRIME*QTCXA)/(QTDISRIME + QTCXA)
972 4510721 : 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 4510721 : (QTDISSTAR/QTCXA)*DTSCAV))
984 : QTRIMECXA = min( QTRIMECXA, &
985 4510721 : ((RNEW*GAREA*DTSCAV)/(CLWX*QM(L)*(FCXA/CFXX(L))))*QTDISSTAR)
986 545797241 : elseif( LICETYP == 2 ) then
987 545797241 : 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 717829456 : QTWASHCXA = zero
999 717829456 : 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 454788428 : if( FCXA > zero ) then
1008 : DELTARIMEMASS = (CLWX*QM(L))*(FCXA/CFXX(L))* &
1009 414562832 : (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 454788428 : if( FCXA > zero ) then
1018 414562832 : DELTARIME = min( RNEW/FCXA,DELTARIMEMASS/(FCXA*GAREA*DTSCAV) ) !GBA
1019 : else
1020 40225596 : DELTARIME = zero
1021 : endif
1022 : !-----------------------------------------------------------------------
1023 : ! If there is some in-cloud precip left, we have new precip formation
1024 : !-----------------------------------------------------------------------
1025 454788428 : RPRECIP = (RNEW-(DELTARIME*FCXA))/CFXX(L) !GBA
1026 :
1027 454788428 : RCXA = RCA + DELTARIME + RPRECIP !kg/m2/s GBA
1028 454788428 : RCXB = RPRECIP !kg/m2/s GBA
1029 454788428 : DCXA = FOUR
1030 454788428 : if( FCXB > zero ) then
1031 113679844 : DCXB = FOUR
1032 : else
1033 341108584 : DCXB = zero
1034 : endif
1035 : !-----------------------------------------------------------------------
1036 : ! RAIN SCAVENGING
1037 : ! For rain, rainout both hno3/aerosols and gases using new precip
1038 : !-----------------------------------------------------------------------
1039 454788428 : if( QTT(L) > zero ) then
1040 454788428 : if( RPRECIP > zero ) then
1041 385610402 : 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 385610402 : QTT(L)*CFXX(L), QTDISCF )
1045 385610402 : call RAINGAS( RRAIN, DTSCAV, CLWX, CFXX(L), &
1046 771220804 : QM(L), QTT(L), QTDISCF, QTRAIN )
1047 385610402 : WRK = QTRAIN/CFXX(L)
1048 385610402 : QTRAINCXA = FCXA*WRK
1049 385610402 : QTRAINCXB = FCXB*WRK
1050 385610402 : 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 454788428 : if( DELTARIME > zero ) then
1065 414562100 : QTCXA = QTT(L)*FCXA
1066 : call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N), &
1067 : HSTAR(L,N), TEM(L), POFL(L), &
1068 414562100 : QM(L), QTCXA, QTDISRIME )
1069 414562100 : QTDISSTAR = (QTDISRIME*QTCXA)/(QTDISRIME + QTCXA)
1070 : QTRIMECXA = QTCXA* &
1071 : (one - exp(-0.24_r8*COLEFFRAIN* &
1072 : ((RCA)**0.75_r8)* & !local
1073 414562100 : (QTDISSTAR/QTCXA)*DTSCAV))
1074 : QTRIMECXA = min( QTRIMECXA, &
1075 414562100 : ((RNEW*GAREA*DTSCAV)/(CLWX*QM(L)*(FCXA/CFXX(L))))*QTDISSTAR)
1076 : else
1077 40226328 : 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 454788428 : if( RCA > zero ) then
1090 414562832 : QTPRECIP = FCXA*QTT(L) - QTDISRIME
1091 414562832 : if( HSTAR(L,N) > 1.e4_r8 ) then
1092 225138859 : if( QTPRECIP > zero ) then
1093 108771753 : QTWASHCXA = QTPRECIP*(one - exp( -0.24_r8*COLEFFAER*((RCA)**0.75_r8)*DTSCAV )) !local
1094 : else
1095 116367106 : QTWASHCXA = zero
1096 : endif
1097 225138859 : QTEVAPCXA = zero
1098 : else
1099 189423973 : RWASH = RCA*GAREA !kg/s local
1100 189423973 : if( QTPRECIP > zero ) then
1101 : call WASHGAS( RWASH, FCA, DTSCAV, QTTOPCA+QTRIMECXA, &
1102 : HSTAR(L,N), TEM(L), POFL(L), &
1103 189423973 : 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 430764798 : CLWX = CLWC(L) + CIWC(L)
1118 430764798 : if( is_hno3 ) then
1119 3530859 : clwx_wrk(l) = clwx
1120 : endif
1121 430764798 : FCXA = FCA
1122 430764798 : FCXB = max( zero,CFXX(L)-FCXA )
1123 430764798 : RCXB = zero
1124 430764798 : DCXB = zero
1125 430764798 : QTRAINCXA = zero
1126 430764798 : QTRAINCXB = zero
1127 430764798 : 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 430764798 : if( FCXA > zero ) then
1136 400512824 : RCXA = min( RCA,RLS(L)/(GAREA*FCXA) ) !kg/m2/s GBA
1137 400512824 : if( FAX > zero .and. ((RCXA+1.e-12_r8) < RLS(L)/(GAREA*FCXA)) ) then
1138 268162954 : RAXADJF = RLS(L)/GAREA - RCXA*FCXA
1139 268162954 : RAMPCT = RAXADJF/(RAX*FAX)
1140 268162954 : FAXADJ = RAMPCT*FAX
1141 268162954 : if( FAXADJ > zero ) then
1142 268162954 : 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 30251974 : RCXA = zero
1153 30251974 : if( FAX > zero ) then
1154 30251974 : RAXADJF = RLS(L)/GAREA
1155 30251974 : RAMPCT = RAXADJF/(RAX*FAX)
1156 30251974 : FAXADJ = RAMPCT*FAX
1157 30251974 : if( FAXADJ > zero ) then
1158 30251974 : 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 430764798 : QTEVAPAXP = min( QTTOPAA,QTTOPAA - (RAMPCT*(QTTOPAA-QTEVAPAXP)) )
1170 430764798 : FAX = FAXADJ
1171 430764798 : RAX = RAXADJ
1172 430764798 : 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 430764798 : if( RCXA <= zero ) then
1184 30251974 : QTEVAPCXA = QTTOPCA
1185 30251974 : RCXA = zero
1186 30251974 : 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 400512824 : if( freezing(l) ) then
1205 164369624 : QTWASHCXA = zero
1206 164369624 : DCXA = ((RCXA/RCA)**VOLPOW)*DCA
1207 164369624 : if( LICETYP == 1 ) then
1208 1347292 : if( TEM(L) <= TMIX ) then
1209 1157445 : 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 1157445 : QM(L), QTT(L), QTEVAPCXA )
1216 1157445 : QTEVAPCXA = min( QTTOPCA,QTEVAPCXA )
1217 : else
1218 189847 : QTEVAPCXA = zero
1219 : endif
1220 163022332 : elseif( LICETYP == 2 ) then
1221 163022332 : QTEVAPCXA = zero
1222 : endif
1223 : else is_freezing_a
1224 236143200 : QTEVAPCXAP = (RCA - RCXA)/RCA*QTTOPCA
1225 236143200 : DCXA = FOUR
1226 236143200 : QTCXA = FCXA*QTT(L)
1227 236143200 : if( HSTAR(L,N) > 1.e4_r8 ) then
1228 126266460 : if( QTT(L) > zero ) then
1229 : call DISGAS( CLWX*(FCXA/CFXX(L)), FCXA, TCMASS(N), &
1230 : HSTAR(L,N), TEM(L), POFL(L), &
1231 126266460 : QM(L), QTCXA, QTDISCXA )
1232 126266460 : if( QTCXA > QTDISCXA ) then
1233 116656047 : QTWASHCXA = (QTCXA - QTDISCXA)*(one - exp( -0.24_r8*COLEFFAER*((RCXA)**0.75_r8)*DTSCAV )) !local
1234 : else
1235 9610413 : QTWASHCXA = zero
1236 : endif
1237 126266460 : QTEVAPCXAW = zero
1238 : else
1239 0 : QTWASHCXA = zero
1240 0 : QTEVAPCXAW = zero
1241 : endif
1242 : else
1243 109876740 : RWASH = RCXA*GAREA !kg/s local
1244 : call WASHGAS( RWASH, FCXA, DTSCAV, QTTOPCA, HSTAR(L,N), &
1245 : TEM(L), POFL(L), QM(L), &
1246 109876740 : QTCXA-QTDISCXA, QTWASHCXA, QTEVAPCXAW )
1247 : endif
1248 236143200 : 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 1603382682 : if( RAX > zero ) then
1259 715549398 : if( .not. freezing(l) ) then
1260 442555122 : QTAX = FAX*QTT(L)
1261 442555122 : if( HSTAR(L,N) > 1.e4_r8 ) then
1262 : QTWASHAX = QTAX* &
1263 : (one - exp(-0.24_r8*COLEFFAER* &
1264 238070534 : ((RAX)**0.75_r8)*DTSCAV)) !local
1265 238070534 : QTEVAPAXW = zero
1266 : else
1267 204484588 : 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 204484588 : QTWASHAX, QTEVAPAXW )
1271 : endif
1272 : else
1273 272994276 : QTEVAPAXW = zero
1274 272994276 : QTWASHAX = zero
1275 : endif
1276 : else
1277 887833284 : QTEVAPAXW = zero
1278 887833284 : QTWASHAX = zero
1279 : endif
1280 1603382682 : 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 1603382682 : 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 1603382682 : if( L > 1 ) then
1300 1496826296 : if( CFR(LM1) >= CFMIN ) then
1301 738091338 : CFXX(LM1) = CFR(LM1)
1302 : else
1303 758734958 : if( adj_factor*RLSOG(LM1) >= (RCXA*FCXA + RCXB*FCXB + RAX*FAX)*(one - EVAPRATE(LM1)) ) then
1304 695210778 : CFXX(LM1) = CFMIN
1305 695210778 : cf_trigger(lm1) = .true.
1306 : else
1307 63524180 : CFXX(LM1) = CFR(LM1)
1308 : endif
1309 758734958 : 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 1496826296 : if( FAX > zero ) then
1321 646090772 : AMPCT = max( zero,min( one,(CFXX(L) + FAX - CFXX(LM1))/FAX ) )
1322 646090772 : AMCLPCT = one - AMPCT
1323 : else
1324 850735524 : AMPCT = zero
1325 850735524 : AMCLPCT = zero
1326 : endif
1327 1496826296 : if( FCXB > zero ) then
1328 371303340 : CLNEWPCT = max( zero,min( (CFXX(LM1) - FCXA)/FCXB,one ) )
1329 371303340 : CLNEWAMPCT = one - CLNEWPCT
1330 : else
1331 1125522956 : CLNEWPCT = zero
1332 1125522956 : CLNEWAMPCT = zero
1333 : endif
1334 1496826296 : if( FCXA > zero ) then
1335 1297295906 : CLOLDPCT = max( zero,min( CFXX(LM1)/FCXA,one ) )
1336 1297295906 : CLOLDAMPCT = one - CLOLDPCT
1337 : else
1338 199530390 : CLOLDPCT = zero
1339 199530390 : CLOLDAMPCT = zero
1340 : endif
1341 : !-----------------------------------------------------------------------
1342 : ! Remix everything for the next level
1343 : !-----------------------------------------------------------------------
1344 1496826296 : FCA = min( CFXX(LM1),FCXA*CLOLDPCT + CLNEWPCT*FCXB + AMCLPCT*FAX )
1345 1496826296 : if( FCA > zero ) then
1346 : !-----------------------------------------------------------------------
1347 : ! Maintain cloud core by reducing NC and AM area going into cloud below
1348 : !-----------------------------------------------------------------------
1349 1466565782 : RCA = (RCXA*FCXA*CLOLDPCT + RCXB*FCXB*CLNEWPCT + RAX*FAX*AMCLPCT)/FCA
1350 1466565782 : 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 1466565782 : if (RCA > zero) then
1364 : DCA = (RCXA*FCXA*CLOLDPCT)/(RCA*FCA)*DCXA + &
1365 : (RCXB*FCXB*CLNEWPCT)/(RCA*FCA)*DCXB + &
1366 1466565782 : (RAX*FAX*AMCLPCT)/(RCA*FCA)*DAX
1367 : else
1368 0 : DCA = zero
1369 0 : FCA = zero
1370 : endif
1371 :
1372 : else
1373 30260514 : FCA = zero
1374 30260514 : DCA = zero
1375 30260514 : RCA = zero
1376 : endif
1377 :
1378 1496826296 : FAMA = FCXA + FCXB + FAX - CFXX(LM1)
1379 1496826296 : if( FAMA > zero ) then
1380 767200904 : RAMA = (RCXA*FCXA*CLOLDAMPCT + RCXB*FCXB*CLNEWAMPCT + RAX*FAX*AMPCT)/FAMA
1381 767200904 : if( RAMA > zero ) then
1382 : DAMA = (RCXA*FCXA*CLOLDAMPCT)/(RAMA*FAMA)*DCXA + &
1383 : (RCXB*FCXB*CLNEWAMPCT)/(RAMA*FAMA)*DCXB + &
1384 762918582 : (RAX*FAX*AMPCT)/(RAMA*FAMA)*DAX
1385 : else
1386 4282322 : FAMA = zero
1387 4282322 : DAMA = zero
1388 : endif
1389 : else
1390 729625392 : FAMA = zero
1391 729625392 : DAMA = zero
1392 729625392 : RAMA = zero
1393 : endif
1394 : else upper_level
1395 106556386 : AMPCT = zero
1396 106556386 : AMCLPCT = zero
1397 106556386 : CLNEWPCT = zero
1398 106556386 : CLNEWAMPCT = zero
1399 106556386 : CLOLDPCT = zero
1400 106556386 : CLOLDAMPCT = zero
1401 : endif upper_level
1402 : else has_rls
1403 2370077286 : RNEW = zero
1404 2370077286 : QTEVAPCXA = QTTOPCA
1405 2370077286 : QTEVAPAX = QTTOPAA
1406 2370077286 : if( L > 1 ) then
1407 2348457544 : if( RLS(LM1) > zero ) then
1408 169117010 : 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 2179340534 : CFXX(LM1) = CFR(LM1)
1416 : endif
1417 : endif
1418 2370077286 : AMPCT = zero
1419 2370077286 : AMCLPCT = zero
1420 2370077286 : CLNEWPCT = zero
1421 2370077286 : CLNEWAMPCT = zero
1422 2370077286 : CLOLDPCT = zero
1423 2370077286 : CLOLDAMPCT = zero
1424 2370077286 : RCA = zero
1425 2370077286 : RAMA = zero
1426 2370077286 : FCA = zero
1427 2370077286 : FAMA = zero
1428 2370077286 : DCA = zero
1429 2370077286 : DAMA = zero
1430 : endif has_rls
1431 :
1432 3973459968 : 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 3973459968 : QTNETLCXA = QTRAINCXA + QTRIMECXA + QTWASHCXA - QTEVAPCXA
1440 3973459968 : QTNETLCXA = min( QTT(L)*FCXA,QTNETLCXA )
1441 :
1442 3973459968 : QTNETLCXB =QTRAINCXB
1443 3973459968 : QTNETLCXB = min( QTT(L)*FCXB,QTNETLCXB )
1444 :
1445 3973459968 : QTNETLAX = QTWASHAX - QTEVAPAX
1446 3973459968 : QTNETLAX = min( QTT(L)*FAX,QTNETLAX )
1447 :
1448 3973459968 : QTTNEW(L) = QTT(L) - (QTNETLCXA + QTNETLCXB + QTNETLAX)
1449 :
1450 3973459968 : 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 3973459968 : 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 3973459968 : 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 3973459968 : QTTOPCAX = (QTTOPCA + QTNETLCXA)*CLOLDPCT + QTNETLCXB*CLNEWPCT + (QTTOPAA + QTNETLAX)*AMCLPCT
1475 3973459968 : QTTOPAAX = (QTTOPCA + QTNETLCXA)*CLOLDAMPCT + QTNETLCXB*CLNEWAMPCT + (QTTOPAA + QTNETLAX)*AMPCT
1476 3973459968 : QTTOPCA = QTTOPCAX
1477 4101636096 : QTTOPAA = QTTOPAAX
1478 : end do level_loop
1479 :
1480 128176128 : 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 4102686720 : QTTJFL(:le,N) = QTTNEW(:le)
1544 :
1545 : end do species_loop
1546 : !
1547 1050624 : return
1548 72960 : end subroutine washo
1549 : !---------------------------------------------------------------------
1550 937295370 : 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 937295370 : if(TM .ge. TICE) then
1574 926438962 : QTDIS=(HSTAR*(QT/(QM*CFX))*0.029_r8*(PR/1.0e3_r8))*(CLWX*QM)
1575 10856408 : elseif (TM .le. TMIX) then
1576 10007666 : MUEMP=exp(-14.2252_r8+(1.55704e-1_r8*TM)-(7.1929e-4_r8*(TM**2.0_r8)))
1577 10007666 : QTDIS=MUEMP*(MOLMASS/18._r8)*(CLWX*QM)
1578 : else
1579 848742 : QTDIS=RETEFF*((HSTAR*(QT/(QM*CFX))*0.029_r8*(PR/1.0e3_r8))*(CLWX*QM))
1580 : endif
1581 :
1582 937295370 : return
1583 : end subroutine DISGAS
1584 :
1585 : !-----------------------------------------------------------------------
1586 390798644 : 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 390798644 : QTDISSTAR=(QTDIS*(QT*CFX))/(QTDIS+(QT*CFX))
1617 :
1618 : !---Tracer Loss frequency (1/s) within cloud fraction:
1619 390798644 : 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 390798644 : QTRAIN = QT*CFX*(1._r8 - exp(-DTSCAV*QTLF))
1624 :
1625 390798644 : return
1626 : end subroutine RAINGAS
1627 :
1628 :
1629 : !-----------------------------------------------------------------------
1630 503785301 : 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 503785301 : if ( BOXF == 0._r8 ) then
1670 233734 : QTWASH = 0._r8
1671 233734 : QTEVAP = 0._r8
1672 233734 : return
1673 : end if
1674 :
1675 : !---effective washout frequency (1/s):
1676 503551567 : 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 503551567 : QTMAX = QT*FWASH*DTSCAV
1679 503551567 : if (QTMAX .gt. QTRTOP) then
1680 : !---more of tracer T can go into rain
1681 384798700 : QTDIF = min (QT, QTMAX-QTRTOP)
1682 384798700 : QTWASH = QTDIF * (1._r8 - exp(-DTSCAV*FWASH))
1683 384798700 : QTEVAP=0._r8
1684 : else
1685 : !--too much of T in rain, must degas/evap T
1686 118752867 : QTWASH = 0._r8
1687 118752867 : QTEVAP = QTRTOP - QTMAX
1688 : endif
1689 :
1690 : return
1691 : end subroutine WASHGAS
1692 :
1693 : !-----------------------------------------------------------------------
1694 1218967148 : 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 1218967148 : RRATEX=RRATE*3600._r8 !mm/hr
1710 1218967148 : WX=CWATER*1.0e3_r8 !g/m3
1711 :
1712 1218967148 : if(RRATEX .gt. 0.04_r8) then
1713 103755022 : THETA=exp(-1.43_r8*dlog10(7._r8*RRATEX))+2.8_r8
1714 : else
1715 1115212126 : THETA=5._r8
1716 : endif
1717 1218967148 : PHI=RRATEX/(3600._r8*10._r8) !cgs units
1718 1218967148 : ETA=exp((3.01_r8*THETA)-10.5_r8)
1719 1218967148 : BETA=THETA/(1._r8+0.638_r8)
1720 1218967148 : ALPHA=exp(4._r8*(BETA-3.5_r8))
1721 1218967148 : BEE=(.638_r8*THETA/(1._r8+.638_r8))-1.0_r8
1722 1218967148 : GAMTHETA = shr_spfn_gamma(THETA)
1723 1218967148 : GAMBETA = shr_spfn_gamma(BETA+1._r8)
1724 : DEMPIRICAL=(((WX*ETA*GAMTHETA)/(1.0e6_r8*ALPHA*PHI*GAMBETA))** &
1725 1218967148 : (-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
|
