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