Original WRF subgrid support version from John Michalakes without fire
[wrffire.git] / wrfv2_fire / chem / module_dep_simple.F
blob07e54ed4111d2607714dfcc15508893a53fc9783
1     MODULE module_dep_simple
3 ! many of these parameters will depend on the RADM mechanism!
4 ! if you change it, lets talk about it and get it done!!!
7       INTEGER, PARAMETER :: dep_seasons = 5 ,   nlu = 25,                   &
8         nseason = 1, nseasons = 2
10 ! following currently hardwired to USGS
12       INTEGER, PARAMETER :: isice_temp=24,iswater_temp=16
13       character, parameter :: mminlu='USGS'
15       INTEGER :: ixxxlu(nlu)
16       REAL :: kpart(nlu),                                                   &
17         rac(nlu,dep_seasons), rclo(nlu,dep_seasons), rcls(nlu,dep_seasons), &
18         rgso(nlu,dep_seasons), rgss(nlu,dep_seasons),                       &
19         ri(nlu,dep_seasons), rlu(nlu,dep_seasons)
21 ! NO MORE THAN 1000 SPECIES FOR DEPOSITION
23       REAL, DIMENSION (1:1000) :: dratio,hstar,hstar4,f0,dhr,scpr23
24 ! .. Default Accessibility ..
25     PUBLIC
26             logical, allocatable :: is_aerosol(:) ! true if field is aerosol (any phase)
27 ! ..
28     CONTAINS
29        subroutine wesely_driver(id,ktau,dtstep,                           &
30                config_flags,                                              &
31                gmt,julday,t_phy,moist,p8w,t8w,                            &
32                p_phy,chem,rho_phy,dz8w,ddvel,aer_res,                     &
33                ivgtyp,tsk,gsw,vegfra,pbl,rmol,ust,znt,xlat,xlong,z,z_at_w,&
34                numgas,                                                    &
35                ids,ide, jds,jde, kds,kde,                                 &
36                ims,ime, jms,jme, kms,kme,                                 &
37                its,ite, jts,jte, kts,kte                                  )
38 !----------------------------------------------------------------------
39   USE module_model_constants 
40   USE module_configure
41   USE module_state_description                       
42   USE module_data_sorgam
43 ! USE module_data_radm2                               
44 ! USE module_data_racm                               
45           
46    INTEGER,      INTENT(IN   ) :: id,julday,                              &
47                                   numgas,                                 &
48                                   ids,ide, jds,jde, kds,kde,              &
49                                   ims,ime, jms,jme, kms,kme,              &
50                                   its,ite, jts,jte, kts,kte     
51    INTEGER,      INTENT(IN   ) ::                                         &
52                                   ktau            
53       REAL,      INTENT(IN   ) ::                                         &
54                              dtstep,gmt
56 ! advected moisture variables
58    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),               &
59          INTENT(IN ) ::                                   moist  
61 ! advected chemical species
63    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                &
64          INTENT(INOUT ) ::                                   chem
66 ! deposition velocities
68    REAL, DIMENSION( its:ite, jts:jte, num_chem ),                         &
69          INTENT(INOUT ) ::   ddvel                     
71 ! input from met model
74    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,              &
75           INTENT(IN   ) ::                                                &
76                                                       t_phy,              &
77                                                       p_phy,              &
78                                                       dz8w,               &
79                                                       z    ,              &
80                                               t8w,p8w,z_at_w ,            &
81                                                     rho_phy
82    INTEGER,DIMENSION( ims:ime , jms:jme )                  ,              &
83           INTENT(IN   ) ::                                                &
84                                                      ivgtyp
85    REAL,  DIMENSION( ims:ime , jms:jme )                   ,              &
86           INTENT(INOUT   ) ::                                                &
87                                                      tsk,                 &
88                                                      gsw,                 &
89                                                   vegfra,                 &
90                                                      pbl,                 &
91                                                      rmol,                &
92                                                      ust,                 &
93                                                      xlat,                &
94                                                      xlong,               &
95                                                      znt
97 !--- deposition and emissions stuff
98 ! ..
99 ! .. Local Scalars ..
100       REAL ::  clwchem,dvfog,dvpart,rad,rhchem,ta,ustar,vegfrac,z1,zntt
101       INTEGER :: iland, iprt, iseason, jce, jcs, n, nr, ipr,jpr,nvr
102       LOGICAL :: highnh3, rainflag, vegflag, wetflag
103 ! ..
104 ! .. Local Arrays ..  
105       REAL :: p(kts:kte-1), srfres(numgas),ddvel0d(numgas)
107 ! necessary for aerosols (module dependent)         
108 !  
109       real :: aer_res(its:ite,jts:jte),rcx(numgas)
110                                                      
111    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
112           
113                                                      
114 ! ..                                                 
115 ! .. Intrinsic Functions ..                       
116       INTRINSIC max, min                             
117 ! ..                                                 
118       CALL wrf_debug(15,'in dry_dep_wesely')
119       iseason = 1
120       if(julday.lt.90.or.julday.gt.270)then
121         iseason=2
122         CALL wrf_debug(15,'setting iseason to 2')
123       endif
124       do 100 j=jts,jte 
125       do 100 i=its,ite 
126       iprt=0
127       iland = ivgtyp(i,j)
128       ta = tsk(i,j)
129       rad = gsw(i,j)
130       vegfrac = vegfra(i,j)
131       pa = .01*p_phy(i,kts,j)
132       clwchem = moist(i,kts,j,p_qc)
133       ustar = ust(i,j)
134       zntt = znt(i,j)
136       z1 = z_at_w(i,kts+1,j)-z_at_w(i,kts,j)
138 !     Set logical default values
139       rainflag = .FALSE.
140       wetflag = .FALSE.
141       highnh3 = .FALSE.
142       if(p_qr.gt.1)then
143          if(moist(i,kts,j,p_qr).gt.0.)rainflag = .true.
144       endif
145       rhchem = MIN( 100.,100. * moist(i,kts,j,p_qv) / &
146                (3.80*exp(17.27*(t_phy(i,kts,j)-273.)/(t_phy(i,kts,j)-36.))/pa))
147       rhchem = MAX(5.,RHCHEM)
148       if (rhchem >= 95.) wetflag = .true.
150       if(chem(i,kts,j,p_nh3).gt.2.*chem(i,kts,j,p_so2))highnh3 = .true.
152 !--- deposition
153       
154 !     if(snowc(i,j).gt.0.)iseason=4
155       CALL rc(rcx,ta,rad,rhchem,iland,iseason,numgas,      &
156         wetflag,rainflag,highnh3,iprt,p_o3,p_so2,p_nh3)
157       srfres=0.
158       DO n = 1, numgas-2
159         srfres(n) = rcx(n)
160       END DO
161       CALL deppart(rmol(i,j),ustar,rhchem,clwchem,iland,dvpart,dvfog)
162       ddvel0d=0.
163       aer_res(i,j)=0.
164       CALL landusevg(ddvel0d,ustar,rmol(i,j),zntt,z1,dvpart,iland,        &
165            numgas,srfres,aer_res(i,j),p_sulf)
167 !wig: CMBZ doe snot have Ho and HO2 last so need to copy all species
168 !      ddvel(i,j,1:numgas-2)=ddvel0d(1:numgas-2)
169       ddvel(i,j,1:numgas)=ddvel0d(1:numgas)
171 100   continue
173 ! For the additional CBMZ species, assign similar RADM counter parts for
174 ! now. Short lived species get a zero velocity since dry dep should be
175 ! unimportant.  **ALSO**, treat p_sulf as h2so4 vapor, not aerosol sulfate
177       if ( (config_flags%chem_opt == CBMZ          ) .or.   &
178            (config_flags%chem_opt == CBMZ_BB       ) .or.   &
179            (config_flags%chem_opt == CBMZ_MOSAIC_4BIN_AQ) .or.   &
180            (config_flags%chem_opt == CBMZ_MOSAIC_8BIN_AQ) .or.   &
181            (config_flags%chem_opt == CBMZ_MOSAIC_4BIN) .or.   &
182            (config_flags%chem_opt == CBMZ_MOSAIC_8BIN) ) then
183          do j=jts,jte
184             do i=its,ite
185                ddvel(i,j,p_sulf)        = ddvel(i,j,p_hno3)
186                ddvel(i,j,p_hcl)         = ddvel(i,j,p_hno3)
187                ddvel(i,j,p_hcl)         = ddvel(i,j,p_hno3)
188                ddvel(i,j,p_ch3o2)       = 0
189                ddvel(i,j,p_ethp)        = 0
190                ddvel(i,j,p_ch3oh)       = ddvel(i,j,p_hcho)
191                ddvel(i,j,p_c2h5oh)      = ddvel(i,j,p_hcho)
192                ddvel(i,j,p_par)         = ddvel(i,j,p_hc5)
193                ddvel(i,j,p_to2)         = 0
194                ddvel(i,j,p_cro)         = 0
195                ddvel(i,j,p_open)        = ddvel(i,j,p_xyl)
196                ddvel(i,j,p_op3)         = ddvel(i,j,p_op2)
197                ddvel(i,j,p_c2o3)        = 0
198                ddvel(i,j,p_ro2)         = 0
199                ddvel(i,j,p_ano2)        = 0
200                ddvel(i,j,p_nap)         = 0
201                ddvel(i,j,p_xo2)         = 0
202                ddvel(i,j,p_xpar)        = 0
203                ddvel(i,j,p_isoprd)      = 0
204                ddvel(i,j,p_isopp)       = 0
205                ddvel(i,j,p_isopn)       = 0
206                ddvel(i,j,p_isopo2)      = 0
207                if( config_flags%chem_opt == CBMZ ) then
208                   ddvel(i,j,p_dms)         = 0
209                   ddvel(i,j,p_msa)         = ddvel(i,j,p_hno3)
210                   ddvel(i,j,p_dmso)        = 0
211                   ddvel(i,j,p_dmso2)       = 0
212                   ddvel(i,j,p_ch3so2h)     = 0
213                   ddvel(i,j,p_ch3sch2oo)   = 0
214                   ddvel(i,j,p_ch3so2)      = 0
215                   ddvel(i,j,p_ch3so3)      = 0
216                   ddvel(i,j,p_ch3so2oo)    = 0
217                   ddvel(i,j,p_ch3so2ch2oo) = 0
218                   ddvel(i,j,p_mtf)         = 0
219                end if
220             end do
221          end do
222       end if
224 END SUBROUTINE wesely_driver
226 ! **********************************************************************
227 ! **************************  SUBROUTINE RC  ***************************
228 ! **********************************************************************
229       SUBROUTINE rc(rcx,t,rad,rh,iland,iseason,numgas,             &
230           wetflag,rainflag,highnh3,iprt,p_o3,p_so2,p_nh3)
231 !     THIS SUBROUTINE CALCULATES SURFACE RESISTENCES ACCORDING
232 !     TO THE MODEL OF
233 !     M. L. WESELY,
234 !     ATMOSPHERIC ENVIRONMENT 23 (1989), 1293-1304
235 !     WITH SOME ADDITIONS ACCORDING TO
236 !     J. W. ERISMAN, A. VAN PUL, AND P. WYERS,
237 !     ATMOSPHERIC ENVIRONMENT 28 (1994), 2595-2607
238 !     WRITTEN BY  WINFRIED SEIDL, APRIL 1997
239 !     MODYFIED BY WINFRIED SEIDL, MARCH 2000
240 !                    FOR MM5 VERSION 3
241 !----------------------------------------------------------------------
242 ! .. Scalar Arguments ..
243         REAL :: rad, rh, t
244         INTEGER :: iland, iseason, numgas
245         LOGICAL :: highnh3, rainflag, wetflag
246         real :: rcx(numgas)
247 ! ..
248         INTEGER :: iprt,p_o3,p_so2,p_nh3
249 ! ..
250 ! .. Local Scalars ..
251         REAL :: rclx, rdc, resice, rgsx, rluo1, rluo2, rlux, rmx, rs, rsmx, &
252           tc, rdtheta, z
253         INTEGER :: n
254 ! ..
255 ! .. Local Arrays ..
256         REAL :: hstary(numgas)
257 ! ..
258 ! .. Intrinsic Functions ..
259         INTRINSIC exp
260 ! ..
261         DO n = 1, numgas
262           rcx(n) = 1.
263         END DO
265         tc = t - 273.15
266         rdtheta = 0.
268         z = 200./(rad+0.1)
270 !!!  HARDWIRE VALUES FOR TESTING
271 !       z=0.4727409
272 !       tc=22.76083
273 !       t=tc+273.15
274 !       rad = 412.8426
275 !       rainflag=.false.
276 !       wetflag=.false.
278         IF ((tc<=0.) .OR. (tc>=40.)) THEN
279           rs = 9999.
280         ELSE
281           rs = ri(iland,iseason)*(1+z*z)*(400./(tc*(40.-tc)))
282         END IF
283         rdc = 100*(1.+1000./(rad+10))/(1+1000.*rdtheta)
284         rluo1 = 1./(1./3000.+1./3./rlu(iland,iseason))
285         rluo2 = 1./(1./1000.+1./3./rlu(iland,iseason))
286         resice = 1000.*exp(-tc-4.)
288         DO n = 1, numgas
289           IF (hstar(n)==0.) GO TO 10
290           hstary(n) = hstar(n)*exp(dhr(n)*(1./t-1./298.))
291           rmx = 1./(hstary(n)/3000.+100.*f0(n))
292           rsmx = rs*dratio(n) + rmx
293           rclx = 1./(hstary(n)/1.E+5/rcls(iland,iseason)+f0(n)/rclo(iland, &
294             iseason)) + resice
295           rgsx = 1./(hstary(n)/1.E+5/rgss(iland,iseason)+f0(n)/rgso(iland, &
296             iseason)) + resice
297           rlux = rlu(iland,iseason)/(1.E-5*hstary(n)+f0(n)) + resice
298           IF (wetflag) THEN
299             rlux = 1./(1./3./rlu(iland,iseason)+1.E-7*hstary(n)+f0(n)/rluo1)
300           END IF
301           IF (rainflag) THEN
302             rlux = 1./(1./3./rlu(iland,iseason)+1.E-7*hstary(n)+f0(n)/rluo2)
303           END IF
304           rcx(n) = 1./(1./rsmx+1./rlux+1./(rdc+rclx)+1./(rac(iland, &
305             iseason)+rgsx))
306           IF (rcx(n)<1.) rcx(n) = 1.
307 10      END DO
309 !     SPECIAL TREATMENT FOR OZONE
310         hstary(p_o3) = hstar(p_o3)*exp(dhr(p_o3)*(1./t-1./298.))
311         rmx = 1./(hstary(p_o3)/3000.+100.*f0(p_o3))
312         rsmx = rs*dratio(p_o3) + rmx
313         rlux = rlu(iland,iseason)/(1.E-5*hstary(p_o3)+f0(p_o3)) + resice
314         rclx = rclo(iland,iseason) + resice
315         rgsx = rgso(iland,iseason) + resice
316         IF (wetflag) rlux = rluo1
317         IF (rainflag) rlux = rluo2
318         rcx(p_o3) = 1./(1./rsmx+1./rlux+1./(rdc+rclx)+1./(rac(iland, &
319           iseason)+rgsx))
320         IF (rcx(p_o3)<1.) rcx(p_o3) = 1.
322 !     SPECIAL TREATMENT FOR SO2 (Wesely)
323 !       HSTARY(P_SO2)=HSTAR(P_SO2)*EXP(DHR(P_SO2)*(1./T-1./298.))
324 !       RMX=1./(HSTARY(P_SO2)/3000.+100.*F0(P_SO2))
325 !       RSMX=RS*DRATIO(P_SO2)+RMX
326 !       RLUX=RLU(ILAND,ISEASON)/(1.E-5*HSTARY(P_SO2)+F0(P_SO2))
327 !    &       +RESICE
328 !       RCLX=RCLS(ILAND,ISEASON)+RESICE
329 !       RGSX=RGSS(ILAND,ISEASON)+RESICE
330 !       IF ((wetflag).OR.(RAINFLAG)) THEN
331 !         IF (ILAND.EQ.1) THEN
332 !           RLUX=50.
333 !         ELSE
334 !           RLUX=100.
335 !         END IF
336 !       END IF
337 !       RCX(P_SO2)=1./(1./RSMX+1./RLUX+1./(RDC+RCLX)
338 !    &                +1./(RAC(ILAND,ISEASON)+RGSX))
339 !       IF (RCX(P_SO2).LT.1.) RCX(P_SO2)=1.
341 !     SO2 according to Erisman et al. 1994
342 !       R_STOM
343         rsmx = rs*dratio(p_so2)
344 !       R_EXT
345         IF (tc>(-1.)) THEN
346           IF (rh<81.3) THEN
347             rlux = 25000.*exp(-0.0693*rh)
348           ELSE
349             rlux = 0.58E12*exp(-0.278*rh)
350           END IF
351         END IF
352         IF (((wetflag) .OR. (rainflag)) .AND. (tc>(-1.))) THEN
353           rlux = 1.
354         END IF
355         IF ((tc>=(-5.)) .AND. (tc<=(-1.))) THEN
356           rlux = 200.
357         END IF
358         IF (tc<(-5.)) THEN
359           rlux = 500.
360         END IF
361 !       INSTEAD OF R_INC R_CL and R_DC of Wesely are used
362         rclx = rcls(iland,iseason)
363 !       DRY SURFACE
364         rgsx = 1000.
365 !       WET SURFACE
366         IF ((wetflag) .OR. (rainflag)) THEN
367           IF (highnh3) THEN
368             rgsx = 0.
369           ELSE
370             rgsx = 500.
371           END IF
372         END IF
373 !       WATER
374         IF (iland==iswater_temp) THEN
375           rgsx = 0.
376         END IF
377 !       SNOW
378         IF ((iseason==4) .OR. (iland==isice_temp)) THEN
379           IF (tc>2.) THEN
380             rgsx = 0.
381           END IF
382           IF ((tc>=(-1.)) .AND. (tc<=2.)) THEN
383             rgsx = 70.*(2.-tc)
384           END IF
385           IF (tc<(-1.)) THEN
386             rgsx = 500.
387           END IF
388         END IF
389 !       TOTAL SURFACE RESISTENCE
390         IF ((iseason/=4) .AND. (ixxxlu(iland)/=1) .AND. (iland/=iswater_temp) .AND. &
391             (iland/=isice_temp)) THEN
392           rcx(p_so2) = 1./(1./rsmx+1./rlux+1./(rclx+rdc+rgsx))
393         ELSE
394           rcx(p_so2) = rgsx
395         END IF
396         IF (rcx(p_so2)<1.) rcx(p_so2) = 1.
397 !     NH3 according to Erisman et al. 1994
398 !       R_STOM
399         rsmx = rs*dratio(p_nh3)
400 !       GRASSLAND (PASTURE DURING GRAZING)
401         IF (ixxxlu(iland)==3) THEN
402           IF (iseason==1) THEN
403 !           SUMMER
404             rcx(p_nh3) = 1000.
405           END IF
406           IF ((iseason==2) .OR. (iseason==3) .OR. (iseason==5)) THEN
407 !           WINTER, NO SNOW
408             IF (tc>-1.) THEN
409               IF (rad/=0.) THEN
410                 rcx(p_nh3) = 50.
411               ELSE
412                 rcx(p_nh3) = 100.
413               END IF
414               IF ((wetflag) .OR. (rainflag)) THEN
415                 rcx(p_nh3) = 20.
416               END IF
417             END IF
418             IF ((tc>=(-5.)) .AND. (tc<=-1.)) THEN
419               rcx(p_nh3) = 200.
420             END IF
421             IF (tc<(-5.)) THEN
422               rcx(p_nh3) = 500.
423             END IF
424           END IF
425         END IF
426 !       AGRICULTURAL LAND (CROPS AND UNGRAZED PASTURE)
427         IF (ixxxlu(iland)==2) THEN
428           IF (iseason==1) THEN
429 !           SUMMER
430             IF (rad/=0.) THEN
431               rcx(p_nh3) = rsmx
432             ELSE
433               rcx(p_nh3) = 200.
434             END IF
435             IF ((wetflag) .OR. (rainflag)) THEN
436               rcx(p_nh3) = 50.
437             END IF
438           END IF
439           IF ((iseason==2) .OR. (iseason==3) .OR. (iseason==5)) THEN
440 !           WINTER, NO SNOW
441             IF (tc>-1.) THEN
442               IF (rad/=0.) THEN
443                 rcx(p_nh3) = rsmx
444               ELSE
445                 rcx(p_nh3) = 300.
446               END IF
447               IF ((wetflag) .OR. (rainflag)) THEN
448                 rcx(p_nh3) = 100.
449               END IF
450             END IF
451             IF ((tc>=(-5.)) .AND. (tc<=-1.)) THEN
452               rcx(p_nh3) = 200.
453             END IF
454             IF (tc<(-5.)) THEN
455               rcx(p_nh3) = 500.
456             END IF
457           END IF
458         END IF
459 !       SEMI-NATURAL ECOSYSTEMS AND FORESTS
460         IF ((ixxxlu(iland)==4) .OR. (ixxxlu(iland)==5) .OR. (ixxxlu( &
461             iland)==6)) THEN
462           IF (rad/=0.) THEN
463             rcx(p_nh3) = 500.
464           ELSE
465             rcx(p_nh3) = 1000.
466           END IF
467           IF ((wetflag) .OR. (rainflag)) THEN
468             IF (highnh3) THEN
469               rcx(p_nh3) = 100.
470             ELSE
471               rcx(p_nh3) = 0.
472             END IF
473           END IF
474           IF ((iseason==2) .OR. (iseason==3) .OR. (iseason==5)) THEN
475 !           WINTER, NO SNOW
476             IF ((tc>=(-5.)) .AND. (tc<=-1.)) THEN
477               rcx(p_nh3) = 200.
478             END IF
479             IF (tc<(-5.)) THEN
480               rcx(p_nh3) = 500.
481             END IF
482           END IF
483         END IF
484 !       WATER
485         IF (iland==iswater_temp) THEN
486           rcx(p_nh3) = 0.
487         END IF
488 !       URBAN AND DESERT (SOIL SURFACES)
489         IF (ixxxlu(iland)==1) THEN
490           IF ( .NOT. wetflag) THEN
491             rcx(p_nh3) = 50.
492           ELSE
493             rcx(p_nh3) = 0.
494           END IF
495         END IF
496 !       SNOW COVERED SURFACES OR PERMANENT ICE
497         IF ((iseason==4) .OR. (iland==isice_temp)) THEN
498           IF (tc>2.) THEN
499             rcx(p_nh3) = 0.
500           END IF
501           IF ((tc>=(-1.)) .AND. (tc<=2.)) THEN
502             rcx(p_nh3) = 70.*(2.-tc)
503           END IF
504           IF (tc<(-1.)) THEN
505             rcx(p_nh3) = 500.
506           END IF
507         END IF
508         IF (rcx(p_nh3)<1.) rcx(p_nh3) = 1.
510       END SUBROUTINE rc
511 ! **********************************************************************
512 ! ************************ SUBROUTINE DEPPART **************************
513 ! **********************************************************************
514       SUBROUTINE deppart(rmol,ustar,rh,clw,iland,dvpart,dvfog)
515 !     THIS SUBROUTINE CALCULATES SURFACE DEPOSITION VELOCITIES
516 !     FOR FINE AEROSOL PARTICLES ACCORDING TO THE MODEL OF
517 !     J. W. ERISMAN, A. VAN PUL, AND P. WYERS,
518 !     ATMOSPHERIC ENVIRONMENT 28 (1994), 2595-2607
519 !     WRITTEN BY WINFRIED SEIDL, APRIL 1997
520 !     MODIFIED BY WINFRIED SEIDL, MARCH 2000
521 !            FOR MM5 VERSION 3
522 ! ----------------------------------------------------------------------
523 ! .. Scalar Arguments ..
524         REAL :: clw, dvfog, dvpart, rh, rmol, ustar
525         INTEGER :: iland
526 ! ..
527 ! .. Intrinsic Functions ..
528         INTRINSIC exp
529 ! ..
530         dvpart = ustar/kpart(iland)
531         IF (rmol<0.) THEN
532 !         INSTABLE LAYERING CORRECTION
533           dvpart = dvpart*(1.+(-300.*rmol)**0.66667)
534         END IF
535         IF (rh>80.) THEN
536 !         HIGH RELATIVE HUMIDITY CORRECTION
537 !         ACCORDING TO J. W. ERISMAN ET AL.
538 !         ATMOSPHERIC ENVIRONMENT 31 (1997), 321-332
539           dvpart = dvpart*(1.+0.37*exp((rh-80.)/20.))
540         END IF
542 !       SEDIMENTATION VELOCITY OF FOG WATER ACCORDING TO
543 !       R. FORKEL, W. SEIDL, R. DLUGI AND E. DEIGELE
544 !       J. GEOPHYS. RES. 95D (1990), 18501-18515
545         dvfog = 0.06*clw
546         IF (ixxxlu(iland)==5) THEN
547 !         TURBULENT DEPOSITION OF FOG WATER IN CONIFEROUS FOREST ACCORDI
548 !         A. T. VERMEULEN ET AL.
549 !         ATMOSPHERIC ENVIRONMENT 31 (1997), 375-386
550           dvfog = dvfog + 0.195*ustar*ustar
551         END IF
553       END SUBROUTINE deppart
554       SUBROUTINE landusevg(vgs,ustar,rmol,z0,zz,dvparx,iland,numgas, &
555           srfres,aer_res,p_sulf)
556 !     This subroutine calculates the species specific deposition velocit
557 !     as a function of the local meteorology and land use.  The depositi
558 !     Velocity is also landuse specific.
559 !     Reference: Hsieh, C.M., Wesely, M.L. and Walcek, C.J. (1986)
560 !                A Dry Deposition Module for Regional Acid Deposition
561 !                EPA report under agreement DW89930060-01
562 !     Revised version by Darrell Winner (January 1991)
563 !        Environmental Engineering Science 138-78
564 !           California Institute of Technology
565 !              Pasadena, CA  91125
566 !     Modified by Winfried Seidl (August 1997)
567 !       Fraunhofer-Institut fuer Atmosphaerische Umweltforschung
568 !                    Garmisch-Partenkirchen, D-82467
569 !          for use of Wesely and Erisman surface resistances
570 !     Inputs:
571 !        Ustar  : The grid average friction velocity (m/s)
572 !        Rmol   : Reciprocal of the Monin-Obukhov length (1/m)
573 !        Z0     : Surface roughness height for the grid square (m)
574 !        SrfRes : Array of landuse/atmospheric/species resistances (s/m)
575 !        Slist  : Array of chemical species codes
576 !        Dvparx : Array of surface deposition velocity of fine aerosol p
577 !     Outputs:
578 !        Vgs    : Array of species and landuse specific deposition
579 !                 velocities (m/s)
580 !        Vg     : Cell-average deposition velocity by species (m/s)
581 !     Variables used:
582 !        SCPR23  : (Schmidt #/Prandtl #)**(2/3) Diffusion correction fac
583 !        Zr      : Reference Height (m)
584 !        Iatmo   : Parameter specifying the stabilty class (Function of
585 !        Z0      : Surface roughness height (m)
586 !        karman  : Von Karman constant (from module_model_constants)
587         USE module_model_constants, only: karman
588 !     Local Variables
589 ! .. Scalar Arguments ..
590         REAL :: dvparx, rmol, ustar, z0, zz
591         real :: aer_res, polint
592         INTEGER :: iland, numgas, p_sulf
593 ! ..
594 ! .. Array Arguments ..
595         REAL :: srfres(numgas), vgs(numgas)
596 ! ..
597 ! .. Local Scalars ..
598         REAL :: vgp, vgpart, zr
599         INTEGER :: jspec
600 ! ..
601 ! .. Local Arrays ..
602         REAL :: vgspec(numgas)
603 ! ..
604 !   Set the reference height (10.0 m)
605         zr = 10.0
607 !   CALCULATE THE DEPOSITION VELOCITY without any surface
608 !   resistance term, i.e. 1 / (ra + rb)
609         CALL depvel(numgas,rmol,zr,z0,ustar,vgspec,vgpart,aer_res)
611 !   Calculate the deposition velocity for each species
612 !   and grid cell by looping through all the possibile combinations
613 !   of the two
615         vgp = 1.0/((1.0/vgpart)+(1.0/dvparx))
617 !   Loop through the various species
619         DO jspec = 1, numgas
621 !   Add in the surface resistance term, rc (SrfRes)
623           vgs(jspec) = 1.0/(1.0/vgspec(jspec)+srfres(jspec))
624         END DO
625         vgs(p_sulf) = vgp
627         CALL cellvg(vgs,ustar,zz,zr,rmol,numgas)
629         RETURN
630       END SUBROUTINE landusevg
632       SUBROUTINE cellvg(vgtemp,ustar,dz,zr,rmol,nspec)
633 !     THIS PROGRAM HAS BEEN DESIGNED TO CALCULATE THE CELL AVERAGE
634 !     DEPOSITION VELOCITY GIVEN THE VALUE OF VG AT SOME REFERENCE
635 !     HEIGHT ZR WHICH IS MUCH SMALLER THAN THE CELL HEIGHT DZ.
636 !       PROGRAM WRITTEN BY GREGORY J.MCRAE (NOVEMBER 1977)
637 !         Modified by Darrell A. Winner    (February 1991)
638 !.....PROGRAM VARIABLES...
639 !     VgTemp   - DEPOSITION VELOCITY AT THE REFERENCE HEIGHT
640 !     USTAR    - FRICTION VELOCITY
641 !     RMOL     - RECIPROCAL OF THE MONIN-OBUKHOV LENGTH
642 !     ZR       - REFERENCE HEIGHT
643 !     DZ       - CELL HEIGHT
644 !     CELLVG   - CELL AVERAGE DEPOSITION VELOCITY
645 !     VK       - VON KARMAN CONSTANT
646         USE module_model_constants, only: karman
647 !     Local Variables
648 ! .. Scalar Arguments ..
649         REAL :: dz, rmol, ustar, zr
650         INTEGER :: nspec
651 ! ..
652 ! .. Array Arguments ..
653         REAL :: vgtemp(nspec)
654 ! ..
655 ! .. Local Scalars ..
656         REAL :: a, fac, pdz, pzr, vk
657         INTEGER :: nss
658 ! ..
659 ! .. Intrinsic Functions ..
660         INTRINSIC alog, sqrt
661 ! ..
662 !     Set the von Karman constant
663         vk = karman
665 !     DETERMINE THE STABILITY BASED ON THE CONDITIONS
666 !             1/L < 0 UNSTABLE
667 !             1/L = 0 NEUTRAL
668 !             1/L > 0 STABLE
671         DO nss = 1, nspec
672           IF (rmol<0) THEN
673             pdz = sqrt(1.0-9.0*dz*rmol)
674             pzr = sqrt(1.0-9.0*zr*rmol)
675             fac = ((pdz-1.0)/(pzr-1.0))*((pzr+1.0)/(pdz+1.0))
676             a = 0.74*dz*alog(fac) + (0.164/rmol)*(pdz-pzr)
677           ELSE IF (rmol==0) THEN
678             a = 0.74*(dz*alog(dz/zr)-dz+zr)
679           ELSE
680             a = 0.74*(dz*alog(dz/zr)-dz+zr) + (2.35*rmol)*(dz-zr)**2
681           END IF
683 !     CALCULATE THE DEPOSITION VELOCITIY
685           vgtemp(nss) = vgtemp(nss)/(1.0+vgtemp(nss)*a/(vk*ustar*(dz-zr)))
686         END DO
688         RETURN
689       END SUBROUTINE cellvg
690       SUBROUTINE depvel(numgas,rmol,zr,z0,ustar,depv,vgpart,aer_res)
691 !     THIS FUNCTION HAS BEEN DESIGNED TO EVALUATE AN UPPER LIMIT
692 !     FOR THE POLLUTANT DEPOSITION VELOCITY AS A FUNCTION OF THE
693 !     SURFACE ROUGHNESS AND METEOROLOGICAL CONDITIONS.
694 !     PROGRAM WRITTEN BY GREGORY J.MCRAE (NOVEMBER 1977)
695 !         Modified by Darrell A. Winner  (Feb. 1991)
696 !                  by Winfried Seidl     (Aug. 1997)
697 !.....PROGRAM VARIABLES...
698 !     RMOL     - RECIPROCAL OF THE MONIN-OBUKHOV LENGTH
699 !     ZR       - REFERENCE HEIGHT
700 !     Z0       - SURFACE ROUGHNESS HEIGHT
701 !     SCPR23   - (Schmidt #/Prandtl #)**(2/3) Diffusion correction fact
702 !     UBAR     - ABSOLUTE VALUE OF SURFACE WIND SPEED
703 !     DEPVEL   - POLLUTANT DEPOSITION VELOCITY
704 !     Vk       - VON KARMAN CONSTANT
705 !     USTAR    - FRICTION VELOCITY U*
706 !     POLINT   - POLLUTANT INTEGRAL
707 !     AER_RES  - AERODYNAMIC RESISTANCE
708 !.....REFERENCES...
709 !     MCRAE, G.J. ET AL. (1983) MATHEMATICAL MODELING OF PHOTOCHEMICAL
710 !       AIR POLLUTION, ENVIRONMENTAL QUALITY LABORATORY REPORT 18,
711 !       CALIFORNIA INSTITUTE OF TECHNOLOGY, PASADENA, CALIFORNIA.
712 !.....RESTRICTIONS...
713 !     1. THE MODEL EDDY DIFFUSIVITIES ARE BASED ON MONIN-OBUKHOV
714 !        SIMILARITY THEORY AND SO ARE ONLY APPLICABLE IN THE
715 !        SURFACE LAYER, A HEIGHT OF O(30M).
716 !     2. ALL INPUT UNITS MUST BE CONSISTENT
717 !     3. THE PHI FUNCTIONS USED TO CALCULATE THE FRICTION
718 !        VELOCITY U* AND THE POLLUTANT INTEGRALS ARE BASED
719 !        ON THE WORK OF BUSINGER ET AL.(1971).
720 !     4. THE MOMENTUM AND POLLUTANT DIFFUSIVITIES ARE NOT
721 !        THE SAME FOR THE CASES L<0 AND L>0.
722         USE module_model_constants, only: karman
723 !     Local Variables
724 ! .. Scalar Arguments ..
725         REAL :: rmol, ustar, vgpart, z0, zr, aer_res
726         INTEGER :: numgas
727 ! ..
728 ! .. Array Arguments ..
729         REAL :: depv(numgas)
730 ! ..
731 ! .. Local Scalars ..
732         REAL :: ao, ar, polint, vk
733         INTEGER :: l
734 ! ..
735 ! .. Intrinsic Functions ..
736         INTRINSIC alog
737 ! ..
738 !     Set the von Karman constant
739         vk = karman
741 !     Calculate the diffusion correction factor
742 !     SCPR23 is calculated as (Sc/Pr)**(2/3) using Sc= 1.15 and Pr= 1.0
743 !     SCPR23 = 1.10
745 !     DETERMINE THE STABILITY BASED ON THE CONDITIONS
746 !             1/L < 0 UNSTABLE
747 !             1/L = 0 NEUTRAL
748 !             1/L > 0 STABLE
750         if(abs(rmol) < 1.E-6 ) rmol = 0.
752         IF (rmol<0) THEN
753           ar = ((1.0-9.0*zr*rmol)**(0.25)+0.001)**2
754           ao = ((1.0-9.0*z0*rmol)**(0.25)+0.001)**2
755           polint = 0.74*(alog((ar-1.0)/(ar+1.0))-alog((ao-1.0)/(ao+1.0)))
756         ELSE IF (rmol==0.) THEN
757           polint = 0.74*alog(zr/z0)
758         ELSE
759           polint = 0.74*alog(zr/z0) + 4.7*rmol*(zr-z0)
760         END IF
762 !     CALCULATE THE Maximum DEPOSITION VELOCITY
764         DO l = 1, numgas
765           depv(l) = ustar*vk/(2.0*scpr23(l)+polint)
766         END DO
767         vgpart = ustar*vk/polint
768         aer_res = polint/(karman*max(ustar,1.0e-4))
770         RETURN
771       END SUBROUTINE depvel
773       SUBROUTINE dep_init(id,config_flags,numgas)
774   USE module_model_constants
775   USE module_configure
776   USE module_state_description                       
777    TYPE (grid_config_rec_type) , INTENT (in) ::     config_flags
778 ! ..
779 ! .. Scalar Arguments ..
780         integer, intent(in) :: id, numgas
782 ! ..
783 ! .. Local Scalars ..
784         REAL :: sc
785         INTEGER :: iland, iseason, l
786         integer :: iprt
787 ! ..
788 ! .. Local Arrays ..
789         REAL :: dat1(nlu,dep_seasons), dat2(nlu,dep_seasons),         &
790                 dat3(nlu,dep_seasons), dat4(nlu,dep_seasons),         &
791                 dat5(nlu,dep_seasons), dat6(nlu,dep_seasons),         &
792                 dat7(nlu,dep_seasons), dvj(numgas)
793 ! ..
794 ! .. Make sure that the model is being run with a soil model. Otherwise,
795 !    iland will be zero in deppart, which will try to pull non-exisant
796 !    array locations.
797         call nl_get_sf_surface_physics(id,l)
798         if( l == 0 ) &
799              call wrf_error_fatal("ERROR: Cannot use dry deposition without using a soil model.")
801 ! ..
802 ! .. Data Statements ..
803 !     RI for stomatal resistance
804 !      data ((ri(ILAND,ISEASON),ILAND=1,nlu),ISEASON=1,dep_seasons)/0.10E+11, &
805         DATA ((dat1(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
806           0.60E+02, 0.60E+02, 0.60E+02, 0.60E+02, 0.70E+02, 0.12E+03, &
807           0.12E+03, 0.12E+03, 0.12E+03, 0.70E+02, 0.13E+03, 0.70E+02, &
808           0.13E+03, 0.10E+03, 0.10E+11, 0.80E+02, 0.10E+03, 0.10E+11, &
809           0.80E+02, 0.10E+03, 0.10E+03, 0.10E+11, 0.10E+11, 0.10E+11, &
810           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
811           0.10E+11, 0.10E+11, 0.10E+11, 0.12E+03, 0.10E+11, 0.10E+11, &
812           0.70E+02, 0.25E+03, 0.50E+03, 0.10E+11, 0.10E+11, 0.50E+03, &
813           0.10E+11, 0.10E+11, 0.50E+03, 0.50E+03, 0.10E+11, 0.10E+11, &
814           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
815           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.12E+03, 0.10E+11, &
816           0.10E+11, 0.70E+02, 0.25E+03, 0.50E+03, 0.10E+11, 0.10E+11, &
817           0.50E+03, 0.10E+11, 0.10E+11, 0.50E+03, 0.50E+03, 0.10E+11, &
818           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
819           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
820           0.10E+11, 0.10E+11, 0.70E+02, 0.40E+03, 0.80E+03, 0.10E+11, &
821           0.10E+11, 0.80E+03, 0.10E+11, 0.10E+11, 0.80E+03, 0.80E+03, &
822           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.12E+03, 0.12E+03, &
823           0.12E+03, 0.12E+03, 0.14E+03, 0.24E+03, 0.24E+03, 0.24E+03, &
824           0.12E+03, 0.14E+03, 0.25E+03, 0.70E+02, 0.25E+03, 0.19E+03, &
825           0.10E+11, 0.16E+03, 0.19E+03, 0.10E+11, 0.16E+03, 0.19E+03, &
826           0.19E+03, 0.10E+11, 0.10E+11, 0.10E+11/
827 ! ..
828         IF (nlu/=25) THEN
829           call wrf_debug(0, 'number of land use classifications not correct ')
830           CALL wrf_error_fatal ( "LAND USE CLASSIFICATIONS NOT 25")
831         END IF
832         IF (dep_seasons/=5) THEN
833           call wrf_debug(0, 'number of dep_seasons not correct ')
834           CALL wrf_error_fatal ( "DEP_SEASONS NOT 5")
835         END IF
837 !     SURFACE RESISTANCE DATA FOR DEPOSITION MODEL OF
838 !     M. L. WESELY, ATMOSPHERIC ENVIRONMENT 23 (1989) 1293-1304
840 !     Seasonal categories:
841 !     1: midsummer with lush vegetation
842 !     2: autumn with unharvested cropland
843 !     3: late autumn with frost, no snow
844 !     4: winter, snow on ground and subfreezing
845 !     5: transitional spring with partially green short annuals
847 !     Land use types:
848 !     USGS type                                Wesely type
849 !      1: Urban and built-up land              1
850 !      2: Dryland cropland and pasture         2
851 !      3: Irrigated cropland and pasture       2
852 !      4: Mix. dry/irrg. cropland and pasture  2
853 !      5: Cropland/grassland mosaic            2
854 !      6: Cropland/woodland mosaic             4
855 !      7: Grassland                            3
856 !      8: Shrubland                            3
857 !      9: Mixed shrubland/grassland            3
858 !     10: Savanna                              3, always summer
859 !     11: Deciduous broadleaf forest           4
860 !     12: Deciduous needleleaf forest          5, autumn and winter modi
861 !     13: Evergreen broadleaf forest           4, always summer
862 !     14: Evergreen needleleaf forest          5
863 !     15: Mixed Forest                         6
864 !     16: Water Bodies                         7
865 !     17: Herbaceous wetland                   9
866 !     18: Wooded wetland                       6
867 !     19: Barren or sparsely vegetated         8
868 !     20: Herbaceous Tundra                    9
869 !     21: Wooded Tundra                        6
870 !     22: Mixed Tundra                         6
871 !     23: Bare Ground Tundra                   8
872 !     24: Snow or Ice                          -, always winter
873 !     25: No data                              8
876 !     Order of data:
877 !      |
878 !      |   seasonal category
879 !     \|/
880 !     ---> landuse type
881 !     1       2       3       4       5       6       7       8       9
882 !     RLU for outer surfaces in the upper canopy
883         DO iseason = 1, dep_seasons
884           DO iland = 1, nlu
885             ri(iland,iseason) = dat1(iland,iseason)
886           END DO
887         END DO
888 !      data ((rlu(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+11, &
889         DATA ((dat2(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
890           0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
891           0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
892           0.20E+04, 0.20E+04, 0.10E+11, 0.25E+04, 0.20E+04, 0.10E+11, &
893           0.25E+04, 0.20E+04, 0.20E+04, 0.10E+11, 0.10E+11, 0.10E+11, &
894           0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, &
895           0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, 0.90E+04, &
896           0.20E+04, 0.40E+04, 0.80E+04, 0.10E+11, 0.90E+04, 0.80E+04, &
897           0.10E+11, 0.90E+04, 0.80E+04, 0.80E+04, 0.10E+11, 0.10E+11, &
898           0.10E+11, 0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, &
899           0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, &
900           0.90E+04, 0.20E+04, 0.40E+04, 0.80E+04, 0.10E+11, 0.90E+04, &
901           0.80E+04, 0.10E+11, 0.90E+04, 0.80E+04, 0.80E+04, 0.10E+11, &
902           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
903           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
904           0.10E+11, 0.10E+11, 0.20E+04, 0.60E+04, 0.90E+04, 0.10E+11, &
905           0.90E+04, 0.90E+04, 0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, &
906           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.40E+04, 0.40E+04, &
907           0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, &
908           0.20E+04, 0.40E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.30E+04, &
909           0.10E+11, 0.40E+04, 0.30E+04, 0.10E+11, 0.40E+04, 0.30E+04, &
910           0.30E+04, 0.10E+11, 0.10E+11, 0.10E+11/
911         DO iseason = 1, dep_seasons
912           DO iland = 1, nlu
913             rlu(iland,iseason) = dat2(iland,iseason)
914           END DO
915         END DO
916 !     RAC for transfer that depends on canopy height and density
917 !      data ((rac(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+03, &
918         DATA ((dat3(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+03, &
919           0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+04, 0.10E+03, &
920           0.10E+03, 0.10E+03, 0.10E+03, 0.20E+04, 0.20E+04, 0.20E+04, &
921           0.20E+04, 0.20E+04, 0.00E+00, 0.30E+03, 0.20E+04, 0.00E+00, &
922           0.30E+03, 0.20E+04, 0.20E+04, 0.00E+00, 0.00E+00, 0.00E+00, &
923           0.10E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+04, &
924           0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.15E+04, 0.20E+04, &
925           0.20E+04, 0.20E+04, 0.17E+04, 0.00E+00, 0.20E+03, 0.17E+04, &
926           0.00E+00, 0.20E+03, 0.17E+04, 0.17E+04, 0.00E+00, 0.00E+00, &
927           0.00E+00, 0.10E+03, 0.10E+02, 0.10E+02, 0.10E+02, 0.10E+02, &
928           0.10E+04, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+04, &
929           0.20E+04, 0.20E+04, 0.20E+04, 0.15E+04, 0.00E+00, 0.10E+03, &
930           0.15E+04, 0.00E+00, 0.10E+03, 0.15E+04, 0.15E+04, 0.00E+00, &
931           0.00E+00, 0.00E+00, 0.10E+03, 0.10E+02, 0.10E+02, 0.10E+02, &
932           0.10E+02, 0.10E+04, 0.10E+02, 0.10E+02, 0.10E+02, 0.10E+02, &
933           0.10E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.15E+04, 0.00E+00, &
934           0.50E+02, 0.15E+04, 0.00E+00, 0.50E+02, 0.15E+04, 0.15E+04, &
935           0.00E+00, 0.00E+00, 0.00E+00, 0.10E+03, 0.50E+02, 0.50E+02, &
936           0.50E+02, 0.50E+02, 0.12E+04, 0.80E+02, 0.80E+02, 0.80E+02, &
937           0.10E+03, 0.12E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.15E+04, &
938           0.00E+00, 0.20E+03, 0.15E+04, 0.00E+00, 0.20E+03, 0.15E+04, &
939           0.15E+04, 0.00E+00, 0.00E+00, 0.00E+00/
940         DO iseason = 1, dep_seasons
941           DO iland = 1, nlu
942             rac(iland,iseason) = dat3(iland,iseason)
943           END DO
944         END DO
945 !     RGSS for ground surface  SO2
946 !      data ((rgss(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.40E+03, &
947         DATA ((dat4(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.40E+03, &
948           0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.50E+03, 0.35E+03, &
949           0.35E+03, 0.35E+03, 0.35E+03, 0.50E+03, 0.50E+03, 0.50E+03, &
950           0.50E+03, 0.10E+03, 0.10E+01, 0.10E+01, 0.10E+03, 0.10E+04, &
951           0.10E+01, 0.10E+03, 0.10E+03, 0.10E+04, 0.10E+03, 0.10E+04, &
952           0.40E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.50E+03, &
953           0.35E+03, 0.35E+03, 0.35E+03, 0.35E+03, 0.50E+03, 0.50E+03, &
954           0.50E+03, 0.50E+03, 0.10E+03, 0.10E+01, 0.10E+01, 0.10E+03, &
955           0.10E+04, 0.10E+01, 0.10E+03, 0.10E+03, 0.10E+04, 0.10E+03, &
956           0.10E+04, 0.40E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, &
957           0.50E+03, 0.35E+03, 0.35E+03, 0.35E+03, 0.35E+03, 0.50E+03, &
958           0.50E+03, 0.50E+03, 0.50E+03, 0.20E+03, 0.10E+01, 0.10E+01, &
959           0.20E+03, 0.10E+04, 0.10E+01, 0.20E+03, 0.20E+03, 0.10E+04, &
960           0.10E+03, 0.10E+04, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, &
961           0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, &
962           0.10E+03, 0.10E+03, 0.50E+03, 0.10E+03, 0.10E+03, 0.10E+01, &
963           0.10E+03, 0.10E+03, 0.10E+04, 0.10E+03, 0.10E+03, 0.10E+03, &
964           0.10E+04, 0.10E+03, 0.10E+04, 0.50E+03, 0.15E+03, 0.15E+03, &
965           0.15E+03, 0.15E+03, 0.50E+03, 0.35E+03, 0.35E+03, 0.35E+03, &
966           0.35E+03, 0.50E+03, 0.50E+03, 0.50E+03, 0.50E+03, 0.20E+03, &
967           0.10E+01, 0.10E+01, 0.20E+03, 0.10E+04, 0.10E+01, 0.20E+03, &
968           0.20E+03, 0.10E+04, 0.10E+03, 0.10E+04/
969         DO iseason = 1, dep_seasons
970           DO iland = 1, nlu
971             rgss(iland,iseason) = dat4(iland,iseason)
972           END DO
973         END DO
974 !     RGSO for ground surface  O3
975 !      data ((rgso(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.30E+03, &
976         DATA ((dat5(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.30E+03, &
977           0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.20E+03, 0.20E+03, &
978           0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
979           0.20E+03, 0.30E+03, 0.20E+04, 0.10E+04, 0.30E+03, 0.40E+03, &
980           0.10E+04, 0.30E+03, 0.30E+03, 0.40E+03, 0.35E+04, 0.40E+03, &
981           0.30E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.20E+03, &
982           0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
983           0.20E+03, 0.20E+03, 0.30E+03, 0.20E+04, 0.80E+03, 0.30E+03, &
984           0.40E+03, 0.80E+03, 0.30E+03, 0.30E+03, 0.40E+03, 0.35E+04, &
985           0.40E+03, 0.30E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, &
986           0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
987           0.20E+03, 0.20E+03, 0.20E+03, 0.30E+03, 0.20E+04, 0.10E+04, &
988           0.30E+03, 0.40E+03, 0.10E+04, 0.30E+03, 0.30E+03, 0.40E+03, &
989           0.35E+04, 0.40E+03, 0.60E+03, 0.35E+04, 0.35E+04, 0.35E+04, &
990           0.35E+04, 0.35E+04, 0.35E+04, 0.35E+04, 0.35E+04, 0.35E+04, &
991           0.35E+04, 0.35E+04, 0.20E+03, 0.35E+04, 0.35E+04, 0.20E+04, &
992           0.35E+04, 0.35E+04, 0.40E+03, 0.35E+04, 0.35E+04, 0.35E+04, &
993           0.40E+03, 0.35E+04, 0.40E+03, 0.30E+03, 0.15E+03, 0.15E+03, &
994           0.15E+03, 0.15E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
995           0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.30E+03, &
996           0.20E+04, 0.10E+04, 0.30E+03, 0.40E+03, 0.10E+04, 0.30E+03, &
997           0.30E+03, 0.40E+03, 0.35E+04, 0.40E+03/
998         DO iseason = 1, dep_seasons
999           DO iland = 1, nlu
1000             rgso(iland,iseason) = dat5(iland,iseason)
1001           END DO
1002         END DO
1003 !     RCLS for exposed surfaces in the lower canopy  SO2
1004 !      data ((rcls(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+11, &
1005         DATA ((dat6(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
1006           0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
1007           0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
1008           0.20E+04, 0.20E+04, 0.10E+11, 0.25E+04, 0.20E+04, 0.10E+11, &
1009           0.25E+04, 0.20E+04, 0.20E+04, 0.10E+11, 0.10E+11, 0.10E+11, &
1010           0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, &
1011           0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, 0.90E+04, &
1012           0.20E+04, 0.20E+04, 0.40E+04, 0.10E+11, 0.90E+04, 0.40E+04, &
1013           0.10E+11, 0.90E+04, 0.40E+04, 0.40E+04, 0.10E+11, 0.10E+11, &
1014           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
1015           0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, &
1016           0.90E+04, 0.20E+04, 0.30E+04, 0.60E+04, 0.10E+11, 0.90E+04, &
1017           0.60E+04, 0.10E+11, 0.90E+04, 0.60E+04, 0.60E+04, 0.10E+11, &
1018           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
1019           0.10E+11, 0.90E+04, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
1020           0.90E+04, 0.90E+04, 0.20E+04, 0.20E+03, 0.40E+03, 0.10E+11, &
1021           0.90E+04, 0.40E+03, 0.10E+11, 0.90E+04, 0.40E+03, 0.40E+03, &
1022           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.40E+04, 0.40E+04, &
1023           0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, &
1024           0.20E+04, 0.40E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.30E+04, &
1025           0.10E+11, 0.40E+04, 0.30E+04, 0.10E+11, 0.40E+04, 0.30E+04, &
1026           0.30E+04, 0.10E+11, 0.10E+11, 0.10E+11/
1027         DO iseason = 1, dep_seasons
1028           DO iland = 1, nlu
1029             rcls(iland,iseason) = dat6(iland,iseason)
1030           END DO
1031         END DO
1032 !     RCLO for exposed surfaces in the lower canopy  O3
1033 !      data ((rclo(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+11, &
1034         DATA ((dat7(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
1035           0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
1036           0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
1037           0.10E+04, 0.10E+04, 0.10E+11, 0.10E+04, 0.10E+04, 0.10E+11, &
1038           0.10E+04, 0.10E+04, 0.10E+04, 0.10E+11, 0.10E+11, 0.10E+11, &
1039           0.10E+11, 0.40E+03, 0.40E+03, 0.40E+03, 0.40E+03, 0.40E+03, &
1040           0.40E+03, 0.40E+03, 0.40E+03, 0.10E+04, 0.40E+03, 0.40E+03, &
1041           0.10E+04, 0.10E+04, 0.60E+03, 0.10E+11, 0.40E+03, 0.60E+03, &
1042           0.10E+11, 0.40E+03, 0.60E+03, 0.60E+03, 0.10E+11, 0.10E+11, &
1043           0.10E+11, 0.10E+11, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
1044           0.40E+03, 0.40E+03, 0.40E+03, 0.40E+03, 0.10E+04, 0.40E+03, &
1045           0.40E+03, 0.10E+04, 0.10E+04, 0.60E+03, 0.10E+11, 0.80E+03, &
1046           0.60E+03, 0.10E+11, 0.80E+03, 0.60E+03, 0.60E+03, 0.10E+11, &
1047           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+04, 0.10E+04, 0.10E+04, &
1048           0.10E+04, 0.40E+03, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
1049           0.40E+03, 0.40E+03, 0.10E+04, 0.15E+04, 0.60E+03, 0.10E+11, &
1050           0.80E+03, 0.60E+03, 0.10E+11, 0.80E+03, 0.60E+03, 0.60E+03, &
1051           0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+04, 0.10E+04, &
1052           0.10E+04, 0.10E+04, 0.50E+03, 0.50E+03, 0.50E+03, 0.50E+03, &
1053           0.10E+04, 0.50E+03, 0.15E+04, 0.10E+04, 0.15E+04, 0.70E+03, &
1054           0.10E+11, 0.60E+03, 0.70E+03, 0.10E+11, 0.60E+03, 0.70E+03, &
1055           0.70E+03, 0.10E+11, 0.10E+11, 0.10E+11/
1056         DO iseason = 1, dep_seasons
1057           DO iland = 1, nlu
1058             rclo(iland,iseason) = dat7(iland,iseason)
1059           END DO
1060         END DO
1061         DO l = 1, numgas
1062           hstar(l) = 0.
1063           hstar4(l) = 0.
1064           dhr(l) = 0.
1065           f0(l) = 0.
1066           dvj(l) = 99.
1067         END DO
1069 !     HENRY''S LAW COEFFICIENTS
1070 !     Effective Henry''s law coefficient at pH 7
1071 !     [KH298]=mole/(l atm)
1073         hstar(p_no2) = 6.40E-3
1074         hstar(p_no) = 1.90E-3
1075         hstar(p_pan) = 2.97E+0
1076         hstar(p_o3) = 1.13E-2
1077         hstar(p_hcho) = 2.97E+3
1078         hstar(p_aco3) = 1.14E+1
1079         hstar(p_tpan) = 2.97E+0
1080         hstar(p_hono) = 3.47E+5
1081         hstar(p_no3) = 1.50E+1
1082         hstar(p_hno4) = 2.00E+13
1083         hstar(p_h2o2) = 7.45E+4
1084         hstar(p_co) = 8.20E-3
1085         hstar(p_ald) = 1.14E+1
1086         hstar(p_op1) = 2.21E+2
1087         hstar(p_op2) = 1.68E+6
1088         hstar(p_paa) = 4.73E+2
1089         hstar(p_ket) = 3.30E+1
1090         hstar(p_gly) = 1.40E+6
1091         hstar(p_mgly) = 3.71E+3
1092         hstar(p_dcb) = 1.40E+6
1093         hstar(p_onit) = 1.13E+0
1094         hstar(p_so2) = 2.53E+5
1095         hstar(p_eth) = 2.00E-3
1096         hstar(p_hc3) = 1.42E-3
1097         hstar(p_hc5) = 1.13E-3
1098         hstar(p_hc8) = 1.42E-3
1099         hstar(p_olt) = 4.76E-3
1100         hstar(p_oli) = 1.35E-3
1101         hstar(p_tol) = 1.51E-1
1102         hstar(p_csl) = 4.00E+5
1103         hstar(p_xyl) = 1.45E-1
1104         hstar(p_iso) = 4.76E-3
1105         hstar(p_hno3) = 2.69E+13
1106         hstar(p_ora1) = 9.85E+6
1107         hstar(p_ora2) = 9.63E+5
1108         hstar(p_nh3) = 1.04E+4
1109         hstar(p_n2o5) = 1.00E+10
1110         if(p_ol2.gt.1)hstar(p_ol2) = 4.67E-3
1112 !  FOLLOWING FOR RACM
1114         if(p_ete.gt.1)then
1115            HSTAR(p_ETE )=4.67E-3
1116            HSTAR(p_API )=4.76E-3
1117            HSTAR(p_LIM )=4.76E-3
1118            HSTAR(p_DIEN)=4.76E-3
1119            HSTAR(p_CH4 )=1.50E-3
1120            HSTAR(p_CO2 )=1.86E-1
1121            HSTAR(p_MACR)=1.14E+1
1122            HSTAR(p_UDD )=1.40E+6
1123            HSTAR(p_HKET)=7.80E+3
1124            DHR(p_ETE )=    0.
1125            DHR(p_API )=    0.
1126            DHR(p_LIM )=    0.
1127            DHR(p_DIEN)=    0.
1128            DHR(p_CH4 )=    0.
1129            DHR(p_CO2 )= 1636.
1130            DHR(p_MACR)= 6266.
1131            DHR(p_UDD )=    0.
1132            DHR(p_HKET)=    0.
1133            F0(p_ETE )=0.
1134            F0(p_API )=0.
1135            F0(p_LIM )=0.
1136            F0(p_DIEN)=0.
1137            F0(p_CH4 )=0.
1138            F0(p_CO2 )=0.
1139            F0(p_MACR)=0.
1140            F0(p_UDD )=0.
1141            F0(p_HKET)=0.
1142            DVJ(p_ETE )=0.189
1143            DVJ(p_API )=0.086
1144            DVJ(p_LIM )=0.086
1145            DVJ(p_DIEN)=0.136
1146            DVJ(p_CH4 )=0.250
1147            DVJ(p_CO2 )=0.151
1148            DVJ(p_MACR)=0.120
1149            DVJ(p_UDD )=0.092
1150            DVJ(p_HKET)=0.116
1151         endif
1152 !     -DH/R (for temperature correction)
1153 !     [-DH/R]=K
1155         dhr(p_no2) = 2500.
1156         dhr(p_no) = 1480.
1157         dhr(p_pan) = 5760.
1158         dhr(p_o3) = 2300.
1159         dhr(p_hcho) = 7190.
1160         dhr(p_aco3) = 6266.
1161         dhr(p_tpan) = 5760.
1162         dhr(p_hono) = 3775.
1163         dhr(p_no3) = 0.
1164         dhr(p_hno4) = 0.
1165         dhr(p_h2o2) = 6615.
1166         dhr(p_co) = 0.
1167         dhr(p_ald) = 6266.
1168         dhr(p_op1) = 5607.
1169         dhr(p_op2) = 10240.
1170         dhr(p_paa) = 6170.
1171         dhr(p_ket) = 5773.
1172         dhr(p_gly) = 0.
1173         dhr(p_mgly) = 7541.
1174         dhr(p_dcb) = 0.
1175         dhr(p_onit) = 5487.
1176         dhr(p_so2) = 5816.
1177         dhr(p_eth) = 0.
1178         dhr(p_hc3) = 0.
1179         dhr(p_hc5) = 0.
1180         dhr(p_hc8) = 0.
1181         dhr(p_olt) = 0.
1182         dhr(p_oli) = 0.
1183         dhr(p_tol) = 0.
1184         dhr(p_csl) = 0.
1185         dhr(p_xyl) = 0.
1186         dhr(p_iso) = 0.
1187         dhr(p_hno3) = 8684.
1188         dhr(p_ora1) = 5716.
1189         dhr(p_ora2) = 8374.
1190         dhr(p_nh3) = 3660.
1191         dhr(p_n2o5) = 0.
1192         if(p_ol2.gt.1)dhr(p_ol2) = 0.
1193 !     REACTIVITY FACTORS
1194 !     [f0]=1
1196         f0(p_no2) = 0.1
1197         f0(p_no) = 0.
1198         f0(p_pan) = 0.1
1199         f0(p_o3) = 1.
1200         f0(p_hcho) = 0.
1201         f0(p_aco3) = 1.
1202         f0(p_tpan) = 0.1
1203         f0(p_hono) = 0.1
1204         f0(p_no3) = 1.
1205         f0(p_hno4) = 0.1
1206         f0(p_h2o2) = 1.
1207         f0(p_co) = 0.
1208         f0(p_ald) = 0.
1209         f0(p_op1) = 0.1
1210         f0(p_op2) = 0.1
1211         f0(p_paa) = 0.1
1212         f0(p_ket) = 0.
1213         f0(p_gly) = 0.
1214         f0(p_mgly) = 0.
1215         f0(p_dcb) = 0.
1216         f0(p_onit) = 0.
1217         f0(p_so2) = 0.
1218         f0(p_eth) = 0.
1219         f0(p_hc3) = 0.
1220         f0(p_hc5) = 0.
1221         f0(p_hc8) = 0.
1222         f0(p_olt) = 0.
1223         f0(p_oli) = 0.
1224         f0(p_tol) = 0.
1225         f0(p_csl) = 0.
1226         f0(p_xyl) = 0.
1227         f0(p_iso) = 0.
1228         f0(p_hno3) = 0.
1229         f0(p_ora1) = 0.
1230         f0(p_ora2) = 0.
1231         f0(p_nh3) = 0.
1232         f0(p_n2o5) = 1.
1233         if(p_ol2.gt.1)f0(p_ol2) = 0.
1234 !     DIFFUSION COEFFICIENTS
1235 !     [DV]=cm2/s (assumed: 1/SQRT(molar mass) when not known)
1237         dvj(p_no2) = 0.147
1238         dvj(p_no) = 0.183
1239         dvj(p_pan) = 0.091
1240         dvj(p_o3) = 0.175
1241         dvj(p_hcho) = 0.183
1242         dvj(p_aco3) = 0.115
1243         dvj(p_tpan) = 0.082
1244         dvj(p_hono) = 0.153
1245         dvj(p_no3) = 0.127
1246         dvj(p_hno4) = 0.113
1247         dvj(p_h2o2) = 0.171
1248         dvj(p_co) = 0.189
1249         dvj(p_ald) = 0.151
1250         dvj(p_op1) = 0.144
1251         dvj(p_op2) = 0.127
1252         dvj(p_paa) = 0.115
1253         dvj(p_ket) = 0.118
1254         dvj(p_gly) = 0.131
1255         dvj(p_mgly) = 0.118
1256         dvj(p_dcb) = 0.107
1257         dvj(p_onit) = 0.092
1258         dvj(p_so2) = 0.126
1259         dvj(p_eth) = 0.183
1260         dvj(p_hc3) = 0.151
1261         dvj(p_hc5) = 0.118
1262         dvj(p_hc8) = 0.094
1263         dvj(p_olt) = 0.154
1264         dvj(p_oli) = 0.121
1265         dvj(p_tol) = 0.104
1266         dvj(p_csl) = 0.096
1267         dvj(p_xyl) = 0.097
1268         dvj(p_iso) = 0.121
1269         dvj(p_hno3) = 0.126
1270         dvj(p_ora1) = 0.153
1271         dvj(p_ora2) = 0.124
1272         dvj(p_nh3) = 0.227
1273         dvj(p_n2o5) = 0.110
1274         dvj(p_ho) = 0.243
1275         dvj(p_ho2) = 0.174
1276         if(p_ol2.gt.1)dvj(p_ol2) = 0.189
1277         DO l = 1, numgas
1278           hstar4(l) = hstar(l) ! preliminary              
1279 ! Correction of diff. coeff
1280           dvj(l) = dvj(l)*(293.15/298.15)**1.75
1281           sc = 0.15/dvj(l) ! Schmidt Number at 20°C   
1282           dratio(l) = 0.242/dvj(l) !                                            ! of water vapor and gas at
1283 ! Ratio of diffusion coeffi
1284           scpr23(l) = (sc/0.72)**(2./3.) ! (Schmidt # / Prandtl #)**
1285         END DO
1289 !     DATA FOR AEROSOL PARTICLE DEPOSITION FOR THE MODEL OF
1290 !     J. W. ERISMAN, A. VAN PUL AND P. WYERS
1291 !     ATMOSPHERIC ENVIRONMENT 28 (1994), 2595-2607
1293 !     vd = (u* / k) * CORRECTION FACTORS
1295 !     CONSTANT K FOR LANDUSE TYPES:
1296 ! urban and built-up land                  
1297         kpart(1) = 500.
1298 ! dryland cropland and pasture             
1299         kpart(2) = 500.
1300 ! irrigated cropland and pasture           
1301         kpart(3) = 500.
1302 ! mixed dryland/irrigated cropland and past
1303         kpart(4) = 500.
1304 ! cropland/grassland mosaic                
1305         kpart(5) = 500.
1306 ! cropland/woodland mosaic                 
1307         kpart(6) = 100.
1308 ! grassland                                
1309         kpart(7) = 500.
1310 ! shrubland                                
1311         kpart(8) = 500.
1312 ! mixed shrubland/grassland                
1313         kpart(9) = 500.
1314 ! savanna                                  
1315         kpart(10) = 500.
1316 ! deciduous broadleaf forest               
1317         kpart(11) = 100.
1318 ! deciduous needleleaf forest              
1319         kpart(12) = 100.
1320 ! evergreen broadleaf forest               
1321         kpart(13) = 100.
1322 ! evergreen needleleaf forest              
1323         kpart(14) = 100.
1324 ! mixed forest                             
1325         kpart(15) = 100.
1326 ! water bodies                             
1327         kpart(16) = 500.
1328 ! herbaceous wetland                       
1329         kpart(17) = 500.
1330 ! wooded wetland                           
1331         kpart(18) = 500.
1332 ! barren or sparsely vegetated             
1333         kpart(19) = 500.
1334 ! herbaceous tundra                        
1335         kpart(20) = 500.
1336 ! wooded tundra                            
1337         kpart(21) = 100.
1338 ! mixed tundra                             
1339         kpart(22) = 500.
1340 ! bare ground tundra                       
1341         kpart(23) = 500.
1342 ! snow or ice                              
1343         kpart(24) = 500.
1344 !     Comments:
1345         kpart(25) = 500.
1346 !     Erisman et al. (1994) give
1347 !     k = 500 for low vegetation and k = 100 for forests.
1349 !     For desert k = 500 is taken according to measurements
1350 !     on bare soil by
1351 !     J. Fontan, A. Lopez, E. Lamaud and A. Druilhet (1997)
1352 !     Vertical Flux Measurements of the Submicronic Aerosol Particles
1353 !     and Parametrisation of the Dry Deposition Velocity
1354 !     in: Biosphere-Atmosphere Exchange of Pollutants
1355 !     and Trace Substances
1356 !     Editor: S. Slanina. Springer-Verlag Berlin, Heidelberg, 1997
1357 !     pp. 381-390
1359 !     For coniferous forest the Erisman value of  k = 100 is taken.
1360 !     Measurements of Erisman et al. (1997) in a coniferous forest
1361 !     in the Netherlands, lead to values of k between 20 and 38
1362 !     (Atmospheric Environment 31 (1997), 321-332).
1363 !     However, these high values of vd may be reached during
1364 !     instable cases. The eddy correlation measurements
1365 !     of Gallagher et al. (1997) made during the same experiment
1366 !     show for stable cases (L>0) values of k between 200 and 250
1367 !     at minimum (Atmospheric Environment 31 (1997), 359-373).
1368 !     Fontan et al. (1997) found k = 250 in a forest
1369 !     of maritime pine in southwestern France.
1371 !     For gras, model calculations of Davidson et al. support
1372 !     the value of 500.
1373 !     C. I. Davidson, J. M. Miller and M. A. Pleskov
1374 !     The Influence of Surface Structure on Predicted Particles
1375 !     Dry Deposition to Natural Gras Canopies
1376 !     Water, Air, and Soil Pollution 18 (1982) 25-43
1378 !     Snow covered surface: The experiment of Ibrahim et al. (1983)
1379 !     gives k = 436 for 0.7 um diameter particles.
1380 !     The deposition velocity of Milford and Davidson (1987)
1381 !     gives k = 154 for continental sulfate aerosol.
1382 !     M. Ibrahim, L. A. Barrie and F. Fanaki
1383 !     Atmospheric Environment 17 (1983), 781-788
1385 !     J. B. Milford and C. I. Davidson
1386 !     The Sizes of Particulate Sulfate and Nitrate in the Atmosphere
1387 !     - A Review
1388 !     JAPCA 37 (1987), 125-134
1389 ! no data                                  
1390 !       WRITE (0,*) ' return from rcread '
1391 !     *********************************************************
1393 !     Simplified landuse scheme for deposition and biogenic emission
1394 !     subroutines
1395 !     (ISWATER and ISICE are already defined elsewhere,
1396 !     therefore water and ice are not considered here)
1398 !     1 urban or bare soil
1399 !     2 agricultural
1400 !     3 grassland
1401 !     4 deciduous forest
1402 !     5 coniferous and mixed forest
1403 !     6 other natural landuse categories
1406         IF (mminlu=='OLD ') THEN
1407           ixxxlu(1) = 1
1408           ixxxlu(2) = 2
1409           ixxxlu(3) = 3
1410           ixxxlu(4) = 4
1411           ixxxlu(5) = 5
1412           ixxxlu(6) = 5
1413           ixxxlu(7) = 0
1414           ixxxlu(8) = 6
1415           ixxxlu(9) = 1
1416           ixxxlu(10) = 6
1417           ixxxlu(11) = 0
1418           ixxxlu(12) = 4
1419           ixxxlu(13) = 6
1420         END IF
1421         IF (mminlu=='USGS') THEN
1422           ixxxlu(1) = 1
1423           ixxxlu(2) = 2
1424           ixxxlu(3) = 2
1425           ixxxlu(4) = 2
1426           ixxxlu(5) = 2
1427           ixxxlu(6) = 4
1428           ixxxlu(7) = 3
1429           ixxxlu(8) = 6
1430           ixxxlu(9) = 3
1431           ixxxlu(10) = 6
1432           ixxxlu(11) = 4
1433           ixxxlu(12) = 5
1434           ixxxlu(13) = 4
1435           ixxxlu(14) = 5
1436           ixxxlu(15) = 5
1437           ixxxlu(16) = 0
1438           ixxxlu(17) = 6
1439           ixxxlu(18) = 4
1440           ixxxlu(19) = 1
1441           ixxxlu(20) = 6
1442           ixxxlu(21) = 4
1443           ixxxlu(22) = 6
1444           ixxxlu(23) = 1
1445           ixxxlu(24) = 0
1446           ixxxlu(25) = 1
1447         END IF
1448         IF (mminlu=='SiB ') THEN
1449           ixxxlu(1) = 4
1450           ixxxlu(2) = 4
1451           ixxxlu(3) = 4
1452           ixxxlu(4) = 5
1453           ixxxlu(5) = 5
1454           ixxxlu(6) = 6
1455           ixxxlu(7) = 3
1456           ixxxlu(8) = 6
1457           ixxxlu(9) = 6
1458           ixxxlu(10) = 6
1459           ixxxlu(11) = 1
1460           ixxxlu(12) = 2
1461           ixxxlu(13) = 6
1462           ixxxlu(14) = 1
1463           ixxxlu(15) = 0
1464           ixxxlu(16) = 0
1465           ixxxlu(17) = 1
1466         END IF
1467         RETURN
1468       END SUBROUTINE dep_init
1469     END MODULE module_dep_simple