Line data Source code
1 : module mo_drydep
2 :
3 : !---------------------------------------------------------------------
4 : ! ... Dry deposition
5 : !---------------------------------------------------------------------
6 :
7 : use shr_kind_mod, only : r8 => shr_kind_r8, shr_kind_cl
8 : use chem_mods, only : gas_pcnst
9 : use pmgrid, only : plev
10 : use spmd_utils, only : masterproc
11 : use ppgrid, only : pcols, begchunk, endchunk
12 : use mo_tracname, only : solsym
13 : use cam_abortutils, only : endrun
14 : use ioFileMod, only : getfil
15 : use pio
16 : use cam_pio_utils, only : cam_pio_openfile, cam_pio_closefile
17 : use cam_logfile, only : iulog
18 : use dyn_grid, only : get_dyn_grid_parm, get_horiz_grid_d
19 : use scamMod, only : single_column
20 :
21 : use shr_drydep_mod, only : nddvels => n_drydep, drydep_list, mapping
22 : use physconst, only : karman
23 :
24 : use infnan, only : nan, assignment(=)
25 :
26 : implicit none
27 :
28 : save
29 :
30 : interface drydep_inti
31 : module procedure dvel_inti_xactive
32 : end interface
33 :
34 : interface drydep
35 : module procedure drydep_fromlnd
36 : end interface
37 :
38 : private
39 :
40 : public :: drydep_inti, drydep, has_drydep
41 : public :: drydep_update
42 : public :: n_land_type, fraction_landuse, drydep_srf_file
43 :
44 : integer :: pan_ndx, mpan_ndx, o3_ndx, ch4_ndx, co_ndx, h2_ndx, ch3cooh_ndx
45 : integer :: sogm_ndx, sogi_ndx, sogt_ndx, sogb_ndx, sogx_ndx
46 :
47 : integer :: so2_ndx, ch3cn_ndx, hcn_ndx, hcooh_ndx
48 :
49 : integer :: o3a_ndx,xpan_ndx,xmpan_ndx
50 :
51 : integer :: cohc_ndx=-1, come_ndx=-1
52 : integer, parameter :: NTAGS = 50
53 : integer :: cotag_ndx(NTAGS)
54 : integer :: tag_cnt
55 :
56 : real(r8), parameter :: small_value = 1.e-36_r8
57 : real(r8), parameter :: large_value = 1.e36_r8
58 : real(r8), parameter :: diffm = 1.789e-5_r8
59 : real(r8), parameter :: diffk = 1.461e-5_r8
60 : real(r8), parameter :: difft = 2.060e-5_r8
61 : real(r8), parameter :: vonkar = karman
62 : real(r8), parameter :: ric = 0.2_r8
63 : real(r8), parameter :: r = 287.04_r8
64 : real(r8), parameter :: cp = 1004._r8
65 : real(r8), parameter :: grav = 9.81_r8
66 : real(r8), parameter :: p00 = 100000._r8
67 : real(r8), parameter :: wh2o = 18.0153_r8
68 : real(r8), parameter :: ph = 1.e-5_r8
69 : real(r8), parameter :: ph_inv = 1._r8/ph
70 : real(r8), parameter :: rovcp = r/cp
71 :
72 : logical, public :: has_dvel(gas_pcnst) = .false.
73 : integer :: map_dvel(gas_pcnst) = 0
74 :
75 : real(r8), protected, allocatable :: fraction_landuse(:,:,:)
76 : real(r8), allocatable, dimension(:,:,:) :: dep_ra ! [s/m] aerodynamic resistance
77 : real(r8), allocatable, dimension(:,:,:) :: dep_rb ! [s/m] resistance across sublayer
78 : integer, parameter :: n_land_type = 11
79 :
80 : integer, allocatable :: spc_ndx(:) ! nddvels
81 : real(r8), public :: crb
82 :
83 : type lnd_dvel_type
84 : real(r8), pointer :: dvel(:,:) ! deposition velocity over land (cm/s)
85 : end type lnd_dvel_type
86 :
87 : type(lnd_dvel_type), allocatable :: lnd(:)
88 : character(len=SHR_KIND_CL) :: drydep_srf_file
89 :
90 : contains
91 :
92 : !---------------------------------------------------------------------------
93 : !---------------------------------------------------------------------------
94 0 : subroutine dvel_inti_fromlnd
95 : use mo_chem_utls, only : get_spc_ndx
96 : use cam_abortutils, only : endrun
97 :
98 : integer :: ispc
99 :
100 0 : allocate(spc_ndx(nddvels))
101 0 : allocate( lnd(begchunk:endchunk) )
102 :
103 0 : do ispc = 1,nddvels
104 :
105 0 : spc_ndx(ispc) = get_spc_ndx(drydep_list(ispc))
106 0 : if (spc_ndx(ispc) < 1) then
107 0 : write(*,*) 'drydep_inti: '//trim(drydep_list(ispc))//' is not included in species set'
108 0 : call endrun('drydep_init: invalid dry deposition species')
109 : endif
110 :
111 : enddo
112 :
113 0 : crb = (difft/diffm)**(2._r8/3._r8) !.666666_r8
114 :
115 0 : endsubroutine dvel_inti_fromlnd
116 :
117 : !-------------------------------------------------------------------------------------
118 : !-------------------------------------------------------------------------------------
119 0 : subroutine drydep_update( state, cam_in )
120 : use physics_types, only : physics_state
121 : use camsrfexch, only : cam_in_t
122 :
123 : type(physics_state), intent(in) :: state ! Physics state variables
124 : type(cam_in_t), intent(in) :: cam_in
125 :
126 0 : if (nddvels<1) return
127 :
128 0 : lnd(state%lchnk)%dvel => cam_in%depvel
129 :
130 0 : end subroutine drydep_update
131 :
132 : !-------------------------------------------------------------------------------------
133 : !-------------------------------------------------------------------------------------
134 0 : subroutine drydep_fromlnd( ocnfrac, icefrac, sfc_temp, pressure_sfc, &
135 : wind_speed, spec_hum, air_temp, pressure_10m, rain, &
136 0 : snow, solar_flux, dvelocity, dflx, mmr, &
137 : tv, ncol, lchnk )
138 :
139 : !-------------------------------------------------------------------------------------
140 : ! combines the deposition velocities provided by the land model with deposition
141 : ! velocities over ocean and sea ice
142 : !-------------------------------------------------------------------------------------
143 :
144 0 : use ppgrid, only : pcols
145 : use chem_mods, only : gas_pcnst
146 :
147 : #if (defined OFFLINE_DYN)
148 : use metdata, only: get_met_fields
149 : #endif
150 :
151 : !-------------------------------------------------------------------------------------
152 : ! ... dummy arguments
153 : !-------------------------------------------------------------------------------------
154 :
155 : real(r8), intent(in) :: icefrac(pcols)
156 : real(r8), intent(in) :: ocnfrac(pcols)
157 : integer, intent(in) :: ncol
158 : integer, intent(in) :: lchnk ! chunk number
159 : real(r8), intent(in) :: sfc_temp(pcols) ! surface temperature (K)
160 : real(r8), intent(in) :: pressure_sfc(pcols) ! surface pressure (Pa)
161 : real(r8), intent(in) :: wind_speed(pcols) ! 10 meter wind speed (m/s)
162 : real(r8), intent(in) :: spec_hum(pcols) ! specific humidity (kg/kg)
163 : real(r8), intent(in) :: air_temp(pcols) ! surface air temperature (K)
164 : real(r8), intent(in) :: pressure_10m(pcols) ! 10 meter pressure (Pa)
165 : real(r8), intent(in) :: rain(pcols)
166 : real(r8), intent(in) :: snow(pcols) ! snow height (m)
167 : real(r8), intent(in) :: solar_flux(pcols) ! direct shortwave radiation at surface (W/m^2)
168 : real(r8), intent(in) :: tv(pcols) ! potential temperature
169 : real(r8), intent(in) :: mmr(pcols,plev,gas_pcnst) ! constituent concentration (kg/kg)
170 : real(r8), intent(out) :: dvelocity(ncol,gas_pcnst) ! deposition velocity (cm/s)
171 : real(r8), intent(inout) :: dflx(pcols,gas_pcnst) ! deposition flux (/cm^2/s)
172 :
173 : !-------------------------------------------------------------------------------------
174 : ! ... local variables
175 : !-------------------------------------------------------------------------------------
176 0 : real(r8) :: ocnice_dvel(ncol,gas_pcnst)
177 : real(r8) :: ocnice_dflx(pcols,gas_pcnst)
178 :
179 0 : real(r8), dimension(ncol) :: term ! work array
180 : integer :: ispec
181 : real(r8) :: lndfrac(pcols)
182 : #if (defined OFFLINE_DYN)
183 : real(r8) :: met_ocnfrac(pcols)
184 : real(r8) :: met_icefrac(pcols)
185 : #endif
186 : integer :: i
187 :
188 0 : lndfrac(:ncol) = 1._r8 - ocnfrac(:ncol) - icefrac(:ncol)
189 :
190 0 : where( lndfrac(:ncol) < 0._r8 )
191 : lndfrac(:ncol) = 0._r8
192 : endwhere
193 :
194 : #if (defined OFFLINE_DYN)
195 : call get_met_fields(lndfrac, met_ocnfrac, met_icefrac, lchnk, ncol)
196 : #endif
197 :
198 : !-------------------------------------------------------------------------------------
199 : ! ... initialize
200 : !-------------------------------------------------------------------------------------
201 : dvelocity(:,:) = 0._r8
202 :
203 : !-------------------------------------------------------------------------------------
204 : ! ... compute the dep velocities over ocean and sea ice
205 : ! land type 7 is used for ocean
206 : ! land type 8 is used for sea ice
207 : !-------------------------------------------------------------------------------------
208 : call drydep_xactive( sfc_temp, pressure_sfc, &
209 : wind_speed, spec_hum, air_temp, pressure_10m, rain, &
210 : snow, solar_flux, ocnice_dvel, ocnice_dflx, mmr, &
211 : tv, ncol, lchnk, &
212 : #if (defined OFFLINE_DYN)
213 : ocnfrc=met_ocnfrac,icefrc=met_icefrac, beglandtype=7, endlandtype=8 )
214 : #else
215 0 : ocnfrc=ocnfrac,icefrc=icefrac, beglandtype=7, endlandtype=8 )
216 : #endif
217 0 : term(:ncol) = 1.e-2_r8 * pressure_10m(:ncol) / (r*tv(:ncol))
218 :
219 0 : do ispec = 1,nddvels
220 : !-------------------------------------------------------------------------------------
221 : ! ... merge the land component with the non-land component
222 : ! ocn and ice already have fractions factored in
223 : !-------------------------------------------------------------------------------------
224 0 : dvelocity(:ncol,spc_ndx(ispec)) = lnd(lchnk)%dvel(:ncol,ispec)*lndfrac(:ncol) &
225 0 : + ocnice_dvel(:ncol,spc_ndx(ispec))
226 : enddo
227 :
228 : !-------------------------------------------------------------------------------------
229 : ! ... special adjustments
230 : !-------------------------------------------------------------------------------------
231 0 : if( mpan_ndx>0 ) then
232 0 : dvelocity(:ncol,mpan_ndx) = dvelocity(:ncol,mpan_ndx)/3._r8
233 : endif
234 0 : if( xmpan_ndx>0 ) then
235 0 : dvelocity(:ncol,xmpan_ndx) = dvelocity(:ncol,xmpan_ndx)/3._r8
236 : endif
237 0 : if( hcn_ndx>0 ) then
238 0 : dvelocity(:ncol,hcn_ndx) = ocnice_dvel(:ncol,hcn_ndx) ! should be zero over land
239 : endif
240 0 : if( ch3cn_ndx>0 ) then
241 0 : dvelocity(:ncol,ch3cn_ndx) = ocnice_dvel(:ncol,ch3cn_ndx) ! should be zero over land
242 : endif
243 :
244 : ! HCOOH, use CH3COOH dep.vel
245 0 : if( hcooh_ndx > 0 .and. ch3cooh_ndx > 0 ) then
246 0 : if( has_dvel(hcooh_ndx) ) then
247 : dvelocity(:ncol,hcooh_ndx) = dvelocity(:ncol,ch3cooh_ndx)
248 : end if
249 : end if
250 :
251 : !-------------------------------------------------------------------------------------
252 : ! ... assign CO tags to CO
253 : ! put this kludge in for now ...
254 : ! -- should be able to set all these via the table mapping in shr_drydep_mod
255 : !-------------------------------------------------------------------------------------
256 0 : if( cohc_ndx>0 .and. co_ndx>0 ) then
257 0 : dvelocity(:ncol,cohc_ndx) = dvelocity(:ncol,co_ndx)
258 : dflx(:ncol,cohc_ndx) = dvelocity(:ncol,co_ndx) * term(:ncol) * mmr(:ncol,plev,cohc_ndx)
259 : endif
260 0 : if( come_ndx>0 .and. co_ndx>0 ) then
261 0 : dvelocity(:ncol,come_ndx) = dvelocity(:ncol,co_ndx)
262 : dflx(:ncol,come_ndx) = dvelocity(:ncol,co_ndx) * term(:ncol) * mmr(:ncol,plev,come_ndx)
263 : endif
264 :
265 0 : if ( co_ndx>0 ) then
266 0 : do i=1,tag_cnt
267 0 : dvelocity(:ncol,cotag_ndx(i)) = dvelocity(:ncol,co_ndx)
268 0 : dflx(:ncol,cotag_ndx(i)) = dvelocity(:ncol,co_ndx) * term(:ncol) * mmr(:ncol,plev,cotag_ndx(i))
269 : enddo
270 : endif
271 :
272 : do ispec = 1,nddvels
273 : !-------------------------------------------------------------------------------------
274 : ! ... compute the deposition flux
275 : !-------------------------------------------------------------------------------------
276 : dflx(:ncol,spc_ndx(ispec)) = dvelocity(:ncol,spc_ndx(ispec)) * term(:ncol) * mmr(:ncol,plev,spc_ndx(ispec))
277 : end do
278 :
279 0 : end subroutine drydep_fromlnd
280 :
281 : !-------------------------------------------------------------------------------------
282 : !-------------------------------------------------------------------------------------
283 0 : subroutine dvel_inti_xactive( depvel_lnd_file )
284 : !-------------------------------------------------------------------------------------
285 : ! ... intialize interactive drydep
286 : !-------------------------------------------------------------------------------------
287 : use dycore, only : dycore_is
288 : use mo_chem_utls, only : get_spc_ndx
289 : use phys_control, only : phys_getopts
290 :
291 : !-------------------------------------------------------------------------------------
292 : ! ... dummy arguments
293 : !-------------------------------------------------------------------------------------
294 : character(len=*), intent(in) :: depvel_lnd_file
295 :
296 : !-------------------------------------------------------------------------------------
297 : ! ... local variables
298 : !-------------------------------------------------------------------------------------
299 : integer :: i
300 : integer :: nlon_veg, nlat_veg, npft_veg
301 : integer :: dimid
302 : integer :: m
303 : integer :: astat
304 : integer :: plon, plat
305 : integer :: ierr, ndx
306 :
307 0 : real(r8), allocatable :: vegetation_map(:,:,:)
308 0 : real(r8), allocatable :: work(:,:)
309 0 : real(r8), allocatable :: landmask(:,:)
310 0 : real(r8), allocatable :: urban(:,:)
311 0 : real(r8), allocatable :: lake(:,:)
312 0 : real(r8), allocatable :: wetland(:,:)
313 0 : real(r8), allocatable :: lon_veg_edge(:)
314 0 : real(r8), allocatable :: lat_veg_edge(:)
315 :
316 : character(len=32) :: test_name
317 : character(len=4) :: tag_name
318 : type(file_desc_t) :: piofile
319 : type(var_desc_t) :: vid
320 :
321 : character(len=shr_kind_cl) :: locfn
322 : logical :: prog_modal_aero
323 :
324 : ! determine if modal aerosols are active so that fraction_landuse array is initialized for modal aerosal dry dep
325 0 : call phys_getopts(prog_modal_aero_out=prog_modal_aero)
326 :
327 0 : call dvel_inti_fromlnd()
328 :
329 0 : if( masterproc ) then
330 0 : write(iulog,*) 'drydep_inti: following species have dry deposition'
331 0 : do i=1,nddvels
332 0 : if( len_trim(drydep_list(i)) > 0 ) then
333 0 : write(iulog,*) 'drydep_inti: '//trim(drydep_list(i))//' is requested to have dry dep'
334 : endif
335 : enddo
336 0 : write(iulog,*) 'drydep_inti:'
337 : endif
338 :
339 : !-------------------------------------------------------------------------------------
340 : ! ... get species indices
341 : !-------------------------------------------------------------------------------------
342 0 : xpan_ndx = get_spc_ndx( 'XPAN' )
343 0 : xmpan_ndx = get_spc_ndx( 'XMPAN' )
344 0 : o3a_ndx = get_spc_ndx( 'O3A' )
345 :
346 0 : ch4_ndx = get_spc_ndx( 'CH4' )
347 0 : h2_ndx = get_spc_ndx( 'H2' )
348 0 : co_ndx = get_spc_ndx( 'CO' )
349 0 : pan_ndx = get_spc_ndx( 'PAN' )
350 0 : mpan_ndx = get_spc_ndx( 'MPAN' )
351 0 : o3_ndx = get_spc_ndx( 'OX' )
352 0 : if( o3_ndx < 0 ) then
353 0 : o3_ndx = get_spc_ndx( 'O3' )
354 : end if
355 0 : so2_ndx = get_spc_ndx( 'SO2' )
356 0 : ch3cooh_ndx = get_spc_ndx( 'CH3COOH')
357 :
358 0 : sogm_ndx = get_spc_ndx( 'SOGM' )
359 0 : sogi_ndx = get_spc_ndx( 'SOGI' )
360 0 : sogt_ndx = get_spc_ndx( 'SOGT' )
361 0 : sogb_ndx = get_spc_ndx( 'SOGB' )
362 0 : sogx_ndx = get_spc_ndx( 'SOGX' )
363 :
364 0 : hcn_ndx = get_spc_ndx( 'HCN')
365 0 : ch3cn_ndx = get_spc_ndx( 'CH3CN')
366 :
367 0 : cohc_ndx = get_spc_ndx( 'COhc' )
368 0 : come_ndx = get_spc_ndx( 'COme' )
369 :
370 0 : tag_cnt=0
371 0 : cotag_ndx(:)=-1
372 0 : do i = 1,NTAGS
373 0 : write(tag_name,'(a2,i2.2)') 'CO',i
374 0 : ndx = get_spc_ndx(tag_name)
375 0 : if (ndx>0) then
376 0 : tag_cnt = tag_cnt+1
377 0 : cotag_ndx(tag_cnt) = ndx
378 : endif
379 : enddo
380 :
381 0 : do i=1,nddvels
382 0 : if ( mapping(i) > 0 ) then
383 0 : test_name = drydep_list(i)
384 0 : m = get_spc_ndx( test_name )
385 0 : has_dvel(m) = .true.
386 : map_dvel(m) = i
387 : endif
388 : enddo
389 :
390 : if( all( .not. has_dvel(:) ) ) then
391 : return
392 : end if
393 :
394 : !---------------------------------------------------------------------------
395 : ! ... allocate module variables
396 : !---------------------------------------------------------------------------
397 : allocate( dep_ra(pcols,n_land_type,begchunk:endchunk),stat=astat )
398 : if( astat /= 0 ) then
399 : write(iulog,*) 'dvel_inti: failed to allocate dep_ra; error = ',astat
400 : call endrun('dvel_inti: failed to allocate dep_ra')
401 : end if
402 : allocate( dep_rb(pcols,n_land_type,begchunk:endchunk),stat=astat )
403 : if( astat /= 0 ) then
404 : write(iulog,*) 'dvel_inti: failed to allocate dep_rb; error = ',astat
405 : call endrun('dvel_inti: failed to allocate dep_rb')
406 : end if
407 :
408 : if (.not.prog_modal_aero) then
409 : return
410 : endif
411 :
412 : allocate( fraction_landuse(pcols,n_land_type, begchunk:endchunk),stat=astat )
413 : if( astat /= 0 ) then
414 : write(iulog,*) 'dvel_inti: failed to allocate fraction_landuse; error = ',astat
415 : call endrun('dvel_inti: failed to allocate fraction_landuse')
416 : end if
417 : fraction_landuse = nan
418 :
419 : plon = get_dyn_grid_parm('plon')
420 : plat = get_dyn_grid_parm('plat')
421 :
422 : if(dycore_is('UNSTRUCTURED') ) then
423 : call get_landuse_and_soilw_from_file()
424 : else
425 : !---------------------------------------------------------------------------
426 : ! ... read landuse map
427 : !---------------------------------------------------------------------------
428 : call getfil (depvel_lnd_file, locfn, 0)
429 : call cam_pio_openfile (piofile, trim(locfn), PIO_NOWRITE)
430 : !---------------------------------------------------------------------------
431 : ! ... get the dimensions
432 : !---------------------------------------------------------------------------
433 : ierr = pio_inq_dimid( piofile, 'lon', dimid )
434 : ierr = pio_inq_dimlen( piofile, dimid, nlon_veg )
435 : ierr = pio_inq_dimid( piofile, 'lat', dimid )
436 : ierr = pio_inq_dimlen( piofile, dimid, nlat_veg )
437 : ierr = pio_inq_dimid( piofile, 'pft', dimid )
438 : ierr = pio_inq_dimlen( piofile, dimid, npft_veg )
439 : !---------------------------------------------------------------------------
440 : ! ... allocate arrays
441 : !---------------------------------------------------------------------------
442 : allocate( vegetation_map(nlon_veg,nlat_veg,npft_veg), work(nlon_veg,nlat_veg), stat=astat )
443 : if( astat /= 0 ) then
444 : write(iulog,*) 'dvel_inti: failed to allocate vegetation_map; error = ',astat
445 : call endrun('dvel_inti: failed to allocate vegetation_map')
446 : end if
447 : allocate( urban(nlon_veg,nlat_veg), lake(nlon_veg,nlat_veg), &
448 : landmask(nlon_veg,nlat_veg), wetland(nlon_veg,nlat_veg), stat=astat )
449 : if( astat /= 0 ) then
450 : write(iulog,*) 'dvel_inti: failed to allocate vegetation_map; error = ',astat
451 : call endrun('dvel_inti: failed to allocate vegetation_map')
452 : end if
453 : allocate( lon_veg_edge(nlon_veg+1), lat_veg_edge(nlat_veg+1), stat=astat )
454 : if( astat /= 0 ) then
455 : write(iulog,*) 'dvel_inti: failed to allocate vegetation lon, lat arrays; error = ',astat
456 : call endrun('dvel_inti: failed to allocate vegetation lon, lat arrays')
457 : end if
458 : !---------------------------------------------------------------------------
459 : ! ... read the vegetation map and landmask
460 : !---------------------------------------------------------------------------
461 : ierr = pio_inq_varid( piofile, 'PCT_PFT', vid )
462 : ierr = pio_get_var( piofile, vid, vegetation_map )
463 :
464 : ierr = pio_inq_varid( piofile, 'LANDMASK', vid )
465 : ierr = pio_get_var( piofile, vid, landmask )
466 :
467 : ierr = pio_inq_varid( piofile, 'PCT_URBAN', vid )
468 : ierr = pio_get_var( piofile, vid, urban )
469 :
470 : ierr = pio_inq_varid( piofile, 'PCT_LAKE', vid )
471 : ierr = pio_get_var( piofile, vid, lake )
472 :
473 : ierr = pio_inq_varid( piofile, 'PCT_WETLAND', vid )
474 : ierr = pio_get_var( piofile, vid, wetland )
475 :
476 : call cam_pio_closefile( piofile )
477 :
478 : !---------------------------------------------------------------------------
479 : ! scale vegetation, urban, lake, and wetland to fraction
480 : !---------------------------------------------------------------------------
481 : vegetation_map(:,:,:) = .01_r8 * vegetation_map(:,:,:)
482 : wetland(:,:) = .01_r8 * wetland(:,:)
483 : lake(:,:) = .01_r8 * lake(:,:)
484 : urban(:,:) = .01_r8 * urban(:,:)
485 : #ifdef DEBUG
486 : if(masterproc) then
487 : write(iulog,*) 'minmax vegetation_map ',minval(vegetation_map),maxval(vegetation_map)
488 : write(iulog,*) 'minmax wetland ',minval(wetland),maxval(wetland)
489 : write(iulog,*) 'minmax landmask ',minval(landmask),maxval(landmask)
490 : end if
491 : #endif
492 : !---------------------------------------------------------------------------
493 : ! ... define lat-lon of vegetation map (1x1)
494 : !---------------------------------------------------------------------------
495 : lat_veg_edge(:) = (/ (-90.0_r8 + (i-1),i=1,nlat_veg+1) /)
496 : lon_veg_edge(:) = (/ ( 0.0_r8 + (i-1),i=1,nlon_veg+1) /)
497 :
498 : !---------------------------------------------------------------------------
499 : ! ... regrid to model grid
500 : !---------------------------------------------------------------------------
501 : call interp_map( plon, plat, nlon_veg, nlat_veg, npft_veg, lat_veg_edge, &
502 : lon_veg_edge, landmask, urban, lake, &
503 : wetland, vegetation_map )
504 :
505 : deallocate( vegetation_map, work, stat=astat )
506 : deallocate( lon_veg_edge, lat_veg_edge, stat=astat )
507 : deallocate( landmask, urban, lake, wetland, stat=astat )
508 : endif ! Unstructured grid
509 :
510 0 : end subroutine dvel_inti_xactive
511 :
512 : !-------------------------------------------------------------------------------------
513 : subroutine get_landuse_and_soilw_from_file()
514 : use ncdio_atm, only : infld
515 :
516 : logical :: readvar
517 :
518 : type(file_desc_t) :: piofile
519 : character(len=shr_kind_cl) :: locfn
520 : logical :: lexist
521 :
522 : if (len_trim(drydep_srf_file) == 0) then
523 : write(iulog,*)'**************************************'
524 : write(iulog,*)' get_landuse_and_soilw_from_file: INFO:'
525 : write(iulog,*)' drydep_srf_file not set:'
526 : write(iulog,*)' setting fraction_landuse to zero'
527 : write(iulog,*)'**************************************'
528 : fraction_landuse = 0._r8
529 : return
530 : end if
531 :
532 : call getfil (drydep_srf_file, locfn, 1, lexist)
533 : if(lexist) then
534 : call cam_pio_openfile(piofile, locfn, PIO_NOWRITE)
535 :
536 : call infld('fraction_landuse', piofile, 'ncol','class',1,pcols,1,n_land_type, begchunk,endchunk, &
537 : fraction_landuse, readvar, gridname='physgrid')
538 : if (.not. readvar) then
539 : write(iulog,*)'**************************************'
540 : write(iulog,*)'get_landuse_and_soilw_from_file: INFO:'
541 : write(iulog,*)' fraction_landuse not read from file: '
542 : write(iulog,*)' ', trim(locfn)
543 : write(iulog,*)' setting all values to zero'
544 : write(iulog,*)'**************************************'
545 : fraction_landuse = 0._r8
546 : end if
547 :
548 : call cam_pio_closefile(piofile)
549 : else
550 : call endrun('Unstructured grids require drydep_srf_file ')
551 : end if
552 :
553 :
554 : end subroutine get_landuse_and_soilw_from_file
555 :
556 : !-------------------------------------------------------------------------------------
557 : subroutine interp_map( plon, plat, nlon_veg, nlat_veg, npft_veg, lat_veg_edge, &
558 : lon_veg_edge, landmask, urban, lake, &
559 : wetland, vegetation_map )
560 :
561 : use mo_constants, only : r2d
562 : use scamMod, only : latiop,loniop,scmlat,scmlon,scm_cambfb_mode
563 : use shr_scam_mod , only: shr_scam_getCloseLatLon ! Standardized system subroutines
564 : use cam_initfiles, only: initial_file_get_id
565 : use dycore, only : dycore_is
566 : use phys_grid, only : get_rlat_all_p, get_rlon_all_p, get_ncols_p
567 :
568 : !-------------------------------------------------------------------------------------
569 : ! ... dummy arguments
570 : !-------------------------------------------------------------------------------------
571 : integer, intent(in) :: plon, plat, nlon_veg, nlat_veg, npft_veg
572 : real(r8), intent(in) :: landmask(nlon_veg,nlat_veg)
573 : real(r8), intent(in) :: urban(nlon_veg,nlat_veg)
574 : real(r8), intent(in) :: lake(nlon_veg,nlat_veg)
575 : real(r8), intent(in) :: wetland(nlon_veg,nlat_veg)
576 : real(r8), intent(in) :: vegetation_map(nlon_veg,nlat_veg,npft_veg)
577 : real(r8), intent(in) :: lon_veg_edge(nlon_veg+1)
578 : real(r8), intent(in) :: lat_veg_edge(nlat_veg+1)
579 :
580 : !-------------------------------------------------------------------------------------
581 : ! ... local variables
582 : !-------------------------------------------------------------------------------------
583 : real(r8) :: closelat,closelon
584 : integer :: latidx,lonidx
585 :
586 : integer, parameter :: veg_ext = 20
587 : type(file_desc_t), pointer :: piofile
588 : integer :: i, j, ii, jj, i_ndx, n
589 : integer, dimension(plon+1) :: ind_lon
590 : integer, dimension(plat+1) :: ind_lat
591 : real(r8) :: total_land
592 : real(r8), dimension(plon+1) :: lon_edge
593 : real(r8), dimension(plat+1) :: lat_edge
594 : real(r8) :: lat1, lon1
595 : real(r8) :: x1, x2, y1, y2, dx, dy
596 : real(r8) :: area, total_area
597 : real(r8), dimension(npft_veg+3) :: fraction
598 : real(r8), dimension(-veg_ext:nlon_veg+veg_ext) :: lon_veg_edge_ext
599 : integer, dimension(-veg_ext:nlon_veg+veg_ext) :: mapping_ext
600 :
601 : real(r8), allocatable :: lam(:), phi(:)
602 :
603 : logical, parameter :: has_npole = .true.
604 : integer :: ploniop,platiop
605 : real(r8) :: tmp_frac_lu(plon,n_land_type,plat)
606 :
607 : real(r8):: rlats(pcols), rlons(pcols)
608 : integer :: lchnk, ncol, icol
609 : logical :: found
610 :
611 : if(dycore_is('UNSTRUCTURED') ) then
612 : call endrun('mo_drydep::interp_map called for UNSTRUCTURED grid')
613 : endif
614 :
615 : allocate(lam(plon), phi(plat))
616 : call get_horiz_grid_d(plat, clat_d_out=phi)
617 : call get_horiz_grid_d(plon, clon_d_out=lam)
618 :
619 : if (single_column) then
620 : if (scm_cambfb_mode) then
621 : piofile => initial_file_get_id()
622 : call shr_scam_getCloseLatLon(piofile,scmlat,scmlon,closelat,closelon,latidx,lonidx)
623 : ploniop=size(loniop)
624 : platiop=size(latiop)
625 : else
626 : latidx=1
627 : lonidx=1
628 : ploniop=1
629 : platiop=1
630 : end if
631 :
632 : lon_edge(1) = loniop(lonidx) * r2d - .5_r8*(loniop(2) - loniop(1)) * r2d
633 :
634 : if (lonidx.lt.ploniop) then
635 : lon_edge(2) = loniop(lonidx+1) * r2d - .5_r8*(loniop(2) - loniop(1)) * r2d
636 : else
637 : lon_edge(2) = lon_edge(1) + (loniop(2) - loniop(1)) * r2d
638 : end if
639 :
640 : lat_edge(1) = latiop(latidx) * r2d - .5_r8*(latiop(2) - latiop(1)) * r2d
641 :
642 : if (latidx.lt.platiop) then
643 : lat_edge(2) = latiop(latidx+1) * r2d - .5_r8*(latiop(2) - latiop(1)) * r2d
644 : else
645 : lat_edge(2) = lat_edge(1) + (latiop(2) - latiop(1)) * r2d
646 : end if
647 : else
648 : do i = 1,plon
649 : lon_edge(i) = lam(i) * r2d - .5_r8*(lam(2) - lam(1)) * r2d
650 : end do
651 : lon_edge(plon+1) = lon_edge(plon) + (lam(2) - lam(1)) * r2d
652 : if( .not. has_npole ) then
653 : do j = 1,plat+1
654 : lat_edge(j) = phi(j) * r2d - .5_r8*(phi(2) - phi(1)) * r2d
655 : end do
656 : else
657 : do j = 1,plat
658 : lat_edge(j) = phi(j) * r2d - .5_r8*(phi(2) - phi(1)) * r2d
659 : end do
660 : lat_edge(plat+1) = lat_edge(plat) + (phi(2) - phi(1)) * r2d
661 : end if
662 : end if
663 : do j = 1,plat+1
664 : lat_edge(j) = min( lat_edge(j), 90._r8 )
665 : lat_edge(j) = max( lat_edge(j),-90._r8 )
666 : end do
667 :
668 : !-------------------------------------------------------------------------------------
669 : ! wrap around the longitudes
670 : !-------------------------------------------------------------------------------------
671 : do i = -veg_ext,0
672 : lon_veg_edge_ext(i) = lon_veg_edge(nlon_veg+i) - 360._r8
673 : mapping_ext (i) = nlon_veg+i
674 : end do
675 : do i = 1,nlon_veg
676 : lon_veg_edge_ext(i) = lon_veg_edge(i)
677 : mapping_ext (i) = i
678 : end do
679 : do i = nlon_veg+1,nlon_veg+veg_ext
680 : lon_veg_edge_ext(i) = lon_veg_edge(i-nlon_veg) + 360._r8
681 : mapping_ext (i) = i-nlon_veg
682 : end do
683 : #ifdef DEBUG
684 : write(iulog,*) 'interp_map : lon_edge ',lon_edge
685 : write(iulog,*) 'interp_map : lat_edge ',lat_edge
686 : write(iulog,*) 'interp_map : mapping_ext ',mapping_ext
687 : #endif
688 : do j = 1,plon+1
689 : lon1 = lon_edge(j)
690 : do i = -veg_ext,nlon_veg+veg_ext
691 : dx = lon_veg_edge_ext(i ) - lon1
692 : dy = lon_veg_edge_ext(i+1) - lon1
693 : if( dx*dy <= 0._r8 ) then
694 : ind_lon(j) = i
695 : exit
696 : end if
697 : end do
698 : end do
699 :
700 : do j = 1,plat+1
701 : lat1 = lat_edge(j)
702 : do i = 1,nlat_veg
703 : dx = lat_veg_edge(i ) - lat1
704 : dy = lat_veg_edge(i+1) - lat1
705 : if( dx*dy <= 0._r8 ) then
706 : ind_lat(j) = i
707 : exit
708 : end if
709 : end do
710 : end do
711 : #ifdef DEBUG
712 : write(iulog,*) 'interp_map : ind_lon ',ind_lon
713 : write(iulog,*) 'interp_map : ind_lat ',ind_lat
714 : #endif
715 : lat_loop : do j = 1,plat
716 : lon_loop : do i = 1,plon
717 : total_area = 0._r8
718 : fraction = 0._r8
719 : do jj = ind_lat(j),ind_lat(j+1)
720 : y1 = max( lat_edge(j),lat_veg_edge(jj) )
721 : y2 = min( lat_edge(j+1),lat_veg_edge(jj+1) )
722 : dy = (y2 - y1)/(lat_veg_edge(jj+1) - lat_veg_edge(jj))
723 : do ii =ind_lon(i),ind_lon(i+1)
724 : i_ndx = mapping_ext(ii)
725 : x1 = max( lon_edge(i),lon_veg_edge_ext(ii) )
726 : x2 = min( lon_edge(i+1),lon_veg_edge_ext(ii+1) )
727 : dx = (x2 - x1)/(lon_veg_edge_ext(ii+1) - lon_veg_edge_ext(ii))
728 : area = dx * dy
729 : total_area = total_area + area
730 : !-----------------------------------------------------------------
731 : ! ... special case for ocean grid point
732 : !-----------------------------------------------------------------
733 : if( nint(landmask(i_ndx,jj)) == 0 ) then
734 : fraction(npft_veg+1) = fraction(npft_veg+1) + area
735 : else
736 : do n = 1,npft_veg
737 : fraction(n) = fraction(n) + vegetation_map(i_ndx,jj,n) * area
738 : end do
739 : fraction(npft_veg+1) = fraction(npft_veg+1) + area * lake (i_ndx,jj)
740 : fraction(npft_veg+2) = fraction(npft_veg+2) + area * wetland(i_ndx,jj)
741 : fraction(npft_veg+3) = fraction(npft_veg+3) + area * urban (i_ndx,jj)
742 : !-----------------------------------------------------------------
743 : ! ... check if land accounts for the whole area.
744 : ! If not, the remaining area is in the ocean
745 : !-----------------------------------------------------------------
746 : total_land = sum(vegetation_map(i_ndx,jj,:)) &
747 : + urban (i_ndx,jj) &
748 : + lake (i_ndx,jj) &
749 : + wetland(i_ndx,jj)
750 : if( total_land < 1._r8 ) then
751 : fraction(npft_veg+1) = fraction(npft_veg+1) + (1._r8 - total_land) * area
752 : end if
753 : end if
754 : end do
755 : end do
756 : !-------------------------------------------------------------------------------------
757 : ! ... divide by total area of grid box
758 : !-------------------------------------------------------------------------------------
759 : fraction(:) = fraction(:)/total_area
760 : !-------------------------------------------------------------------------------------
761 : ! ... make sure we don't have too much or too little
762 : !-------------------------------------------------------------------------------------
763 : if( abs( sum(fraction) - 1._r8) > .001_r8 ) then
764 : fraction(:) = fraction(:)/sum(fraction)
765 : end if
766 : !-------------------------------------------------------------------------------------
767 : ! ... map to Wesely land classification
768 : !-------------------------------------------------------------------------------------
769 : tmp_frac_lu(i, 1, j) = fraction(20)
770 : tmp_frac_lu(i, 2, j) = sum(fraction(16:17))
771 : tmp_frac_lu(i, 3, j) = sum(fraction(13:15))
772 : tmp_frac_lu(i, 4, j) = sum(fraction( 5: 9))
773 : tmp_frac_lu(i, 5, j) = sum(fraction( 2: 4))
774 : tmp_frac_lu(i, 6, j) = fraction(19)
775 : tmp_frac_lu(i, 7, j) = fraction(18)
776 : tmp_frac_lu(i, 8, j) = fraction( 1)
777 : tmp_frac_lu(i, 9, j) = 0._r8
778 : tmp_frac_lu(i,10, j) = 0._r8
779 : tmp_frac_lu(i,11, j) = sum(fraction(10:12))
780 : end do lon_loop
781 : end do lat_loop
782 :
783 : do lchnk = begchunk, endchunk
784 : ncol = get_ncols_p(lchnk)
785 : call get_rlat_all_p(lchnk, ncol, rlats(:ncol))
786 : call get_rlon_all_p(lchnk, ncol, rlons(:ncol))
787 : do icol= 1,ncol
788 : found=.false.
789 : find_col: do j = 1,plat
790 : do i = 1,plon
791 : if (rlats(icol)==phi(j) .and. rlons(icol)==lam(i)) then
792 : found=.true.
793 : exit find_col
794 : endif
795 : enddo
796 : enddo find_col
797 :
798 : if (.not.found) call endrun('mo_drydep::interp_map not able find physics column coordinate')
799 : fraction_landuse(icol,1:n_land_type,lchnk) = tmp_frac_lu(i,1:n_land_type,j)
800 :
801 : end do
802 :
803 : !-------------------------------------------------------------------------------------
804 : ! ... make sure there are no out of range values
805 : !-------------------------------------------------------------------------------------
806 : where (fraction_landuse(:ncol,:n_land_type,lchnk) < 0._r8) fraction_landuse(:ncol,:n_land_type,lchnk) = 0._r8
807 : where (fraction_landuse(:ncol,:n_land_type,lchnk) > 1._r8) fraction_landuse(:ncol,:n_land_type,lchnk) = 1._r8
808 : end do
809 :
810 : end subroutine interp_map
811 :
812 : !-------------------------------------------------------------------------------------
813 : !-------------------------------------------------------------------------------------
814 0 : subroutine drydep_xactive( sfc_temp, pressure_sfc, &
815 : wind_speed, spec_hum, air_temp, pressure_10m, rain, &
816 0 : snow, solar_flux, dvel, dflx, mmr, &
817 : tv, ncol, lchnk, &
818 : ocnfrc, icefrc, beglandtype, endlandtype )
819 : !-------------------------------------------------------------------------------------
820 : ! code based on wesely (atmospheric environment, 1989, vol 23, p. 1293-1304) for
821 : ! calculation of r_c, and on walcek et. al. (atmospheric enviroment, 1986,
822 : ! vol. 20, p. 949-964) for calculation of r_a and r_b
823 : !
824 : ! as suggested in walcek (u_i)(u*_i) = (u_a)(u*_a)
825 : ! is kept constant where i represents a subgrid environment and a the
826 : ! grid average environment. thus the calculation proceeds as follows:
827 : ! va the grid averaged wind is calculated on dots
828 : ! z0(i) the grid averaged roughness coefficient is calculated
829 : ! ri(i) the grid averaged richardson number is calculated
830 : ! --> the grid averaged (u_a)(u*_a) is calculated
831 : ! --> subgrid scale u*_i is calculated assuming (u_i) given as above
832 : ! --> final deposotion velocity is weighted average of subgrid scale velocities
833 : !
834 : ! code written by P. Hess, rewritten in fortran 90 by JFL (August 2000)
835 : ! modified by JFL to be used in MOZART-2 (October 2002)
836 : !-------------------------------------------------------------------------------------
837 :
838 : use shr_drydep_mod, only: z0, rgso, rgss, ri, rclo, rcls, rlu, rac
839 : use shr_drydep_mod, only: shr_drydep_setHCoeff, foxd, drat
840 : use physconst, only: tmelt
841 :
842 : !-------------------------------------------------------------------------------------
843 : ! ... dummy arguments
844 : !-------------------------------------------------------------------------------------
845 : integer, intent(in) :: ncol
846 : real(r8), intent(in) :: sfc_temp(pcols) ! surface temperature (K)
847 : real(r8), intent(in) :: pressure_sfc(pcols) ! surface pressure (Pa)
848 : real(r8), intent(in) :: wind_speed(pcols) ! 10 meter wind speed (m/s)
849 : real(r8), intent(in) :: spec_hum(pcols) ! specific humidity (kg/kg)
850 : real(r8), intent(in) :: air_temp(pcols) ! surface air temperature (K)
851 : real(r8), intent(in) :: pressure_10m(pcols) ! 10 meter pressure (Pa)
852 : real(r8), intent(in) :: rain(pcols)
853 : real(r8), intent(in) :: snow(pcols) ! snow height (m)
854 :
855 : real(r8), intent(in) :: solar_flux(pcols) ! direct shortwave radiation at surface (W/m^2)
856 : real(r8), intent(in) :: tv(pcols) ! potential temperature
857 : real(r8), intent(in) :: mmr(pcols,plev,gas_pcnst) ! constituent concentration (kg/kg)
858 : real(r8), intent(out) :: dvel(ncol,gas_pcnst) ! deposition velocity (cm/s)
859 : real(r8), intent(inout) :: dflx(pcols,gas_pcnst) ! deposition flux (/cm^2/s)
860 :
861 : integer, intent(in) :: lchnk ! chunk number
862 :
863 : integer, intent(in), optional :: beglandtype
864 : integer, intent(in), optional :: endlandtype
865 :
866 : real(r8), intent(in), optional :: ocnfrc(pcols)
867 : real(r8), intent(in), optional :: icefrc(pcols)
868 :
869 : !-------------------------------------------------------------------------------------
870 : ! ... local variables
871 : !-------------------------------------------------------------------------------------
872 : real(r8), parameter :: scaling_to_cm_per_s = 100._r8
873 : real(r8), parameter :: rain_threshold = 1.e-7_r8 ! of the order of 1cm/day expressed in m/s
874 :
875 : integer :: i, ispec, lt, m
876 : integer :: sndx
877 :
878 : real(r8) :: slope = 0._r8
879 : real(r8) :: z0water ! revised z0 over water
880 : real(r8) :: p ! pressure at midpoint first layer
881 : real(r8) :: pg ! surface pressure
882 : real(r8) :: es ! saturation vapor pressure
883 : real(r8) :: ws ! saturation mixing ratio
884 : real(r8) :: hvar ! constant to compute xmol
885 : real(r8) :: h ! constant to compute xmol
886 : real(r8) :: psih ! stability correction factor
887 : real(r8) :: rs ! constant for calculating rsmx
888 : real(r8) :: rmx ! resistance by vegetation
889 : real(r8) :: zovl ! ratio of z to m-o length
890 : real(r8) :: cvarb ! cvar averaged over landtypes
891 : real(r8) :: bb ! b averaged over landtypes
892 : real(r8) :: ustarb ! ustar averaged over landtypes
893 0 : real(r8) :: tc(ncol) ! temperature in celsius
894 0 : real(r8) :: cts(ncol) ! correction to rlu rcl and rgs for frost
895 :
896 : !-------------------------------------------------------------------------------------
897 : ! local arrays: dependent on location and species
898 : !-------------------------------------------------------------------------------------
899 0 : real(r8), dimension(ncol,nddvels) :: heff
900 :
901 : !-------------------------------------------------------------------------------------
902 : ! local arrays: dependent on location only
903 : !-------------------------------------------------------------------------------------
904 0 : integer :: index_season(ncol,n_land_type)
905 0 : real(r8), dimension(ncol) :: tha ! atmospheric virtual potential temperature
906 0 : real(r8), dimension(ncol) :: thg ! ground virtual potential temperature
907 0 : real(r8), dimension(ncol) :: z ! height of lowest level
908 0 : real(r8), dimension(ncol) :: va ! magnitude of v on cross points
909 0 : real(r8), dimension(ncol) :: ribn ! richardson number
910 0 : real(r8), dimension(ncol) :: qs ! saturation specific humidity
911 0 : real(r8), dimension(ncol) :: crs ! multiplier to calculate crs
912 0 : real(r8), dimension(ncol) :: rdc ! part of lower canopy resistance
913 0 : real(r8), dimension(ncol) :: uustar ! u*ustar (assumed constant over grid)
914 0 : real(r8), dimension(ncol) :: z0b ! average roughness length over grid
915 0 : real(r8), dimension(ncol) :: wrk ! work array
916 0 : real(r8), dimension(ncol) :: term ! work array
917 0 : real(r8), dimension(ncol) :: resc ! work array
918 0 : real(r8), dimension(ncol) :: lnd_frc ! work array
919 0 : logical, dimension(ncol) :: unstable
920 0 : logical, dimension(ncol) :: has_rain
921 0 : logical, dimension(ncol) :: has_dew
922 :
923 : !-------------------------------------------------------------------------------------
924 : ! local arrays: dependent on location and landtype
925 : !-------------------------------------------------------------------------------------
926 0 : real(r8), dimension(ncol,n_land_type) :: rds ! resistance for deposition of sulfate
927 0 : real(r8), dimension(ncol,n_land_type) :: b ! buoyancy parameter for unstable conditions
928 0 : real(r8), dimension(ncol,n_land_type) :: cvar ! height parameter
929 0 : real(r8), dimension(ncol,n_land_type) :: ustar ! friction velocity
930 0 : real(r8), dimension(ncol,n_land_type) :: xmol ! monin-obukhov length
931 :
932 : !-------------------------------------------------------------------------------------
933 : ! local arrays: dependent on location, landtype and species
934 : !-------------------------------------------------------------------------------------
935 0 : real(r8), dimension(ncol,n_land_type,gas_pcnst) :: rsmx ! vegetative resistance (plant mesophyll)
936 0 : real(r8), dimension(ncol,n_land_type,gas_pcnst) :: rclx ! lower canopy resistance
937 0 : real(r8), dimension(ncol,n_land_type,gas_pcnst) :: rlux ! vegetative resistance (upper canopy)
938 0 : real(r8), dimension(ncol,n_land_type) :: rlux_o3 ! vegetative resistance (upper canopy)
939 0 : real(r8), dimension(ncol,n_land_type,gas_pcnst) :: rgsx ! ground resistance
940 : real(r8) :: vds
941 0 : logical :: fr_lnduse(ncol,n_land_type) ! wrking array
942 : real(r8) :: dewm ! multiplier for rs when dew occurs
943 :
944 0 : real(r8) :: lcl_frc_landuse(ncol,n_land_type)
945 :
946 : integer :: beglt, endlt
947 :
948 : !-------------------------------------------------------------------------------------
949 : ! jfl : mods for PAN
950 : !-------------------------------------------------------------------------------------
951 : real(r8) :: dv_pan
952 : real(r8) :: c0_pan(11) = (/ 0.000_r8, 0.006_r8, 0.002_r8, 0.009_r8, 0.015_r8, &
953 : 0.006_r8, 0.000_r8, 0.000_r8, 0.000_r8, 0.002_r8, 0.002_r8 /)
954 : real(r8) :: k_pan (11) = (/ 0.000_r8, 0.010_r8, 0.005_r8, 0.004_r8, 0.003_r8, &
955 : 0.005_r8, 0.000_r8, 0.000_r8, 0.000_r8, 0.075_r8, 0.002_r8 /)
956 :
957 : if (present( beglandtype)) then
958 0 : beglt = beglandtype
959 : else
960 : beglt = 1
961 : endif
962 : if (present( endlandtype)) then
963 0 : endlt = endlandtype
964 : else
965 : endlt = n_land_type
966 : endif
967 :
968 : !-------------------------------------------------------------------------------------
969 : ! initialize
970 : !-------------------------------------------------------------------------------------
971 : do m = 1,gas_pcnst
972 : dvel(:,m) = 0._r8
973 : end do
974 :
975 : if( all( .not. has_dvel(:) ) ) then
976 : return
977 : end if
978 :
979 : !-------------------------------------------------------------------------------------
980 : ! define species-dependent parameters (temperature dependent)
981 : !-------------------------------------------------------------------------------------
982 : call shr_drydep_setHCoeff( ncol, sfc_temp, heff )
983 :
984 : do lt = 1,n_land_type
985 : dep_ra (:,lt,lchnk) = 0._r8
986 : dep_rb (:,lt,lchnk) = 0._r8
987 : rds(:,lt) = 0._r8
988 : end do
989 :
990 : !-------------------------------------------------------------------------------------
991 : ! season index only for ocn and sea ice
992 : !-------------------------------------------------------------------------------------
993 : index_season = 4
994 : !-------------------------------------------------------------------------------------
995 : ! special case for snow covered terrain
996 : !-------------------------------------------------------------------------------------
997 : do i = 1,ncol
998 : if( snow(i) > .01_r8 ) then
999 : index_season(i,:) = 4
1000 : end if
1001 : end do
1002 : !-------------------------------------------------------------------------------------
1003 : ! scale rain and define logical arrays
1004 : !-------------------------------------------------------------------------------------
1005 : has_rain(:ncol) = rain(:ncol) > rain_threshold
1006 :
1007 : !-------------------------------------------------------------------------------------
1008 : ! loop over longitude points
1009 : !-------------------------------------------------------------------------------------
1010 : col_loop : do i = 1,ncol
1011 : p = pressure_10m(i)
1012 : pg = pressure_sfc(i)
1013 : !-------------------------------------------------------------------------------------
1014 : ! potential temperature
1015 : !-------------------------------------------------------------------------------------
1016 : tha(i) = air_temp(i) * (p00/p )**rovcp * (1._r8 + .61_r8*spec_hum(i))
1017 : thg(i) = sfc_temp(i) * (p00/pg)**rovcp * (1._r8 + .61_r8*spec_hum(i))
1018 : !-------------------------------------------------------------------------------------
1019 : ! height of 1st level
1020 : !-------------------------------------------------------------------------------------
1021 : z(i) = - r/grav * air_temp(i) * (1._r8 + .61_r8*spec_hum(i)) * log(p/pg)
1022 : !-------------------------------------------------------------------------------------
1023 : ! wind speed
1024 : !-------------------------------------------------------------------------------------
1025 : va(i) = max( .01_r8,wind_speed(i) )
1026 : !-------------------------------------------------------------------------------------
1027 : ! Richardson number
1028 : !-------------------------------------------------------------------------------------
1029 : ribn(i) = z(i) * grav * (tha(i) - thg(i))/thg(i) / (va(i)*va(i))
1030 : ribn(i) = min( ribn(i),ric )
1031 : unstable(i) = ribn(i) < 0._r8
1032 : !-------------------------------------------------------------------------------------
1033 : ! saturation vapor pressure (Pascals)
1034 : ! saturation mixing ratio
1035 : ! saturation specific humidity
1036 : !-------------------------------------------------------------------------------------
1037 : es = 611._r8*exp( 5414.77_r8*(sfc_temp(i) - tmelt)/(tmelt*sfc_temp(i)) )
1038 : ws = .622_r8*es/(pg - es)
1039 : qs(i) = ws/(1._r8 + ws)
1040 : has_dew(i) = .false.
1041 : if( qs(i) <= spec_hum(i) ) then
1042 : has_dew(i) = .true.
1043 : end if
1044 : if( sfc_temp(i) < tmelt ) then
1045 : has_dew(i) = .false.
1046 : end if
1047 : !-------------------------------------------------------------------------------------
1048 : ! constant in determining rs
1049 : !-------------------------------------------------------------------------------------
1050 : tc(i) = sfc_temp(i) - tmelt
1051 : if( sfc_temp(i) > tmelt .and. sfc_temp(i) < 313.15_r8 ) then
1052 : crs(i) = (1._r8 + (200._r8/(solar_flux(i) + .1_r8))**2) * (400._r8/(tc(i)*(40._r8 - tc(i))))
1053 : else
1054 : crs(i) = large_value
1055 : end if
1056 : !-------------------------------------------------------------------------------------
1057 : ! rdc (lower canopy res)
1058 : !-------------------------------------------------------------------------------------
1059 : rdc(i) = 100._r8*(1._r8 + 1000._r8/(solar_flux(i) + 10._r8))/(1._r8 + 1000._r8*slope)
1060 : end do col_loop
1061 :
1062 : !-------------------------------------------------------------------------------------
1063 : ! ... form working arrays
1064 : !-------------------------------------------------------------------------------------
1065 : lcl_frc_landuse(:,:) = 0._r8
1066 :
1067 : if ( present(ocnfrc) .and. present(icefrc) ) then
1068 : do i=1,ncol
1069 : ! land type 7 is used for ocean
1070 : ! land type 8 is used for sea ice
1071 : lcl_frc_landuse(i,7) = ocnfrc(i)
1072 : lcl_frc_landuse(i,8) = icefrc(i)
1073 : enddo
1074 : endif
1075 : do lt = 1,n_land_type
1076 : do i=1,ncol
1077 : fr_lnduse(i,lt) = lcl_frc_landuse(i,lt) > 0._r8
1078 : enddo
1079 : end do
1080 :
1081 : !-------------------------------------------------------------------------------------
1082 : ! find grid averaged z0: z0bar (the roughness length) z_o=exp[S(f_i*ln(z_oi))]
1083 : ! this is calculated so as to find u_i, assuming u*u=u_i*u_i
1084 : !-------------------------------------------------------------------------------------
1085 : z0b(:) = 0._r8
1086 : do lt = 1,n_land_type
1087 : do i = 1,ncol
1088 : if( fr_lnduse(i,lt) ) then
1089 : z0b(i) = z0b(i) + lcl_frc_landuse(i,lt) * log( z0(index_season(i,lt),lt) )
1090 : end if
1091 : end do
1092 : end do
1093 :
1094 : !-------------------------------------------------------------------------------------
1095 : ! find the constant velocity uu*=(u_i)(u*_i)
1096 : !-------------------------------------------------------------------------------------
1097 : do i = 1,ncol
1098 : z0b(i) = exp( z0b(i) )
1099 : cvarb = vonkar/log( z(i)/z0b(i) )
1100 : !-------------------------------------------------------------------------------------
1101 : ! unstable and stable cases
1102 : !-------------------------------------------------------------------------------------
1103 : if( unstable(i) ) then
1104 : bb = 9.4_r8*(cvarb**2)*sqrt( abs(ribn(i))*z(i)/z0b(i) )
1105 : ustarb = cvarb * va(i) * sqrt( 1._r8 - (9.4_r8*ribn(i)/(1._r8 + 7.4_r8*bb)) )
1106 : else
1107 : ustarb = cvarb * va(i)/(1._r8 + 4.7_r8*ribn(i))
1108 : end if
1109 : uustar(i) = va(i)*ustarb
1110 : end do
1111 :
1112 : !-------------------------------------------------------------------------------------
1113 : ! calculate the friction velocity for each land type u_i=uustar/u*_i
1114 : !-------------------------------------------------------------------------------------
1115 : do lt = beglt,endlt
1116 : do i = 1,ncol
1117 : if( fr_lnduse(i,lt) ) then
1118 : if( unstable(i) ) then
1119 : cvar(i,lt) = vonkar/log( z(i)/z0(index_season(i,lt),lt) )
1120 : b(i,lt) = 9.4_r8*(cvar(i,lt)**2)* sqrt( abs(ribn(i))*z(i)/z0(index_season(i,lt),lt) )
1121 : ustar(i,lt) = sqrt( cvar(i,lt)*uustar(i)*sqrt( 1._r8 - (9.4_r8*ribn(i)/(1._r8 + 7.4_r8*b(i,lt))) ) )
1122 : else
1123 : cvar(i,lt) = vonkar/log( z(i)/z0(index_season(i,lt),lt) )
1124 : ustar(i,lt) = sqrt( cvar(i,lt)*uustar(i)/(1._r8 + 4.7_r8*ribn(i)) )
1125 : end if
1126 : end if
1127 : end do
1128 : end do
1129 :
1130 : !-------------------------------------------------------------------------------------
1131 : ! revise calculation of friction velocity and z0 over water
1132 : !-------------------------------------------------------------------------------------
1133 : lt = 7
1134 : do i = 1,ncol
1135 : if( fr_lnduse(i,lt) ) then
1136 : if( unstable(i) ) then
1137 : z0water = (.016_r8*(ustar(i,lt)**2)/grav) + diffk/(9.1_r8*ustar(i,lt))
1138 : cvar(i,lt) = vonkar/(log( z(i)/z0water ))
1139 : b(i,lt) = 9.4_r8*(cvar(i,lt)**2)*sqrt( abs(ribn(i))*z(i)/z0water )
1140 : ustar(i,lt) = sqrt( cvar(i,lt)*uustar(i)* sqrt( 1._r8 - (9.4_r8*ribn(i)/(1._r8+ 7.4_r8*b(i,lt))) ) )
1141 : else
1142 : z0water = (.016_r8*(ustar(i,lt)**2)/grav) + diffk/(9.1_r8*ustar(i,lt))
1143 : cvar(i,lt) = vonkar/(log(z(i)/z0water))
1144 : ustar(i,lt) = sqrt( cvar(i,lt)*uustar(i)/(1._r8 + 4.7_r8*ribn(i)) )
1145 : end if
1146 : end if
1147 : end do
1148 :
1149 : !-------------------------------------------------------------------------------------
1150 : ! compute monin-obukhov length for unstable and stable conditions/ sublayer resistance
1151 : !-------------------------------------------------------------------------------------
1152 : do lt = beglt,endlt
1153 : do i = 1,ncol
1154 : if( fr_lnduse(i,lt) ) then
1155 : hvar = (va(i)/0.74_r8) * (tha(i) - thg(i)) * (cvar(i,lt)**2)
1156 : if( unstable(i) ) then ! unstable
1157 : h = hvar*(1._r8 - (9.4_r8*ribn(i)/(1._r8 + 5.3_r8*b(i,lt))))
1158 : else
1159 : h = hvar/((1._r8+4.7_r8*ribn(i))**2)
1160 : end if
1161 : xmol(i,lt) = thg(i) * ustar(i,lt) * ustar(i,lt) / (vonkar * grav * h)
1162 : end if
1163 : end do
1164 : end do
1165 :
1166 : !-------------------------------------------------------------------------------------
1167 : ! psih
1168 : !-------------------------------------------------------------------------------------
1169 : do lt = beglt,endlt
1170 : do i = 1,ncol
1171 : if( fr_lnduse(i,lt) ) then
1172 : if( xmol(i,lt) < 0._r8 ) then
1173 : zovl = z(i)/xmol(i,lt)
1174 : zovl = max( -1._r8,zovl )
1175 : psih = exp( .598_r8 + .39_r8*log( -zovl ) - .09_r8*(log( -zovl ))**2 )
1176 : vds = 2.e-3_r8*ustar(i,lt) * (1._r8 + (300/(-xmol(i,lt)))**0.666_r8)
1177 : else
1178 : zovl = z(i)/xmol(i,lt)
1179 : zovl = min( 1._r8,zovl )
1180 : psih = -5._r8 * zovl
1181 : vds = 2.e-3_r8*ustar(i,lt)
1182 : end if
1183 : dep_ra (i,lt,lchnk) = (vonkar - psih*cvar(i,lt))/(ustar(i,lt)*vonkar*cvar(i,lt))
1184 : dep_rb (i,lt,lchnk) = (2._r8/(vonkar*ustar(i,lt))) * crb
1185 : rds(i,lt) = 1._r8/vds
1186 : end if
1187 : end do
1188 : end do
1189 :
1190 : !-------------------------------------------------------------------------------------
1191 : ! surface resistance : depends on both land type and species
1192 : ! land types are computed seperately, then resistance is computed as average of values
1193 : ! following wesely rc=(1/(rs+rm) + 1/rlu +1/(rdc+rcl) + 1/(rac+rgs))**-1
1194 : !
1195 : ! compute rsmx = 1/(rs+rm) : multiply by 3 if surface is wet
1196 : !-------------------------------------------------------------------------------------
1197 : species_loop1 : do ispec = 1,gas_pcnst
1198 : if( has_dvel(ispec) ) then
1199 : m = map_dvel(ispec)
1200 : do lt = beglt,endlt
1201 : do i = 1,ncol
1202 : if( fr_lnduse(i,lt) ) then
1203 : sndx = index_season(i,lt)
1204 : if( ispec == o3_ndx .or. ispec == o3a_ndx .or. ispec == so2_ndx ) then
1205 : rmx = 0._r8
1206 : else
1207 : rmx = 1._r8/(heff(i,m)/3000._r8 + 100._r8*foxd(m))
1208 : end if
1209 : cts(i) = 1000._r8*exp( - tc(i) - 4._r8 ) ! correction for frost
1210 : rgsx(i,lt,ispec) = cts(i) + 1._r8/((heff(i,m)/(1.e5_r8*rgss(sndx,lt))) + (foxd(m)/rgso(sndx,lt)))
1211 : !-------------------------------------------------------------------------------------
1212 : ! special case for H2 and CO;; CH4 is set ot a fraction of dv(H2)
1213 : !-------------------------------------------------------------------------------------
1214 : if( ispec == h2_ndx .or. ispec == co_ndx .or. ispec == ch4_ndx ) then
1215 : !-------------------------------------------------------------------------------------
1216 : ! no deposition on snow, ice, desert, and water
1217 : !-------------------------------------------------------------------------------------
1218 : if( lt == 1 .or. lt == 7 .or. lt == 8 .or. sndx == 4 ) then
1219 : rgsx(i,lt,ispec) = large_value
1220 : end if
1221 : end if
1222 : if( lt == 7 ) then
1223 : rclx(i,lt,ispec) = large_value
1224 : rsmx(i,lt,ispec) = large_value
1225 : rlux(i,lt,ispec) = large_value
1226 : else
1227 : rs = ri(sndx,lt)*crs(i)
1228 : if ( has_dew(i) .or. has_rain(i) ) then
1229 : dewm = 3._r8
1230 : else
1231 : dewm = 1._r8
1232 : end if
1233 : rsmx(i,lt,ispec) = (dewm*rs*drat(m) + rmx)
1234 : !-------------------------------------------------------------------------------------
1235 : ! jfl : special case for PAN
1236 : !-------------------------------------------------------------------------------------
1237 : if( ispec == pan_ndx .or. ispec == xpan_ndx ) then
1238 : dv_pan = c0_pan(lt) * (1._r8 - exp( -k_pan(lt)*(dewm*rs*drat(m))*1.e-2_r8 ))
1239 : if( dv_pan > 0._r8 .and. sndx /= 4 ) then
1240 : rsmx(i,lt,ispec) = ( 1._r8/dv_pan )
1241 : end if
1242 : end if
1243 : rclx(i,lt,ispec) = cts(i) + 1._r8/((heff(i,m)/(1.e5_r8*rcls(sndx,lt))) + (foxd(m)/rclo(sndx,lt)))
1244 : rlux(i,lt,ispec) = cts(i) + rlu(sndx,lt)/(1.e-5_r8*heff(i,m) + foxd(m))
1245 : end if
1246 : end if
1247 : end do
1248 : end do
1249 : end if
1250 : end do species_loop1
1251 :
1252 : do lt = beglt,endlt
1253 : if( lt /= 7 ) then
1254 : do i = 1,ncol
1255 : if( fr_lnduse(i,lt) ) then
1256 : sndx = index_season(i,lt)
1257 : !-------------------------------------------------------------------------------------
1258 : ! ... no effect if sfc_temp < O C
1259 : !-------------------------------------------------------------------------------------
1260 : if( sfc_temp(i) > tmelt ) then
1261 : if( has_dew(i) ) then
1262 : rlux_o3(i,lt) = 3000._r8*rlu(sndx,lt)/(1000._r8 + rlu(sndx,lt))
1263 : if( o3_ndx > 0 ) then
1264 : rlux(i,lt,o3_ndx) = rlux_o3(i,lt)
1265 : endif
1266 : if( o3a_ndx > 0 ) then
1267 : rlux(i,lt,o3a_ndx) = rlux_o3(i,lt)
1268 : endif
1269 : end if
1270 : if( has_rain(i) ) then
1271 : ! rlux(i,lt,o3_ndx) = 1./(1.e-3 + (1./(3.*rlu(sndx,lt))))
1272 : rlux_o3(i,lt) = 3000._r8*rlu(sndx,lt)/(1000._r8 + 3._r8*rlu(sndx,lt))
1273 : if( o3_ndx > 0 ) then
1274 : rlux(i,lt,o3_ndx) = rlux_o3(i,lt)
1275 : endif
1276 : if( o3a_ndx > 0 ) then
1277 : rlux(i,lt,o3a_ndx) = rlux_o3(i,lt)
1278 : endif
1279 : end if
1280 : end if
1281 :
1282 : if ( o3_ndx > 0 ) then
1283 : rclx(i,lt,o3_ndx) = cts(i) + rclo(index_season(i,lt),lt)
1284 : rlux(i,lt,o3_ndx) = cts(i) + rlux(i,lt,o3_ndx)
1285 : end if
1286 : if ( o3a_ndx > 0 ) then
1287 : rclx(i,lt,o3a_ndx) = cts(i) + rclo(index_season(i,lt),lt)
1288 : rlux(i,lt,o3a_ndx) = cts(i) + rlux(i,lt,o3a_ndx)
1289 : end if
1290 :
1291 : end if
1292 : end do
1293 : end if
1294 : end do
1295 :
1296 : species_loop2 : do ispec = 1,gas_pcnst
1297 : m = map_dvel(ispec)
1298 : if( has_dvel(ispec) ) then
1299 : if( ispec /= o3_ndx .and. ispec /= o3a_ndx .and. ispec /= so2_ndx ) then
1300 : do lt = beglt,endlt
1301 : if( lt /= 7 ) then
1302 : do i = 1,ncol
1303 : if( fr_lnduse(i,lt) ) then
1304 : !-------------------------------------------------------------------------------------
1305 : ! no effect if sfc_temp < O C
1306 : !-------------------------------------------------------------------------------------
1307 : if( sfc_temp(i) > tmelt ) then
1308 : if( has_dew(i) ) then
1309 : rlux(i,lt,ispec) = 1._r8/((1._r8/(3._r8*rlux(i,lt,ispec))) &
1310 : + 1.e-7_r8*heff(i,m) + foxd(m)/rlux_o3(i,lt))
1311 : end if
1312 : end if
1313 :
1314 : end if
1315 : end do
1316 : end if
1317 : end do
1318 : else if( ispec == so2_ndx ) then
1319 : do lt = beglt,endlt
1320 : if( lt /= 7 ) then
1321 : do i = 1,ncol
1322 : if( fr_lnduse(i,lt) ) then
1323 : !-------------------------------------------------------------------------------------
1324 : ! no effect if sfc_temp < O C
1325 : !-------------------------------------------------------------------------------------
1326 : if( sfc_temp(i) > tmelt ) then
1327 : if( qs(i) <= spec_hum(i) ) then
1328 : rlux(i,lt,ispec) = 100._r8
1329 : end if
1330 : if( has_rain(i) ) then
1331 : ! rlux(i,lt,ispec) = 1./(2.e-4 + (1./(3.*rlu(index_season(i,lt),lt))))
1332 : rlux(i,lt,ispec) = 15._r8*rlu(index_season(i,lt),lt)/(5._r8 + 3.e-3_r8*rlu(index_season(i,lt),lt))
1333 : end if
1334 : end if
1335 : rclx(i,lt,ispec) = cts(i) + rcls(index_season(i,lt),lt)
1336 : rlux(i,lt,ispec) = cts(i) + rlux(i,lt,ispec)
1337 :
1338 : end if
1339 : end do
1340 : end if
1341 : end do
1342 : do i = 1,ncol
1343 : if( fr_lnduse(i,1) .and. (has_dew(i) .or. has_rain(i)) ) then
1344 : rlux(i,1,ispec) = 50._r8
1345 : end if
1346 : end do
1347 : end if
1348 : end if
1349 : end do species_loop2
1350 :
1351 : !-------------------------------------------------------------------------------------
1352 : ! compute rc
1353 : !-------------------------------------------------------------------------------------
1354 : term(:ncol) = 1.e-2_r8 * pressure_10m(:ncol) / (r*tv(:ncol))
1355 : species_loop3 : do ispec = 1,gas_pcnst
1356 : if( has_dvel(ispec) ) then
1357 : wrk(:) = 0._r8
1358 : lt_loop: do lt = beglt,endlt
1359 : do i = 1,ncol
1360 : if (fr_lnduse(i,lt)) then
1361 : resc(i) = 1._r8/( 1._r8/rsmx(i,lt,ispec) + 1._r8/rlux(i,lt,ispec) &
1362 : + 1._r8/(rdc(i) + rclx(i,lt,ispec)) &
1363 : + 1._r8/(rac(index_season(i,lt),lt) + rgsx(i,lt,ispec)))
1364 :
1365 : resc(i) = max( 10._r8,resc(i) )
1366 :
1367 : lnd_frc(i) = lcl_frc_landuse(i,lt)
1368 : endif
1369 : enddo
1370 : !-------------------------------------------------------------------------------------
1371 : ! ... compute average deposition velocity
1372 : !-------------------------------------------------------------------------------------
1373 : select case( solsym(ispec) )
1374 : case( 'SO2' )
1375 : if( lt == 7 ) then
1376 : where( fr_lnduse(:ncol,lt) )
1377 : ! assume no surface resistance for SO2 over water`
1378 : wrk(:) = wrk(:) + lnd_frc(:)/(dep_ra(:ncol,lt,lchnk) + dep_rb(:ncol,lt,lchnk))
1379 : endwhere
1380 : else
1381 : where( fr_lnduse(:ncol,lt) )
1382 : wrk(:) = wrk(:) + lnd_frc(:)/(dep_ra(:ncol,lt,lchnk) + dep_rb(:ncol,lt,lchnk) + resc(:))
1383 : endwhere
1384 : end if
1385 :
1386 : ! JFL - increase in dry deposition of SO2 to improve bias over US/Europe
1387 : wrk(:) = wrk(:) * 2._r8
1388 :
1389 : case( 'SO4' )
1390 : where( fr_lnduse(:ncol,lt) )
1391 : wrk(:) = wrk(:) + lnd_frc(:)/(dep_ra(:ncol,lt,lchnk) + rds(:,lt))
1392 : endwhere
1393 : case( 'NH4', 'NH4NO3', 'XNH4NO3' )
1394 : where( fr_lnduse(:ncol,lt) )
1395 : wrk(:) = wrk(:) + lnd_frc(:)/(dep_ra(:ncol,lt,lchnk) + 0.5_r8*rds(:,lt))
1396 : endwhere
1397 :
1398 : !-------------------------------------------------------------------------------------
1399 : ! ... special case for Pb (for consistency with offline code)
1400 : !-------------------------------------------------------------------------------------
1401 : case( 'Pb' )
1402 : if( lt == 7 ) then
1403 : where( fr_lnduse(:ncol,lt) )
1404 : wrk(:) = wrk(:) + lnd_frc(:) * 0.05e-2_r8
1405 : endwhere
1406 : else
1407 : where( fr_lnduse(:ncol,lt) )
1408 : wrk(:ncol) = wrk(:ncol) + lnd_frc(:ncol) * 0.2e-2_r8
1409 : endwhere
1410 : end if
1411 :
1412 : !-------------------------------------------------------------------------------------
1413 : ! ... special case for carbon aerosols
1414 : !-------------------------------------------------------------------------------------
1415 : case( 'CB1', 'CB2', 'OC1', 'OC2', 'SOAM', 'SOAI', 'SOAT', 'SOAB','SOAX' )
1416 : where( fr_lnduse(:ncol,lt) )
1417 : wrk(:ncol) = wrk(:ncol) + lnd_frc(:ncol) * 0.10e-2_r8
1418 : endwhere
1419 :
1420 : !-------------------------------------------------------------------------------------
1421 : ! deposition over ocean for HCN, CH3CN
1422 : ! velocity estimated from aircraft measurements (E.Apel, INTEX-B)
1423 : !-------------------------------------------------------------------------------------
1424 : case( 'HCN','CH3CN' )
1425 : if( lt == 7 ) then ! over ocean only
1426 : where( fr_lnduse(:ncol,lt) .and. snow(:ncol) < 0.01_r8 )
1427 : wrk(:ncol) = wrk(:ncol) + lnd_frc(:ncol) * 0.2e-2_r8
1428 : endwhere
1429 : end if
1430 : case default
1431 : where( fr_lnduse(:ncol,lt) )
1432 : wrk(:ncol) = wrk(:ncol) + lnd_frc(:ncol)/(dep_ra(:ncol,lt,lchnk) + dep_rb(:ncol,lt,lchnk) + resc(:ncol))
1433 : endwhere
1434 : end select
1435 : end do lt_loop
1436 : dvel(:ncol,ispec) = wrk(:ncol) * scaling_to_cm_per_s
1437 : dflx(:ncol,ispec) = term(:ncol) * dvel(:ncol,ispec) * mmr(:ncol,plev,ispec)
1438 : end if
1439 :
1440 : end do species_loop3
1441 :
1442 : if ( beglt > 1 ) return
1443 :
1444 : !-------------------------------------------------------------------------------------
1445 : ! ... special adjustments
1446 : !-------------------------------------------------------------------------------------
1447 : if( mpan_ndx > 0 ) then
1448 : if( has_dvel(mpan_ndx) ) then
1449 : dvel(:ncol,mpan_ndx) = dvel(:ncol,mpan_ndx)/3._r8
1450 : dflx(:ncol,mpan_ndx) = term(:ncol) * dvel(:ncol,mpan_ndx) * mmr(:ncol,plev,mpan_ndx)
1451 : end if
1452 : end if
1453 : if( xmpan_ndx > 0 ) then
1454 : if( has_dvel(xmpan_ndx) ) then
1455 : dvel(:ncol,xmpan_ndx) = dvel(:ncol,xmpan_ndx)/3._r8
1456 : dflx(:ncol,xmpan_ndx) = term(:ncol) * dvel(:ncol,xmpan_ndx) * mmr(:ncol,plev,xmpan_ndx)
1457 : end if
1458 : end if
1459 :
1460 : ! HCOOH, use CH3COOH dep.vel
1461 : if( hcooh_ndx > 0) then
1462 : if( has_dvel(hcooh_ndx) ) then
1463 : dvel(:ncol,hcooh_ndx) = dvel(:ncol,ch3cooh_ndx)
1464 : dflx(:ncol,hcooh_ndx) = term(:ncol) * dvel(:ncol,hcooh_ndx) * mmr(:ncol,plev,hcooh_ndx)
1465 : end if
1466 : end if
1467 : !
1468 : ! SOG species
1469 : !
1470 : if( sogm_ndx > 0) then
1471 : if( has_dvel(sogm_ndx) ) then
1472 : dvel(:ncol,sogm_ndx) = dvel(:ncol,ch3cooh_ndx)
1473 : dflx(:ncol,sogm_ndx) = term(:ncol) * dvel(:ncol,sogm_ndx) * mmr(:ncol,plev,sogm_ndx)
1474 : end if
1475 : end if
1476 : if( sogi_ndx > 0) then
1477 : if( has_dvel(sogi_ndx) ) then
1478 : dvel(:ncol,sogi_ndx) = dvel(:ncol,ch3cooh_ndx)
1479 : dflx(:ncol,sogi_ndx) = term(:ncol) * dvel(:ncol,sogi_ndx) * mmr(:ncol,plev,sogi_ndx)
1480 : end if
1481 : end if
1482 : if( sogt_ndx > 0) then
1483 : if( has_dvel(sogt_ndx) ) then
1484 : dvel(:ncol,sogt_ndx) = dvel(:ncol,ch3cooh_ndx)
1485 : dflx(:ncol,sogt_ndx) = term(:ncol) * dvel(:ncol,sogt_ndx) * mmr(:ncol,plev,sogt_ndx)
1486 : end if
1487 : end if
1488 : if( sogb_ndx > 0) then
1489 : if( has_dvel(sogb_ndx) ) then
1490 : dvel(:ncol,sogb_ndx) = dvel(:ncol,ch3cooh_ndx)
1491 : dflx(:ncol,sogb_ndx) = term(:ncol) * dvel(:ncol,sogb_ndx) * mmr(:ncol,plev,sogb_ndx)
1492 : end if
1493 : end if
1494 : if( sogx_ndx > 0) then
1495 : if( has_dvel(sogx_ndx) ) then
1496 : dvel(:ncol,sogx_ndx) = dvel(:ncol,ch3cooh_ndx)
1497 : dflx(:ncol,sogx_ndx) = term(:ncol) * dvel(:ncol,sogx_ndx) * mmr(:ncol,plev,sogx_ndx)
1498 : end if
1499 : end if
1500 :
1501 0 : end subroutine drydep_xactive
1502 :
1503 : !-------------------------------------------------------------------------------------
1504 : !-------------------------------------------------------------------------------------
1505 0 : function has_drydep( name )
1506 :
1507 : character(len=*), intent(in) :: name
1508 :
1509 : logical :: has_drydep
1510 : integer :: i
1511 :
1512 0 : has_drydep = .false.
1513 :
1514 0 : do i=1,nddvels
1515 0 : if ( trim(name) == trim(drydep_list(i)) ) then
1516 0 : has_drydep = .true.
1517 0 : exit
1518 : endif
1519 : enddo
1520 :
1521 0 : endfunction has_drydep
1522 :
1523 0 : end module mo_drydep
|
