| blob | 2e1e91261ad142b7e3ea4e3293898027a44387dc |
1 MODULE module_cam_mam_addemiss
2 !WRF:MODEL_LAYER:CHEMICS
4 !----------------------------------------------------------------------
5 ! module_cam_mam_addemiss.F
6 !
7 ! adapted from module_mosaic_addemiss.F in nov, 2010 by r.c.easter
8 !
9 !----------------------------------------------------------------------
11 #include "MODAL_AERO_CPP_DEFINES.h"
13 ! rce 2005-feb-18 - one fix for indices of volumcen_sect, [now (isize,itype)]
14 ! rce 2005-jan-14 - added subr mosaic_seasalt_emiss (and a call to it)
15 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
16 ! variables in module_data_mosaic_asect
19 private
20 public :: cam_mam_addemiss
23 CONTAINS
27 !----------------------------------------------------------------------
28 subroutine cam_mam_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, &
29 config_flags, chem, slai, ust, smois, ivgtyp, isltyp, &
30 emis_ant,ebio_iso,ebio_olt,ebio_oli,rho_phy, &
31 dust_emiss_active, seasalt_emiss_active, &
32 ids,ide, jds,jde, kds,kde, &
33 ims,ime, jms,jme, kms,kme, &
34 its,ite, jts,jte, kts,kte )
35 !
36 ! adds emissions for cam-mam aerosol species
37 ! (i.e., emissions tendencies over time dtstep are applied
38 ! to the aerosol mass and number mixing ratios)
39 !
40 ! currently this routine
41 ! - applies emissions held in the emis_ant array
42 ! - does online sea-salt emissions using a modified version of the
43 ! mosaic sea-salt emissions routine
44 ! - does online dust emissions using a modified version of the
45 ! mosaic dust emissions routine
46 !
47 ! still to do
48 ! - for emis_ant emissions, implement emitted sizes similar to those used in cam
49 ! - implement the cam sea-salt and dust emissions routines
50 !
52 USE module_configure, only: grid_config_rec_type
53 USE module_state_description
54 ! USE module_state_description, only: num_chem, param_first_scalar, &
55 ! emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_cm,num_emis_ant
57 USE module_data_cam_mam_asect, only: ai_phase, &
58 dens_so4_aer, dens_nh4_aer, dens_pom_aer, &
59 dens_bc_aer, dens_dust_aer, dens_seas_aer, &
60 lptr_so4_aer, lptr_nh4_aer, lptr_pom_aer, &
61 lptr_bc_aer, lptr_dust_aer, lptr_seas_aer, &
62 maxd_asize, maxd_atype, nsize_aer, ntype_aer, numptr_aer
64 USE modal_aero_data, only: &
65 modeptr_accum, modeptr_aitken, modeptr_coarse, &
66 modeptr_coardust, modeptr_finedust, modeptr_pcarbon, &
67 modeptr_fineseas, modeptr_coarseas
69 USE module_cam_mam_init, only: &
70 pom_emit_1p4_factor, so4_emit_1p2_factor
72 USE module_data_sorgam, only: &
73 dgvem_i, dgvem_j, dgvem_c, sgem_i, sgem_j, sgem_c
75 IMPLICIT NONE
77 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
79 INTEGER, INTENT(IN ) :: id, &
80 ids,ide, jds,jde, kds,kde, &
81 ims,ime, jms,jme, kms,kme, &
82 its,ite, jts,jte, kts,kte
84 INTEGER, INTENT(IN) :: dust_emiss_active, seasalt_emiss_active
86 REAL, INTENT(IN ) :: dtstep
88 ! 10-m wind speed components (m/s)
89 REAL, DIMENSION( ims:ime , jms:jme ) , &
90 INTENT(IN ) :: u10, v10, xland, slai, ust
91 INTEGER, DIMENSION( ims:ime , jms:jme ) , &
92 INTENT(IN ) :: ivgtyp, isltyp
94 ! trace species mixing ratios (aerosol mass = ug/kg-air; number = #/kg-air)
95 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
96 INTENT(INOUT ) :: chem
97 !
98 ! aerosol emissions arrays ((ug/m3)*m/s)
99 !
100 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
101 REAL, DIMENSION( ims:ime, 1:config_flags%kemit, jms:jme,num_emis_ant ), &
102 INTENT(IN ) :: &
103 emis_ant
104 !jdf
105 REAL, DIMENSION( ims:ime , jms:jme ) , &
106 INTENT(IN ) :: ebio_iso, ebio_olt, ebio_oli
107 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
108 INTENT(IN ) :: rho_phy
109 !jdf
111 ! 1/(dry air density) and layer thickness (m)
112 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
113 INTENT(IN ) :: alt, dz8w
115 REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , &
116 INTENT(IN ) :: smois
118 ! local variables
120 ! anth_emiss_size_opt = 1 --> sizes follow those used in cam5
121 ! anth_emiss_size_opt = 2 --> sizes follow those used in sorgam
122 integer, parameter :: anth_emiss_size_opt = 1 ! user can change this value
123 ! dust_emiss_size_opt = 1 --> sizes follow those used in cam5
124 ! dust_emiss_size_opt = 2 --> sizes follow those used in sorgam
125 integer, parameter :: dust_emiss_size_opt = 1 ! user can change this value
126 ! seas_emiss_size_opt = 1 --> sizes follow those used in cam5
127 ! seas_emiss_size_opt = 2 --> sizes follow those used in sorgam
128 integer, parameter :: seas_emiss_size_opt = 1 ! user can change this value
130 integer :: i, iphase, isize, itype
131 integer :: itype_accum, itype_aitken, itype_pcarbon, &
132 itype_finedust, itype_coardust
133 integer :: j, k
134 integer :: l, lemit, lemit_type, lnumb
135 integer :: m, n
136 integer :: p1st
138 real, parameter :: pi = 3.1415926536
139 real :: dp_anth_emiss_aitken, dp_anth_emiss_accum, dp_anth_emiss_coarse
140 real :: dlo_seas_emiss( maxd_asize, maxd_atype ), &
141 dhi_seas_emiss( maxd_asize, maxd_atype )
142 real :: fact_mass, factbb_mass, factbb_numb
143 real :: tmp_dens, tmp_dp, tmp_vol1p
144 real :: total_soag, fact_soag
147 character(len=100) :: msg
150 iphase = ai_phase
153 !
154 ! code applies emissions in the emis_ant array,
155 ! then does "online" sea-salt and dust emissions
156 !
157 ! currently the code only works with
158 ! config_flags%emiss_inpt_opt = emiss_inpt_default OR emiss_inpt_pnnl_rs
159 !
161 emiss_inpt_select_1: SELECT CASE( config_flags%emiss_inpt_opt )
163 !jdf
164 CASE( emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_mam )
165 write(*,*) 'cam_mam_addemiss DOING emiss - emiss_inpt_opt =', &
166 config_flags%emiss_inpt_opt
167 !jdf
168 CASE DEFAULT
169 write(*,*) 'cam_mam_addemiss SKIPPING emiss - emiss_inpt_opt =', &
170 config_flags%emiss_inpt_opt
171 return
173 END SELECT emiss_inpt_select_1
176 p1st = param_first_scalar
178 itype_accum = modeptr_accum
179 itype_aitken = modeptr_aitken
180 #if ( defined MODAL_AERO_3MODE )
181 itype_pcarbon = modeptr_accum
182 itype_finedust = modeptr_accum
183 itype_coardust = modeptr_coarse
184 #elif ( defined MODAL_AERO_7MODE )
185 if ( (config_flags%chem_opt /= CBMZ_CAM_MAM7_NOAQ) .and. &
186 (config_flags%chem_opt /= CBMZ_CAM_MAM7_AQ ) ) then
187 call wrf_error_fatal( 'emission package for MAM7 is not implemented yet' )
188 end if
189 #else
190 itype_pcarbon = modeptr_pcarbon
191 itype_finedust = modeptr_finedust
192 itype_coardust = modeptr_coardust
193 #endif
195 ! set emitted sizes
196 if (anth_emiss_size_opt == 1) then
197 ! volume-mean diameter for emitted particles (cm)
198 dp_anth_emiss_aitken = 0.0504e-4
199 dp_anth_emiss_accum = 0.134e-4
200 dp_anth_emiss_coarse = 2.06e-4
201 else if (anth_emiss_size_opt == 2) then
202 ! volume-mean diameter for emitted particles (cm)
203 dp_anth_emiss_aitken = 1.0e2 * dgvem_i / exp( 1.5 * (log(sgem_i)**2) )
204 dp_anth_emiss_accum = 1.0e2 * dgvem_j / exp( 1.5 * (log(sgem_j)**2) )
205 dp_anth_emiss_coarse = 1.0e2 * dgvem_c / exp( 1.5 * (log(sgem_c)**2) )
206 else
207 write(msg,'(2a,i7)') 'subr cam_mam_addemiss', &
208 ' - illegal anth_emiss_size_opt = ', anth_emiss_size_opt
209 call wrf_error_fatal( msg )
210 end if
213 !jdf add case statements here to switch between 2 emission file options
214 emiss_inpt_select_2: SELECT CASE( config_flags%emiss_inpt_opt )
216 CASE( emiss_inpt_pnnl_rs )
217 !jdf
218 do 1900 lemit_type = 1, 11
220 iphase = 1
221 isize = 1
223 ! following if/then/else blocks determine
224 ! which chem species will receive a particular emis_ant(...,p_e_...)
225 ! additional mass factor for consistency with cam_mam assumptions
226 ! particle size and density for converting mass to number emissions
227 !
228 ! oc (org), ec, and pm25 - both the "i" and "j" go to the same
229 ! cam_mam mode
230 ! pm25 and pm_10 - go to the cam_mam dust species
231 ! no3 - ignore cam_mam does not treat no3
232 !
233 if (lemit_type == 1) then
234 lemit = p_e_so4i
235 itype = itype_aitken
236 l = lptr_so4_aer(isize,itype,iphase)
237 tmp_dens = dens_so4_aer
238 factbb_mass = so4_emit_1p2_factor
239 tmp_dp = dp_anth_emiss_aitken
240 else if (lemit_type == 2) then
241 lemit = p_e_so4j
242 itype = itype_accum
243 l = lptr_so4_aer(isize,itype,iphase)
244 tmp_dens = dens_so4_aer
245 ! so4_emit_1p2_factor is 115/96 if so4 (in 3-mode) is being treated
246 ! as ammonium bisulfate, or 1.0 if treated as pure sulfate
247 factbb_mass = so4_emit_1p2_factor
248 tmp_dp = dp_anth_emiss_accum
250 else if (lemit_type == 3) then
251 lemit = p_e_no3i
252 cycle ! cam_mam does not treat no3
253 else if (lemit_type == 4) then
254 lemit = p_e_no3j
255 cycle ! cam_mam does not treat no3
257 else if (lemit_type == 5) then
258 lemit = p_e_orgi
259 itype = itype_pcarbon
260 l = lptr_pom_aer(isize,itype,iphase)
261 tmp_dens = dens_pom_aer
262 ! pom_emit_1p4_factor is 1.4 if pom emissions are carbon-only emissions,
263 ! or 1.0 if they are organic matter emissions
264 factbb_mass = pom_emit_1p4_factor
265 tmp_dp = dp_anth_emiss_accum
266 else if (lemit_type == 6) then
267 lemit = p_e_orgj
268 itype = itype_pcarbon
269 l = lptr_pom_aer(isize,itype,iphase)
270 tmp_dens = dens_pom_aer
271 factbb_mass = pom_emit_1p4_factor
272 tmp_dp = dp_anth_emiss_accum
274 else if (lemit_type == 7) then
275 lemit = p_e_eci
276 itype = itype_pcarbon
277 l = lptr_bc_aer(isize,itype,iphase)
278 tmp_dens = dens_bc_aer
279 factbb_mass = 1.0
280 tmp_dp = dp_anth_emiss_accum
281 else if (lemit_type == 8) then
282 lemit = p_e_ecj
283 itype = itype_pcarbon
284 l = lptr_bc_aer(isize,itype,iphase)
285 tmp_dens = dens_bc_aer
286 factbb_mass = 1.0
287 tmp_dp = dp_anth_emiss_accum
289 else if (lemit_type == 9) then
290 lemit = p_e_pm25i
291 itype = itype_finedust
292 l = lptr_dust_aer(isize,itype,iphase)
293 tmp_dens = dens_dust_aer
294 factbb_mass = 1.0
295 tmp_dp = dp_anth_emiss_accum
296 else if (lemit_type == 10) then
297 lemit = p_e_pm25j
298 itype = itype_finedust
299 l = lptr_dust_aer(isize,itype,iphase)
300 tmp_dens = dens_dust_aer
301 factbb_mass = 1.0
302 tmp_dp = dp_anth_emiss_accum
303 else if (lemit_type == 11) then
304 lemit = p_e_pm_10
305 itype = itype_coardust
306 l = lptr_dust_aer(isize,itype,iphase)
307 tmp_dens = dens_dust_aer
308 factbb_mass = 1.0
309 tmp_dp = dp_anth_emiss_coarse
310 else
311 cycle
312 end if
314 if ((l < p1st) .or. (l > num_chem)) cycle
316 lnumb = numptr_aer(isize,itype,iphase)
317 if ((lnumb < p1st) .or. (lnumb > num_chem)) lnumb = -999888777
319 tmp_vol1p = (pi/6.0)*(tmp_dp**3) ! mean volume of emitted particles (cm3)
320 ! factbb_numb convert mass emissions (ug/m2/s) to number emissions (#/m2/s)
321 factbb_numb = 1.0e-6/(tmp_dens*tmp_vol1p)
323 do j = jts, jte
324 do k = 1, min(config_flags%kemit,kte)
325 do i = its, ite
327 ! emissions fluxes are ug/m2/s
328 ! mult by tmp_fact_mass to get mixing ratio change (ug/kg) over dtstep
329 fact_mass = (dtstep/dz8w(i,k,j))*alt(i,k,j)*factbb_mass
331 chem(i,k,j,l) = chem(i,k,j,l) + emis_ant(i,k,j,lemit)*fact_mass
333 if (lnumb > 0) then
334 chem(i,k,j,lnumb) = chem(i,k,j,lnumb) + &
335 emis_ant(i,k,j,lemit)*fact_mass*factbb_numb
336 end if
338 end do ! i
339 end do ! k
340 end do ! j
342 1900 continue
343 !
344 ! now do emissions for e_soag
345 ! add biogenic isoprene for soag_isoprene
346 ! no monoterpenes in CBMZ at present for soag_terpene
347 ! lump anthropogenic ald, hc3, hc5, hc8, ket, oli, olt, and ora2 which is par in CBMZ for soag_bigalk
348 ! lump biogenic olet and olei for soag_bigene
349 ! lump anthropogenic tol and xyl for soag_toluene
350 !
351 do j = jts, jte
352 do k = 1, 1
353 do i = its, ite
354 fact_soag = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
355 total_soag = (emis_ant(i,k,j,p_e_tol) * (92.14*0.15/12.0)) + &
356 (emis_ant(i,k,j,p_e_xyl) * (106.16*0.15/12.0)) + &
357 (emis_ant(i,k,j,p_e_hc5) * (77.6*0.05/12.0)) + &
358 (emis_ant(i,k,j,p_e_olt) * (72.34*0.05/12.0)) + &
359 (emis_ant(i,k,j,p_e_oli) * (75.78*0.05/12.0)) + &
360 (ebio_iso(i,j) * (68.11*0.04/12.0))
361 ! 0.4*emis_ant(i,k,j,p_e_ald) + 2.9*emis_ant(i,k,j,p_e_hc3) + &
362 ! 4.8*emis_ant(i,k,j,p_e_hc5) + 7.9*emis_ant(i,k,j,p_e_hc8) + &
363 ! 0.9*emis_ant(i,k,j,p_e_ket) + 2.8*emis_ant(i,k,j,p_e_oli) + &
364 ! 1.8*emis_ant(i,k,j,p_e_olt) + 1.0*emis_ant(i,k,j,p_e_ora2) + &
365 ! (ebio_iso(i,j) * (68.11*0.04/12.0)) + &
366 ! ebio_olt(i,j) + &
367 ! ebio_oli(i,j)
368 chem(i,k,j,p_soag) = chem(i,k,j,p_soag) + total_soag*fact_soag
369 end do ! i
370 end do ! k
371 end do ! j
372 do j = jts, jte
373 do k = 2, min(config_flags%kemit,kte)
374 do i = its, ite
375 fact_soag = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
376 total_soag = (emis_ant(i,k,j,p_e_tol) * (92.14*0.15/12.0)) + &
377 (emis_ant(i,k,j,p_e_xyl) * (106.16*0.15/12.0)) + &
378 (emis_ant(i,k,j,p_e_hc5) * (77.6*0.05/12.0)) + &
379 (emis_ant(i,k,j,p_e_olt) * (72.34*0.05/12.0)) + &
380 (emis_ant(i,k,j,p_e_oli) * (75.78*0.05/12.0))
381 ! emis_ant(i,k,j,p_e_xyl) + &
382 ! 0.4*emis_ant(i,k,j,p_e_ald) + 2.9*emis_ant(i,k,j,p_e_hc3) + &
383 ! 4.8*emis_ant(i,k,j,p_e_hc5) + 7.9*emis_ant(i,k,j,p_e_hc8) + &
384 ! 0.9*emis_ant(i,k,j,p_e_ket) + 2.8*emis_ant(i,k,j,p_e_oli) + &
385 ! 1.8*emis_ant(i,k,j,p_e_olt) + 1.0*emis_ant(i,k,j,p_e_ora2)
386 chem(i,k,j,p_soag) = chem(i,k,j,p_soag) + total_soag*fact_soag
387 end do ! i
388 end do ! k
389 end do ! j
391 CASE( emiss_inpt_pnnl_mam )
393 do 1910 lemit_type = 1, 14
395 iphase = 1
396 isize = 1
397 if (lemit_type == 1) then
398 lemit = p_e_so4i
399 itype = itype_aitken
400 l = lptr_so4_aer(isize,itype,iphase)
401 factbb_mass = so4_emit_1p2_factor
402 else if (lemit_type == 2) then
403 lemit = p_e_so4j
404 itype = itype_accum
405 l = lptr_so4_aer(isize,itype,iphase)
406 ! so4_emit_1p2_factor is 115/96 if so4 (in 3-mode) is being treated
407 ! as ammonium bisulfate, or 1.0 if treated as pure sulfate
408 factbb_mass = so4_emit_1p2_factor
409 else if (lemit_type == 3) then
410 lemit = p_e_orgj
411 itype = itype_pcarbon
412 l = lptr_pom_aer(isize,itype,iphase)
413 factbb_mass = pom_emit_1p4_factor
414 else if (lemit_type == 4) then
415 lemit = p_e_ecj
416 itype = itype_pcarbon
417 l = lptr_bc_aer(isize,itype,iphase)
418 factbb_mass = 1.0
419 else if (lemit_type == 5) then
420 lemit = p_e_dust_a1
421 itype = itype_finedust
422 l = lptr_dust_aer(isize,itype,iphase)
423 factbb_mass = 1.0
424 else if (lemit_type == 6) then
425 lemit = p_e_dust_a3
426 itype = itype_coardust
427 l = lptr_dust_aer(isize,itype,iphase)
428 factbb_mass = 1.0
429 else if (lemit_type == 7) then
430 lemit = p_e_ncl_a1
431 itype = itype_accum
432 l = lptr_seas_aer(isize,itype,iphase)
433 factbb_mass = 1.0
434 else if (lemit_type == 8) then
435 lemit = p_e_ncl_a2
436 itype = itype_aitken
437 l = lptr_seas_aer(isize,itype,iphase)
438 factbb_mass = 1.0
439 else if (lemit_type == 9) then
440 lemit = p_e_ncl_a3
441 itype = itype_coardust
442 l = lptr_seas_aer(isize,itype,iphase)
443 factbb_mass = 1.0
444 !PMA
445 !Currently my pre-processing script put num_a1 and num_a2 emissions from
446 !sea salt and dust into so4j_num and so4i_num
447 !num_a3 emissions are solely from dust and sea salt
448 !
450 else if (lemit_type == 10) then
451 lemit = p_e_so4i_num
452 itype = itype_aitken
453 l = numptr_aer(isize,itype,iphase)
454 else if (lemit_type == 11) then
455 lemit = p_e_so4j_num
456 itype = itype_accum
457 l = numptr_aer(isize,itype,iphase)
458 else if (lemit_type == 12) then
459 lemit = p_e_orgj_num
460 itype = itype_accum
461 l = numptr_aer(isize,itype,iphase)
462 else if (lemit_type == 13) then
463 lemit = p_e_ecj_num
464 itype = itype_accum
465 l = numptr_aer(isize,itype,iphase)
466 else if (lemit_type == 14) then
467 lemit = p_e_num_a3
468 itype = itype_coardust
469 l = numptr_aer(isize,itype,iphase)
470 else
471 cycle
472 end if
474 if ((l < p1st) .or. (l > num_chem)) cycle
476 do j = jts, jte
477 do k = 1, min(config_flags%kemit,kte)
478 do i = its, ite
480 ! emissions fluxes are ug/m2/s
481 ! mult by tmp_fact_mass to get mixing ratio change (ug/kg) over dtstep
482 fact_mass = (dtstep/dz8w(i,k,j))*alt(i,k,j)*factbb_mass
483 !PMA
484 factbb_numb = (dtstep/dz8w(i,k,j))*alt(i,k,j)
485 if (lemit_type < 10) then
486 chem(i,k,j,l) = chem(i,k,j,l) + emis_ant(i,k,j,lemit)*fact_mass
487 else
488 chem(i,k,j,l) = chem(i,k,j,l) + emis_ant(i,k,j,lemit)*factbb_numb
489 end if
491 end do ! i
492 end do ! k
493 end do ! j
495 1910 continue
496 !
497 ! units for e_soag* assumed to be ug/m2/s for now
498 !
499 do j = jts, jte
500 do k = 1, min(config_flags%kemit,kte)
501 do i = its, ite
502 fact_soag = (dtstep/dz8w(i,k,j))*alt(i,k,j)/1000.0
503 total_soag = emis_ant(i,k,j,p_e_soag_bigalk) + &
504 emis_ant(i,k,j,p_e_soag_bigene) + &
505 emis_ant(i,k,j,p_e_soag_isoprene) + &
506 emis_ant(i,k,j,p_e_soag_terpene) + &
507 emis_ant(i,k,j,p_e_soag_toluene)
508 chem(i,k,j,p_soag) = chem(i,k,j,p_soag) + total_soag*fact_soag
509 end do ! i
510 end do ! k
511 end do ! j
514 END SELECT emiss_inpt_select_2
515 !jdf end of add case statements here to switch between 2 emission file options
517 !
518 ! seasalt emissions
519 !
521 ! set emitted sizes
522 dlo_seas_emiss(:,:) = 0.0
523 dhi_seas_emiss(:,:) = 0.0
525 if (seas_emiss_size_opt == 1) then
526 #if ( defined MODAL_AERO_3MODE )
527 ! cut sizes are 0.02, 0.08, 1.0, 10.0 um
528 dlo_seas_emiss(1,modeptr_aitken ) = 0.02e-4
529 dhi_seas_emiss(1,modeptr_aitken ) = 0.08e-4
530 dlo_seas_emiss(1,modeptr_accum ) = 0.08e-4
531 dhi_seas_emiss(1,modeptr_accum ) = 1.00e-4
532 dlo_seas_emiss(1,modeptr_coarse ) = 1.00e-4
533 dhi_seas_emiss(1,modeptr_coarse ) = 10.0e-4
534 #else
535 ! cut sizes are 0.02, 0.08, 0.3, 1.0, 10.0 um
536 dlo_seas_emiss(1,modeptr_aitken ) = 0.02e-4
537 dhi_seas_emiss(1,modeptr_aitken ) = 0.08e-4
538 dlo_seas_emiss(1,modeptr_accum ) = 0.08e-4
539 dhi_seas_emiss(1,modeptr_accum ) = 0.30e-4
540 dlo_seas_emiss(1,modeptr_fineseas) = 0.30e-4
541 dhi_seas_emiss(1,modeptr_fineseas) = 1.00e-4
542 dlo_seas_emiss(1,modeptr_coarseas) = 1.00e-4
543 dhi_seas_emiss(1,modeptr_coarseas) = 10.0e-4
544 #endif
546 else if (seas_emiss_size_opt == 2) then
547 #if ( defined MODAL_AERO_3MODE )
548 ! cut sizes are 0.1, 1.25, 10.0 um and no aitken
549 dlo_seas_emiss(1,modeptr_accum ) = 0.10e-4
550 dhi_seas_emiss(1,modeptr_accum ) = 1.25e-4
551 dlo_seas_emiss(1,modeptr_coarse ) = 1.25e-4
552 dhi_seas_emiss(1,modeptr_coarse ) = 10.0e-4
553 #else
554 ! cut sizes are 0.1, 0.3, 1.25, 10.0 um and no aitken
555 dlo_seas_emiss(1,modeptr_accum ) = 0.10e-4
556 dhi_seas_emiss(1,modeptr_accum ) = 0.30e-4
557 dlo_seas_emiss(1,modeptr_fineseas) = 0.30e-4
558 dhi_seas_emiss(1,modeptr_fineseas) = 1.25e-4
559 dlo_seas_emiss(1,modeptr_coarseas) = 1.25e-4
560 dhi_seas_emiss(1,modeptr_coarseas) = 10.0e-4
561 #endif
563 else
564 write(msg,'(2a,i7)') 'subr cam_mam_addemiss', &
565 ' - illegal seas_emiss_size_opt = ', seas_emiss_size_opt
566 call wrf_error_fatal( msg )
567 end if
569 ! now do the sea-salt emissions
570 if (config_flags%seas_opt == 2) then
571 call cam_mam_mosaic_seasalt_emiss( &
572 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
573 dlo_seas_emiss, dhi_seas_emiss, &
574 ids,ide, jds,jde, kds,kde, &
575 ims,ime, jms,jme, kms,kme, &
576 its,ite, jts,jte, kts,kte )
577 end if
580 !
581 ! dust emissions
582 !
584 ! jdf: preliminary version that has not been made generic for situation
585 if (config_flags%dust_opt == 2) then
586 call cam_mam_mosaic_dust_emiss( &
587 slai, ust, smois, ivgtyp, isltyp, &
588 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
589 dust_emiss_size_opt, &
590 ids,ide, jds,jde, kds,kde, &
591 ims,ime, jms,jme, kms,kme, &
592 its,ite, jts,jte, kts,kte )
593 end if
596 return
599 END subroutine cam_mam_addemiss
603 !----------------------------------------------------------------------
604 subroutine cam_mam_mosaic_seasalt_emiss( &
605 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
606 dlo_seas_emiss, dhi_seas_emiss, &
607 ids,ide, jds,jde, kds,kde, &
608 ims,ime, jms,jme, kms,kme, &
609 its,ite, jts,jte, kts,kte )
610 !
611 ! adds seasalt emissions for mosaic aerosol species
612 ! (i.e., seasalt emissions tendencies over time dtstep are applied
613 ! to the aerosol mixing ratios)
614 !
616 USE module_configure, only: grid_config_rec_type
617 USE module_state_description, only: num_chem, param_first_scalar
618 USE module_data_cam_mam_asect, only: ai_phase, &
619 lptr_seas_aer, &
620 maxd_asize, maxd_atype, nsize_aer, ntype_aer, numptr_aer
621 ! USE module_mosaic_addemiss, only: seasalt_emitfactors_1bin
623 IMPLICIT NONE
625 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
627 INTEGER, INTENT(IN ) :: id, &
628 ids,ide, jds,jde, kds,kde, &
629 ims,ime, jms,jme, kms,kme, &
630 its,ite, jts,jte, kts,kte
632 REAL, INTENT(IN ) :: dtstep
634 ! 10-m wind speed components (m/s)
635 REAL, DIMENSION( ims:ime , jms:jme ), &
636 INTENT(IN ) :: u10, v10, xland
638 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
639 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
640 INTENT(INOUT ) :: chem
642 ! alt = 1.0/(dry air density) in (m3/kg)
643 ! dz8w = layer thickness in (m)
644 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
645 INTENT(IN ) :: alt, dz8w
647 ! emissions cut sizes (cm)
648 REAL, DIMENSION( maxd_asize, maxd_atype ), &
649 INTENT(IN ) :: dlo_seas_emiss, dhi_seas_emiss
651 ! local variables
652 integer i, j, k, l, l_seas, n
653 integer iphase, itype
654 integer p1st
656 real dum, dumdlo, dumdhi, dumoceanfrac, dumspd10
657 real factaa, factbb
659 real :: ssemfact_numb( maxd_asize, maxd_atype )
660 real :: ssemfact_mass( maxd_asize, maxd_atype )
662 write(*,*) 'in subr cam_mam_mosaic_seasalt_emiss'
663 p1st = PARAM_FIRST_SCALAR
665 ! for now just do itype=1
666 iphase = ai_phase
668 ! compute emissions factors for each size bin
669 ! (limit emissions to dp > 0.1 micrometer)
670 ssemfact_mass(:,:) = 0.0
671 ssemfact_numb(:,:) = 0.0
672 do itype = 1, ntype_aer
673 do n = 1, nsize_aer(itype)
674 dumdlo = max( dlo_seas_emiss(n,itype), 0.1e-4 )
675 dumdhi = max( dhi_seas_emiss(n,itype), 0.1e-4 )
676 if (dumdlo >= dumdhi) cycle
677 call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, &
678 ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) )
680 ! convert mass emissions factor from (g/m2/s) to (ug/m2/s)
681 ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6
682 end do
683 end do
686 ! loop over i,j and apply seasalt emissions
687 k = kts
688 do 1830 j = jts, jte
689 do 1820 i = its, ite
691 !Skip this point if over land. xland=1 for land and 2 for water.
692 !Also, there is no way to differentiate fresh from salt water.
693 !Currently, this assumes all water is salty.
694 if( xland(i,j) < 1.5 ) cycle
696 !wig: As far as I can tell, only real.exe knows the fractional breakdown
697 ! of land use. So, in wrf.exe, dumoceanfrac will always be 1.
698 dumoceanfrac = 1. !fraction of grid i,j that is salt water
699 dumspd10 = dumoceanfrac* &
700 ( (u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5*3.41) )
702 ! factaa is (s*m2/kg-air)
703 ! factaa*ssemfact_mass(n) is (s*m2/kg-air)*(ug/m2/s) = ug/kg-air
704 ! factaa*ssemfact_numb(n) is (s*m2/kg-air)*( #/m2/s) = #/kg-air
705 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
707 factbb = factaa * dumspd10
709 do 1815 itype = 1, ntype_aer
710 do 1810 n = 1, nsize_aer(itype)
711 if (ssemfact_mass(n,itype) <= 0.0) cycle
713 ! only apply emissions if bin has both na and cl species
714 l_seas = lptr_seas_aer(n,itype,iphase)
715 if (l_seas < p1st) cycle
717 chem(i,k,j,l_seas) = chem(i,k,j,l_seas) + &
718 factbb * ssemfact_mass(n,itype)
720 l = numptr_aer(n,itype,iphase)
721 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
722 factbb * ssemfact_numb(n,itype)
724 1810 continue
725 1815 continue
727 1820 continue
728 1830 continue
730 return
732 END subroutine cam_mam_mosaic_seasalt_emiss
736 !c----------------------------------------------------------------------
737 !c following is from gong06b.f in
738 !c /net/cirrus/files1/home/rce/oldfiles1/box/seasaltg
739 !c----------------------------------------------------------------------
740 subroutine seasalt_emitfactors_1bin( ireduce_smallr_emit, &
741 dpdrylo_cm, dpdryhi_cm, &
742 emitfact_numb, emitfact_surf, emitfact_mass )
743 !c
744 !c computes seasalt emissions factors for a specifed
745 !c dry particle size range
746 !c dpdrylo_cm = lower dry diameter (cm)
747 !c dpdryhi_cm = upper dry diameter (cm)
748 !c
749 !c number and mass emissions are then computed as
750 !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41)
751 !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41)
752 !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41)
753 !c
754 !c where spd10 = 10 m windspeed in m/s
755 !c
756 !c uses bubble emissions formula (eqn 5a) from
757 !c Gong et al. [JGR, 1997, p 3805-3818]
758 !c
759 !c *** for rdry < rdry_star, this formula overpredicts emissions.
760 !c A strictly ad hoc correction is applied to the formula,
761 !c based on sea-salt size measurements of
762 !c O'Dowd et al. [Atmos Environ, 1997, p 73-80]
763 !c
764 !c *** the correction is only applied when ireduce_smallr_emit > 0
765 !c
766 implicit none
768 !c subr arguments
769 integer ireduce_smallr_emit
770 real dpdrylo_cm, dpdryhi_cm, &
771 emitfact_numb, emitfact_surf, emitfact_mass
773 !c local variables
774 integer isub_bin, nsub_bin
776 real alnrdrylo
777 real drydens, drydens_43pi_em12, x_4pi_em8
778 real dum, dumadjust, dumb, dumexpb
779 real dumsum_na, dumsum_ma, dumsum_sa
780 real drwet, dlnrdry
781 real df0drwet, df0dlnrdry, df0dlnrdry_star
782 real relhum
783 real rdry, rdrylo, rdryhi, rdryaa, rdrybb
784 real rdrylowermost, rdryuppermost, rdry_star
785 real rwet, rwetaa, rwetbb
786 real rdry_cm, rwet_cm
787 real sigmag_star
788 real xmdry, xsdry
790 real pi
791 parameter (pi = 3.1415936536)
793 !c c1-c4 are constants for seasalt hygroscopic growth parameterization
794 !c in Eqn 3 and Table 2 of Gong et al. [1997]
795 real c1, c2, c3, c4, onethird
796 parameter (c1 = 0.7674)
797 parameter (c2 = 3.079)
798 parameter (c3 = 2.573e-11)
799 parameter (c4 = -1.424)
800 parameter (onethird = 1.0/3.0)
803 !c dry particle density (g/cm3)
804 drydens = 2.165
805 !c factor for radius (micrometers) to mass (g)
806 drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12
807 !c factor for radius (micrometers) to surface (cm2)
808 x_4pi_em8 = 4.0*pi*1.0e-8
809 !c bubble emissions formula assume 80% RH
810 relhum = 0.80
812 !c rdry_star = dry radius (micrometers) below which the
813 !c dF0/dr emission formula is adjusted downwards
814 rdry_star = 0.1
815 if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20
816 !c sigmag_star = geometric standard deviation used for
817 !c rdry < rdry_star
818 sigmag_star = 1.9
820 !c initialize sums
821 dumsum_na = 0.0
822 dumsum_sa = 0.0
823 dumsum_ma = 0.0
825 !c rdrylowermost, rdryuppermost = lower and upper
826 !c dry radii (micrometers) for overall integration
827 rdrylowermost = dpdrylo_cm*0.5e4
828 rdryuppermost = dpdryhi_cm*0.5e4
830 !c
831 !c "section 1"
832 !c integrate over rdry > rdry_star, where the dF0/dr emissions
833 !c formula is applicable
834 !c (when ireduce_smallr_emit <= 0, rdry_star = -1.0e20,
835 !c and the entire integration is done here)
836 !c
837 if (rdryuppermost .le. rdry_star) goto 2000
839 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
840 !c for this part of the integration
841 rdrylo = max( rdrylowermost, rdry_star )
842 rdryhi = rdryuppermost
844 nsub_bin = 1000
846 alnrdrylo = log( rdrylo )
847 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
849 !c compute rdry, rwet (micrometers) at lowest size
850 rdrybb = exp( alnrdrylo )
851 rdry_cm = rdrybb*1.0e-4
852 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
853 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
854 rwetbb = rwet_cm*1.0e4
856 do 1900 isub_bin = 1, nsub_bin
858 !c rdry, rwet at sub_bin lower boundary are those
859 !c at upper boundary of previous sub_bin
860 rdryaa = rdrybb
861 rwetaa = rwetbb
863 !c compute rdry, rwet (micrometers) at sub_bin upper boundary
864 dum = alnrdrylo + isub_bin*dlnrdry
865 rdrybb = exp( dum )
867 rdry_cm = rdrybb*1.0e-4
868 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
869 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
870 rwetbb = rwet_cm*1.0e4
872 !c geometric mean rdry, rwet (micrometers) for sub_bin
873 rdry = sqrt(rdryaa * rdrybb)
874 rwet = sqrt(rwetaa * rwetbb)
875 drwet = rwetbb - rwetaa
877 !c xmdry is dry mass in g
878 xmdry = drydens_43pi_em12 * (rdry**3.0)
880 !c xsdry is dry surface in cm2
881 xsdry = x_4pi_em8 * (rdry**2.0)
883 !c dumb is "B" in Gong's Eqn 5a
884 !c df0drwet is "dF0/dr" in Gong's Eqn 5a
885 dumb = ( 0.380 - log10(rwet) ) / 0.650
886 dumexpb = exp( -dumb*dumb)
887 df0drwet = 1.373 * (rwet**(-3.0)) * &
888 (1.0 + 0.057*(rwet**1.05)) * &
889 (10.0**(1.19*dumexpb))
891 dumsum_na = dumsum_na + drwet*df0drwet
892 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
893 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
895 1900 continue
898 !c
899 !c "section 2"
900 !c integrate over rdry < rdry_star, where the dF0/dr emissions
901 !c formula is just an extrapolation and predicts too many emissions
902 !c
903 !c 1. compute dF0/dln(rdry) = (dF0/drwet)*(drwet/dlnrdry)
904 !c at rdry_star
905 !c 2. for rdry < rdry_star, assume dF0/dln(rdry) is lognormal,
906 !c with the same lognormal parameters observed in
907 !c O'Dowd et al. [1997]
908 !c
909 2000 if (rdrylowermost .ge. rdry_star) goto 3000
911 !c compute dF0/dln(rdry) at rdry_star
912 rdryaa = 0.99*rdry_star
913 rdry_cm = rdryaa*1.0e-4
914 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
915 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
916 rwetaa = rwet_cm*1.0e4
918 rdrybb = 1.01*rdry_star
919 rdry_cm = rdrybb*1.0e-4
920 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
921 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
922 rwetbb = rwet_cm*1.0e4
924 rwet = 0.5*(rwetaa + rwetbb)
925 dumb = ( 0.380 - log10(rwet) ) / 0.650
926 dumexpb = exp( -dumb*dumb)
927 df0drwet = 1.373 * (rwet**(-3.0)) * &
928 (1.0 + 0.057*(rwet**1.05)) * &
929 (10.0**(1.19*dumexpb))
931 drwet = rwetbb - rwetaa
932 dlnrdry = log( rdrybb/rdryaa )
933 df0dlnrdry_star = df0drwet * (drwet/dlnrdry)
936 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
937 !c for this part of the integration
938 rdrylo = rdrylowermost
939 rdryhi = min( rdryuppermost, rdry_star )
941 nsub_bin = 1000
943 alnrdrylo = log( rdrylo )
944 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
946 do 2900 isub_bin = 1, nsub_bin
948 !c geometric mean rdry (micrometers) for sub_bin
949 dum = alnrdrylo + (isub_bin-0.5)*dlnrdry
950 rdry = exp( dum )
952 !c xmdry is dry mass in g
953 xmdry = drydens_43pi_em12 * (rdry**3.0)
955 !c xsdry is dry surface in cm2
956 xsdry = x_4pi_em8 * (rdry**2.0)
958 !c dumadjust is adjustment factor to reduce dF0/dr
959 dum = log( rdry/rdry_star ) / log( sigmag_star )
960 dumadjust = exp( -0.5*dum*dum )
962 df0dlnrdry = df0dlnrdry_star * dumadjust
964 dumsum_na = dumsum_na + dlnrdry*df0dlnrdry
965 dumsum_ma = dumsum_ma + dlnrdry*df0dlnrdry*xmdry
966 dumsum_sa = dumsum_sa + dlnrdry*df0dlnrdry*xsdry
968 2900 continue
971 !c
972 !c all done
973 !c
974 3000 emitfact_numb = dumsum_na
975 emitfact_mass = dumsum_ma
976 emitfact_surf = dumsum_sa
978 return
979 end subroutine seasalt_emitfactors_1bin
983 !----------------------------------------------------------------------
984 subroutine cam_mam_mosaic_dust_emiss( slai,ust, smois, ivgtyp, isltyp, &
985 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
986 dust_emiss_size_opt, &
987 ids,ide, jds,jde, kds,kde, &
988 ims,ime, jms,jme, kms,kme, &
989 its,ite, jts,jte, kts,kte )
990 !
991 ! adds dust emissions for mosaic aerosol species (i.e. emission tendencies
992 ! over time dtstep are applied to the aerosol mixing ratios)
993 !
994 ! This is a simple dust scheme based on Shaw et al. (2008) to appear in
995 ! Atmospheric Environment, recoded by Jerome Fast
996 !
997 ! NOTE:
998 ! 1) This version only works with the 8-bin version of MOSAIC.
999 ! 2) Dust added to MOSAIC's other inorganic specie, OIN. If Ca and CO3 are
1000 ! activated in the Registry, a small fraction also added to Ca and CO3.
1001 ! 3) The main departure from Shaw et al., is now alphamask is computed since
1002 ! the land-use categories in that paper and in WRF differ. WRF currently
1003 ! does not have that many land-use categories and adhoc assumptions had to
1004 ! be made. This version was tested for Mexico in the dry season. The main
1005 ! land-use categories in WRF that are likely dust sources are grass, shrub,
1006 ! and savannna (that WRF has in the desert regions of NW Mexico). Having
1007 ! dust emitted from these types for other locations and other times of the
1008 ! year is not likely to be valid.
1009 ! 4) An upper bound on ustar was placed because the surface parameterizations
1010 ! in WRF can produce unrealistically high values that lead to very high
1011 ! dust emission rates.
1012 ! 5) Other departures' from Shaw et al. noted below, but are probably not as
1013 ! important as 2) and 3).
1014 !
1016 USE module_configure, only: grid_config_rec_type
1017 USE module_state_description, only: num_chem, param_first_scalar
1019 USE module_data_cam_mam_asect, only: ai_phase, &
1020 lptr_dust_aer, ntype_aer, numptr_aer
1022 USE modal_aero_data, only: &
1023 modeptr_accum, modeptr_coarse, modeptr_coardust, modeptr_finedust
1025 IMPLICIT NONE
1027 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1029 INTEGER, INTENT(IN ) :: id, &
1030 dust_emiss_size_opt, &
1031 ids,ide, jds,jde, kds,kde, &
1032 ims,ime, jms,jme, kms,kme, &
1033 its,ite, jts,jte, kts,kte
1035 REAL, INTENT(IN ) :: dtstep
1037 ! 10-m wind speed components (m/s)
1038 REAL, DIMENSION( ims:ime , jms:jme ), &
1039 INTENT(IN ) :: u10, v10, xland, slai, ust
1040 INTEGER, DIMENSION( ims:ime , jms:jme ), &
1041 INTENT(IN ) :: ivgtyp, isltyp
1043 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
1044 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
1045 INTENT(INOUT ) :: chem
1047 ! alt = 1.0/(dry air density) in (m3/kg)
1048 ! dz8w = layer thickness in (m)
1049 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1050 INTENT(IN ) :: alt, dz8w
1052 REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , &
1053 INTENT(IN ) :: smois
1055 ! local variables
1056 integer, parameter :: max_dust_mode = 2
1057 integer :: i, ii, j, k, l, l_dust, l_numb, n
1058 integer :: imode, iphase, isize, itype, izob
1059 integer :: p1st
1060 integer :: isize_dust_mode(max_dust_mode)
1061 integer :: itype_dust_mode(max_dust_mode)
1062 integer :: n_dust_mode
1064 real :: dum, dumdlo, dumdhi, dumlandfrac, dumspd10
1065 real :: factaa, factbb, factcc, fracoin, fracca, fracco3, fractot
1066 real :: ustart, ustar1, ustart0
1067 real :: alphamask, f8, f50, f51, f52, wetfactor, sumdelta, ftot
1068 real :: smois_grav, wp, pclay
1069 real :: beta(4,7)
1070 real :: gamma(4), delta(4)
1071 real :: sz(8)
1072 real :: dustflux, densdust
1073 real :: mass1part_mos8bin(8)
1074 real :: sz_wght_dust_mode(8,max_dust_mode)
1075 real :: dcen_mos8bin(8)
1077 character(len=100) :: msg
1079 write(*,*) 'in subr cam_mam_mosaic_dust_emiss'
1080 p1st = param_first_scalar
1083 ! the mass emission fractions [sz(1:8)] were derived for the
1084 ! mosaic 8-bin size bins
1085 ! the emissions for each mosaic bin are assigned to the cam-mam fine or coarse mode
1086 ! (e.g, the fine mode may get bins 1-5 and the coarse mode bins 6-8)
1087 ! for calculating number emissions from mass emissions, use the mosaic bin
1088 ! sizes (and corresponding 1-particle masses) in order to get number emissions
1089 ! consistent with mosaic
1090 !
1091 ! note: the sz() values were recommend by Natalie Mahowold, so probably follow the
1092 ! assumptions used in the "dead" dust module.
1093 ! applying those assumptions along with the cut sizes might be a better approach
1094 ! than adapting the mosaic 8-bin sz() values
1096 #if (defined MODAL_AERO)
1097 ! 3 modes -- dust is in accum and coarse modes
1098 n_dust_mode = 2
1099 itype_dust_mode(1) = modeptr_accum
1100 itype_dust_mode(2) = modeptr_coarse
1101 isize_dust_mode(:) = 1
1102 if (dust_emiss_size_opt == 1) then
1103 ! cam5 cut sizes for dust emiss are 0.1, 1.0, 10.0 um
1104 ! fine mode gets bins 1-4 and 60% of bin 5
1105 sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 0.6, 0., 0., 0. /)
1106 sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0.4, 1., 1., 1. /)
1107 else if (dust_emiss_size_opt == 2) then
1108 ! sorgam cut sizes for dust emiss are 0.1, 2.5, 10.0 um
1109 ! fine mode gets bins 1-6
1110 sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 1., 1., 0., 0. /)
1111 sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0., 0., 1., 1. /)
1112 else
1113 write(msg,'(2a,i7)') 'subr cam_mam_mosaic_dust_emiss', &
1114 ' - illegal dust_emiss_size_opt = ', dust_emiss_size_opt
1115 call wrf_error_fatal( msg )
1116 end if
1117 #else
1118 ! 7 modes -- dust is in fine-dust and coarse-dust modes
1119 n_dust_mode = 2
1120 itype_dust_mode(1) = modeptr_finedust
1121 itype_dust_mode(2) = modeptr_coardust
1122 isize_dust_mode(:) = 1
1123 if (dust_emiss_size_opt == 1) then
1124 ! cam5 cut sizes for dust emiss are 0.1, 2.0, 10.0 um
1125 ! fine mode gets bins 1-5 and 60% of bin 6
1126 ! (for 7-mode cam, the fine-dust mode is larger than the accumulation mode,
1127 ! so the cut size is larger than with 3-mode cam)
1128 sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 1., 0.6, 0., 0. /)
1129 sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0., 0.4, 1., 1. /)
1130 else if (dust_emiss_size_opt == 2) then
1131 ! sorgam cut sizes for dust emiss are 0.1, 2.5, 10.0 um
1132 ! fine mode gets bins 1-6
1133 sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 1., 1., 0., 0. /)
1134 sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0., 0., 1., 1. /)
1135 else
1136 write(msg,'(2a,i7)') 'subr cam_mam_mosaic_dust_emiss', &
1137 'illegal dust_emiss_size_opt = ', dust_emiss_size_opt
1138 call wrf_error_fatal( msg )
1139 end if
1140 #endif
1142 ! bin center diameters for mosaic-8bin (cm)
1143 dcen_mos8bin(8) = 1.0e-4*(10.0/sqrt(2.0))
1144 do n = 7, 1, -1
1145 dcen_mos8bin(n) = dcen_mos8bin(n+1)*0.5
1146 end do
1147 ! mass1part_mos8bin is mass of a single particle in ug,
1148 ! assuming density of dust ~2.5 g cm-3
1149 densdust=2.5
1150 do n = 1, 8
1151 mass1part_mos8bin(n)=0.523598*(dcen_mos8bin(n)**3)*densdust*1.0e06
1152 end do
1155 !
1156 ! from Nickovic et al., JGR, 2001 and Shaw et al. 2007
1157 ! beta: fraction of clay, small silt, large silt, and sand correcsponding to Zobler class (7)
1158 ! beta (1,*) for 0.5-1 um
1159 ! beta (2,*) for 1-10 um
1160 ! beta (3,*) for 10-25 um
1161 ! beta (4,*) for 25-50 um
1162 !
1163 beta(1,1)=0.12
1164 beta(2,1)=0.04
1165 beta(3,1)=0.04
1166 beta(4,1)=0.80
1167 beta(1,2)=0.34
1168 beta(2,2)=0.28
1169 beta(3,2)=0.28
1170 beta(4,2)=0.10
1171 beta(1,3)=0.45
1172 beta(2,3)=0.15
1173 beta(3,3)=0.15
1174 beta(4,3)=0.25
1175 beta(1,4)=0.12
1176 beta(2,4)=0.09
1177 beta(3,4)=0.09
1178 beta(4,4)=0.70
1179 beta(1,5)=0.40
1180 beta(2,5)=0.05
1181 beta(3,5)=0.05
1182 beta(4,5)=0.50
1183 beta(1,6)=0.34
1184 beta(2,6)=0.18
1185 beta(3,6)=0.18
1186 beta(4,6)=0.30
1187 beta(1,7)=0.22
1188 beta(2,7)=0.09
1189 beta(3,7)=0.09
1190 beta(4,7)=0.60
1191 gamma(1)=0.08
1192 gamma(2)=1.00
1193 gamma(3)=1.00
1194 gamma(4)=0.12
1195 !
1196 ! * Mass fractions for each size bin. These values were recommended by
1197 ! Natalie Mahowold, with bins 7 and 8 the same as bins 3 and 4 from CAM.
1198 ! * Changed slightly since Natelie's estimates do not add up to 1.0
1199 ! * This would need to be made more generic for other bin sizes.
1200 ! sz(1)=0
1201 ! sz(2)=1.78751e-06
1202 ! sz(3)=0.000273786
1203 ! sz(4)=0.00847978
1204 ! sz(5)=0.056055
1205 ! sz(6)=0.0951896
1206 ! sz(7)=0.17
1207 ! sz(8)=0.67
1208 !jdf sz(1)=0.0
1209 !jdf sz(2)=0.0
1210 !jdf sz(3)=0.0005
1211 !jdf sz(4)=0.0095
1212 !jdf sz(5)=0.03
1213 !jdf sz(6)=0.10
1214 !jdf sz(7)=0.18
1215 !jdf sz(8)=0.68
1216 sz(1)=0.0
1217 sz(2)=0.0
1218 sz(3)=0.0005
1219 sz(4)=0.0095
1220 sz(5)=0.01
1221 sz(6)=0.06
1222 sz(7)=0.20
1223 sz(8)=0.72
1225 iphase = ai_phase
1227 ! loop over i,j and apply dust emissions
1228 k = kts
1229 do 1830 j = jts, jte
1230 do 1820 i = its, ite
1232 if( xland(i,j) > 1.5 ) cycle
1234 ! compute wind speed anyway, even though ustar is used below
1236 dumlandfrac = 1.
1237 dumspd10=(u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5)
1238 if(dumspd10 >= 5.0) then
1239 dumspd10 = dumlandfrac* &
1240 ( dumspd10*dumspd10*(dumspd10-5.0))
1241 else
1242 dumspd10=0.
1243 endif
1245 ! part1 - compute vegetation mask
1246 !
1247 ! * f8, f50, f51, f52 refer to vegetation classes from the Olsen categories
1248 ! for desert, sand desert, grass semi-desert, and shrub semi-desert
1249 ! * in WRF, vegetation types 7, 8 and 10 are grassland, shrubland, and savanna
1250 ! that are dominate types in Mexico and probably have some erodable surface
1251 ! during the dry season
1252 ! * currently modified these values so that only a small fraction of cell
1253 ! area is erodable
1254 ! * these values are highly tuneable!
1256 alphamask=0.001
1257 if (ivgtyp(i,j) .eq. 7) then
1258 f8=0.005
1259 f50=0.00
1260 f51=0.10
1261 f51=0.066
1262 f52=0.00
1263 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
1264 endif
1265 if (ivgtyp(i,j) .eq. 8) then
1266 f8=0.010
1267 f50=0.00
1268 f51=0.00
1269 f52=0.15
1270 f52=0.10
1271 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
1272 endif
1273 if (ivgtyp(i,j) .eq. 10) then
1274 f8=0.00
1275 f50=0.00
1276 f51=0.01
1277 f52=0.00
1278 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
1279 endif
1281 ! part2 - zobler
1282 !
1283 ! * in Shaw's paper, dust is computed for 4 size ranges:
1284 ! 0.5-1 um
1285 ! 1-10 um
1286 ! 10-25 um
1287 ! 25-50 um
1288 ! * Shaw's paper also accounts for sub-grid variability in soil
1289 ! texture, but here we just assume the same soil texture for each
1290 ! grid cell
1291 ! * since MOSAIC is currently has a maximum size range up to 10 um,
1292 ! neglect upper 2 size ranges and lowest size range (assume small)
1293 ! * map WRF soil classes arbitrarily to Zolber soil textural classes
1294 ! * skip dust computations for WRF soil classes greater than 13, i.e.
1295 ! do not compute dust over water, bedrock, and other surfaces
1296 ! * should be skipping for water surface at this point anyway
1297 !
1298 izob=0
1299 if(isltyp(i,j).eq.1) izob=1
1300 if(isltyp(i,j).eq.2) izob=1
1301 if(isltyp(i,j).eq.3) izob=4
1302 if(isltyp(i,j).eq.4) izob=2
1303 if(isltyp(i,j).eq.5) izob=2
1304 if(isltyp(i,j).eq.6) izob=2
1305 if(isltyp(i,j).eq.7) izob=7
1306 if(isltyp(i,j).eq.8) izob=2
1307 if(isltyp(i,j).eq.9) izob=6
1308 if(isltyp(i,j).eq.10) izob=5
1309 if(isltyp(i,j).eq.11) izob=2
1310 if(isltyp(i,j).eq.12) izob=3
1311 if(isltyp(i,j).ge.13) izob=0
1312 if(izob.eq.0) goto 1820
1313 !
1314 ! part3 - dustprod
1315 !
1316 do ii=1,4
1317 delta(ii)=0.0
1318 enddo
1319 sumdelta=0.0
1320 do ii=1,4
1321 delta(ii)=beta(ii,izob)*gamma(ii)
1322 if(ii.lt.4) then
1323 sumdelta=sumdelta+delta(ii)
1324 endif
1325 enddo
1326 do ii=1,4
1327 delta(ii)=delta(ii)/sumdelta
1328 enddo
1330 ! part4 - wetness
1331 !
1332 ! * assume dry for now, have passed in soil moisture to this routine
1333 ! but needs to be included here
1334 ! * wetfactor less than 1 would reduce dustflux
1335 ! * convert model soil moisture (m3/m3) to gravimetric soil moisture
1336 ! (mass of water / mass of soil in %) assuming a constant density
1337 ! for soil
1338 pclay=beta(1,izob)*100.
1339 wp=0.0014*pclay*pclay+0.17*pclay
1340 smois_grav=(smois(i,1,j)/2.6)*100.
1341 if(smois_grav.gt.wp) then
1342 wetfactor=sqrt(1.0+1.21*(smois_grav-wp)**0.68)
1343 else
1344 wetfactor=1.0
1345 endif
1346 ! wetfactor=1.0
1348 ! part5 - dustflux
1349 ! lower bound on ustar = 20 cm/s as in Shaw et al, but set upper
1350 ! bound to 100 cm/s
1352 ustar1=ust(i,j)*100.0
1353 if(ustar1.gt.100.0) ustar1=100.0
1354 ustart0=20.0
1355 ustart=ustart0*wetfactor
1356 if(ustar1.le.ustart) then
1357 dustflux=0.0
1358 else
1359 dustflux=1.0e-14*(ustar1**4)*(1.0-(ustart/ustar1))
1360 endif
1361 dustflux=dustflux*10.0
1362 ! units kg m-2 s-1
1363 ftot=0.0
1364 do ii=1,2
1365 ftot=ftot+dustflux*alphamask*delta(ii)
1366 enddo
1367 ! convert to ug m-2 s-1
1368 ftot=ftot*1.0e+09
1370 ! apportion other inorganics only
1371 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
1372 factbb = factaa * ftot
1373 fractot = 1.0
1375 do imode = 1, n_dust_mode
1376 ! loop over the cam-mam modes that have dust
1377 itype = itype_dust_mode(imode)
1378 isize = isize_dust_mode(imode)
1380 l_dust = lptr_dust_aer(isize,itype,iphase)
1381 if ((l_dust < p1st) .or. (l_dust > num_chem)) cycle
1382 l_numb = numptr_aer(isize,itype,iphase)
1383 if ((l_numb < p1st) .or. (l_numb > num_chem)) l_numb = -1
1385 ! if (ivgtyp(i,j) .eq. 8) print*,'jdf',i,j,ustar1,ustart0,factaa,ftot
1386 ! do 1810 n = 1, nsize_aer(itype)
1387 do n = 1, 8
1388 ! calculate contribution of each mosaic-8bin bin to the current cam-mam mode
1389 if (sz_wght_dust_mode(n,imode) <= 0.0) cycle
1390 factcc = factbb * sz(n) * fractot * sz_wght_dust_mode(n,imode)
1392 chem(i,k,j,l_dust) = chem(i,k,j,l_dust) + factcc
1394 if (l_numb >= p1st) &
1395 chem(i,k,j,l_numb) = chem(i,k,j,l_numb) + factcc/mass1part_mos8bin(n)
1397 end do ! n
1399 end do ! imode
1402 1820 continue
1403 1830 continue
1405 return
1407 END subroutine cam_mam_mosaic_dust_emiss
1410 !----------------------------------------------------------------------
1411 END MODULE module_cam_mam_addemiss