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