| blob | 19c00787836cc5371c2515c4296e6e5ce5992f0f |
1 MODULE module_gocart_drydep
3 CONTAINS
5 subroutine gocart_drydep_driver(dtstep, &
6 config_flags,numgas, &
7 t_phy,moist,p8w,t8w,rmol,aer_res, &
8 p_phy,chem,rho_phy,dz8w,ddvel,xland,hfx, &
9 ivgtyp,tsk,vegfra,pbl,ust,znt,xlat,xlong, &
10 dustdrydep_1,dustdrydep_2,dustdrydep_3, &
11 dustdrydep_4,dustdrydep_5, &
12 depvelocity, &
13 ids,ide, jds,jde, kds,kde, &
14 ims,ime, jms,jme, kms,kme, &
15 its,ite, jts,jte, kts,kte )
17 USE module_model_constants
18 USE module_configure
19 USE module_state_description
21 IMPLICIT NONE
23 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
24 ims,ime, jms,jme, kms,kme, &
25 its,ite, jts,jte, kts,kte,numgas
26 INTEGER,DIMENSION( ims:ime , jms:jme ) , &
27 INTENT(IN ) :: &
28 ivgtyp
30 REAL, INTENT(IN ) :: dtstep
31 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
32 INTENT(IN ) :: moist
33 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
34 INTENT(INOUT ) :: chem
35 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
36 INTENT(IN ) :: &
37 t_phy, &
38 p_phy, &
39 dz8w, &
40 t8w,p8w,rho_phy
41 REAL, DIMENSION( its:ite, jts:jte, num_chem ), &
42 INTENT(INOUT) :: ddvel
43 REAL, DIMENSION( ims:ime , jms:jme ) , &
44 INTENT(IN) :: &
45 tsk, &
46 vegfra, &
47 pbl, &
48 ust, &
49 xlat, &
50 xlong, &
51 rmol,xland,znt,hfx
52 REAL, DIMENSION( its:ite, jts:jte ), &
53 INTENT(IN) :: aer_res
54 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
55 dustdrydep_1, dustdrydep_2, dustdrydep_3, &
56 dustdrydep_4, dustdrydep_5, depvelocity
58 !! .. Local Scalars ..
59 INTEGER :: iland, iprt, iseason, jce, jcs, &
60 n, nr, ipr, jpr, nvr, &
61 idrydep_onoff,imx,jmx,lmx
63 integer :: ii,jj,kk,i,j,k,nv
64 integer,dimension (1,1) :: ilwi,ireg
66 REAL :: clwchem, dvfog, dvpart, &
67 rad, rhchem, ta, vegfrac, z1,zntt
68 ! real*8, DIMENSION (1,1,1,3) :: erodin
69 ! real*8, DIMENSION (5) :: tc,bems
70 ! real*8, dimension (1,1) :: z0,w10m,gwet,airden,airmas,delz_sfc,hflux,ts,pblz,ustar,ps
71 real*8, dimension (1,1) :: z0,airden,delz_sfc,hflux,ts,pblz,ustar,ps
72 REAL*8 :: dvel(1,1), drydf(1,1)
73 LOGICAL :: highnh3, rainflag, vegflag, wetflag
75 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
77 do nv=numgas+1,num_chem
78 do j=jts,jte
79 do i=its,ite
80 ddvel(i,j,nv) = 0.
81 enddo
82 enddo
83 enddo
85 imx = 1
86 jmx = 1
87 lmx = 1
88 do j = max(jds+1,jts),min(jde-1,jte)
89 do i = max(ids+1,its),min(ide-1,ite)
90 ipr = 0
91 dvel(1,1) = 0._8
92 if( xland(i,j) > 1.5 ) then
93 ilwi(1,1) = 1
94 else
95 ilwi(1,1) = 0
96 endif
97 ! for aerosols, ii=1 or ii=2
98 ii = 1
99 ! if(ivgtyp(i,j).eq.19.or.ivgtyp(i,j).eq.23)ii=1
100 ireg(1,1) = 1
101 airden(1,1) = real( rho_phy(i,kts,j),kind=8 )
102 delz_sfc(1,1) = real( dz8w(i,kts,j),kind=8 )
103 ustar(1,1) = max( 1.e-1_8,real(ust(i,j),kind=8) )
104 hflux(1,1) = real( hfx(i,j),kind=8 )
105 pblz(1,1) = real( pbl(i,j),kind=8 )
106 ps(1,1) = real(p8w(i,kts,j),kind=8)*.01_8
107 z0(1,1) = real( znt(i,j),kind=8 )
108 ts(1,1) = real( tsk(i,j),kind=8 )
109 ! if(i.eq.23.and.j.eq.74)ipr=1
110 call depvel_gocart(config_flags,ipr,ii,imx,jmx,lmx,&
111 airden, delz_sfc, pblz, ts, ustar, hflux, ilwi, &
112 ps, z0, dvel, drydf,g,rmol(i,j),aer_res(i,j))
113 do nv = p_p25,num_chem
114 ddvel(i,j,nv) = real( dvel(1,1),kind=4 )
115 enddo
116 ddvel(i,j,p_sulf) = real( dvel(1,1),kind=4 )
117 ddvel(i,j,p_msa) = real( dvel(1,1),kind=4 )
119 depvelocity(i,j) = ddvel(i,j,p_dust_5)
120 !drydep [kg/m2] = drydep [kg/m2]+1.e-9*dt[s]*dvel [m/s] * chem [ug/kg] * airden [kg/m3]
121 dustdrydep_1(i,j)=dustdrydep_1(i,j)+1.e-9*dtstep*dvel(1,1)*chem(i,1,j,p_dust_1)*airden(1,1)
122 dustdrydep_2(i,j)=dustdrydep_2(i,j)+1.e-9*dtstep*dvel(1,1)*chem(i,1,j,p_dust_2)*airden(1,1)
123 dustdrydep_3(i,j)=dustdrydep_3(i,j)+1.e-9*dtstep*dvel(1,1)*chem(i,1,j,p_dust_3)*airden(1,1)
124 dustdrydep_4(i,j)=dustdrydep_4(i,j)+1.e-9*dtstep*dvel(1,1)*chem(i,1,j,p_dust_4)*airden(1,1)
125 dustdrydep_5(i,j)=dustdrydep_5(i,j)+1.e-9*dtstep*dvel(1,1)*chem(i,1,j,p_dust_5)*airden(1,1)
127 enddo
128 enddo
130 end subroutine gocart_drydep_driver
132 SUBROUTINE depvel_gocart( config_flags,ipr,ii,imx,jmx,lmx, &
133 airden, delz_sfc, pblz, ts, ustar, hflux, ilwi, &
134 ps, z0, dvel, drydf,g0,rmol,aer_res)
136 ! ****************************************************************************
137 ! * *
138 ! * Calculate dry deposition velocity. *
139 ! * *
140 ! * Input variables: *
141 ! * AEROSOL(k) - Logical, T = aerosol species, F = gas species *
142 ! * IREG(i,j) - # of landtypes in grid square *
143 ! * ILAND(i,j,ldt) - Land type ID for element ldt =1,IREG(i,j) *
144 ! * IUSE(i,j,ldt) - Fraction of gridbox area occupied by land type *
145 ! * element ldt *
146 ! * USTAR(i,j) - Friction velocity (m s-1) *
147 ! * DELZ_SFC(i,j) - Thickness of layer above surface *
148 ! * PBLZ(i,j) - Mixing depth (m) *
149 ! * Z0(i,j) - Roughness height (m) *
150 ! * *
151 ! * Determined in this subroutine (local): *
152 ! * OBK - Monin-Obukhov length (m): set to 1.E5 m under *
153 ! * neutral conditions *
154 ! * Rs(ldt) - Bulk surface resistance(s m-1) for species k to *
155 ! * surface ldt *
156 ! * Ra - Aerodynamic resistance. *
157 ! * Rb - Sublayer resistance. *
158 ! * Rs - Surface resistance. *
159 ! * Rttl - Total deposition resistance (s m-1) for species k *
160 ! * Rttl(k) = Ra + Rb + Rs. *
161 ! * *
162 ! * Returned: *
163 ! * DVEL(i,j,k) - Deposition velocity (m s-1) of species k *
164 ! * DRYDf(i,j,k) - Deposition frequency (s-1) of species k, *
165 ! * = DVEL / DELZ_SFC *
166 ! * *
167 ! ****************************************************************************
169 USE module_configure
171 IMPLICIT NONE
173 REAL*8, INTENT(IN) :: airden(imx,jmx), delz_sfc(imx,jmx)
174 REAL*8, INTENT(IN) :: hflux(imx,jmx), ts(imx,jmx)
175 REAL*8, INTENT(IN) :: ustar(imx,jmx), pblz(imx,jmx)
176 REAL*8, INTENT(IN) :: ps(imx,jmx)
177 INTEGER, INTENT(IN) :: ilwi(imx,jmx)
178 INTEGER, INTENT(IN) :: imx,jmx,lmx
179 REAL*8, INTENT(IN) :: z0(imx,jmx)
180 REAL, INTENT(IN) :: g0,rmol,aer_res
181 REAL*8, INTENT(OUT) :: dvel(imx,jmx), drydf(imx,jmx)
183 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
185 REAL*8 :: obk, vds, czh, rttl, frac, logmfrac, psi_h, cz, eps
186 REAL*8 :: vd, ra, rb, rs
187 INTEGER :: ipr,i, j, k, ldt, iolson, ii
188 CHARACTER(LEN=50) :: msg
189 REAL*8 :: prss, tempk, tempc, xnu, ckustr, reyno, aird, diam, xm, z
190 REAL*8 :: frpath, speed, dg, dw, rt
191 REAL*8 :: rad0, rix, gfact, gfaci, rdc, rixx, rluxx, rgsx, rclx
192 REAL*8 :: dtmp1, dtmp2, dtmp3, dtmp4
193 REAL*8 :: biofit,vk
195 ! executable statements
196 j_loop: DO j = 1,jmx
197 i_loop: DO i = 1,imx
198 vk = .4_8
199 vd = 0._8
200 ra = 0._8
201 rb = 0._8 ! only required for gases (SO2)
202 rs = 0.0_8
204 ! ****************************************************************************
205 ! * Compute the the Monin-Obhukov length. *
206 ! * The direct computation of the Monin-Obhukov length is: *
207 ! * *
208 ! * - Air density * Cp * T(surface air) * Ustar^3 *
209 ! * OBK = ---------------------------------------------- *
210 ! * vK * g * Sensible Heat flux *
211 ! * *
212 ! * Cp = 1000 J/kg/K = specific heat at constant pressure *
213 ! * vK = 0.4 = von Karman's constant *
214 ! ****************************************************************************
216 IF (abs(hflux(i,j)) <= 1.e-5_8) THEN
217 obk = 1.0E5_8
218 ELSE
219 ! MINVAL(hflux), MINVAL(airden), MINVAL(ustar) =??
220 obk = -airden(i,j) * 1000.0_8 * ts(i,j) * (ustar(i,j))**3 &
221 / (vk * real(g0,kind=8) * hflux(i,j))
222 ! -- debug:
223 IF ( obk == 0.0_8 ) WRITE(*,211) obk, i, j
224 211 FORMAT(1X,'OBK=', E11.2, 1X,' i,j = ', 2I4)
226 END IF
227 ! write(0,*)1./obk,rmol
229 if(rmol.ne.0.)then
230 obk=1._8/real(rmol,kind=8)
231 else
232 obk=1.e5_8
233 endif
235 ! cz = delz_sfc(i,j) / 2.0_8 ! center of the grid box above surface
236 cz = 2._8
238 ! ****************************************************************************
239 ! * (1) Aerosodynamic resistance Ra and sublayer resistance Rb. *
240 ! * *
241 ! * The Reynolds number REYNO diagnoses whether a surface is *
242 ! * aerodynamically rough (REYNO > 10) or smooth. Surface is *
243 ! * rough in all cases except over water with low wind speeds. *
244 ! * *
245 ! * For gas species over land and ice (REYNO >= 10) and for aerosol *
246 ! * species for all surfaces: *
247 ! * *
248 ! * Ra = 1./VT (VT from GEOS Kzz at L=1, m/s). *
249 ! * *
250 ! * The following equations are from Walcek et al, 1986: *
251 ! * *
252 ! * For gas species when REYNO < 10 (smooth), Ra and Rb are combined *
253 ! * as Ra: *
254 ! * *
255 ! * Ra = { ln(ku* z1/Dg) - Sh } / ku* eq.(13) *
256 ! * *
257 ! * where z1 is the altitude at the center of the lowest model layer *
258 ! * (CZ); *
259 ! * Sh is a stability correction function; *
260 ! * k is the von Karman constant (0.4, vK); *
261 ! * u* is the friction velocity (USTAR). *
262 ! * *
263 ! * Sh is computed as a function of z1 and L eq ( 4) and (5)): *
264 ! * *
265 ! * 0 < z1/L <= 1: Sh = -5 * z1/L *
266 ! * z1/L < 0: Sh = exp{ 0.598 + 0.39*ln(E) - 0.09(ln(E))^2 } *
267 ! * where E = min(1,-z1/L) (Balkanski, thesis). *
268 ! * *
269 ! * For gas species when REYNO >= 10, *
270 ! * *
271 ! * Rb = 2/ku* (Dair/Dg)**(2/3) eq.(12) *
272 ! * where Dg is the gas diffusivity, and *
273 ! * Dair is the air diffusivity. *
274 ! * *
275 ! * For aerosol species, Rb is combined with surface resistance as Rs. *
276 ! * *
277 ! ****************************************************************************
279 frac = cz / obk
280 IF (frac > 1.0_8) frac = 1.0_8
281 IF (frac > 0.0_8 .AND. frac <= 1.0_8) THEN
282 psi_h = -5.0_8*frac
283 ELSE IF (frac < 0.0_8) THEN
284 eps = MIN(1.0_8, -frac)
285 logmfrac = LOG(eps)
286 psi_h = EXP( 0.598_8 + 0.39_8 * logmfrac - 0.09_8 * (logmfrac)**2 )
287 END IF
288 !--------------------------------------------------------------
289 ! Aerosol species, Rs here is the combination of Rb and Rs.
291 ra = (LOG(cz/z0(i,j)) - psi_h) / (vk*ustar(i,j))
293 vds = 0.002_8*ustar(i,j)
294 IF (obk < 0.0_8) vds = vds * (1.0_8+(-300.0_8/obk)**0.6667_8)
296 czh = pblz(i,j)/obk
297 IF (czh < -30.0_8) vds = 0.0009_8*ustar(i,j)*(-czh)**0.6667_8
298 if(ipr.eq.1) write(0,*)ra,aer_res,vds
300 IF( config_flags%chem_opt /= CHEM_VASH .and. &
301 config_flags%chem_opt /= DUST )THEN
302 ra = real(aer_res,kind=8)
303 ENDIF
305 ! --Set Vds to be less than VDSMAX (entry in input file divided --
306 ! by 1.E4). VDSMAX is taken from Table 2 of Walcek et al. [1986].
307 ! Invert to get corresponding R
308 ! if(ii.eq.1) then
309 ! rs=1.0_8/MIN(vds,2.0e-2_8)
310 ! else
311 rs=1.0_8/MIN(vds,2.0e-3_8)
312 ! endif
315 ! ------ Set max and min values for bulk surface resistances ------
317 rs= MAX(1.0_8, MIN(rs, 9.9990e+3_8))
319 ! ****************************************************************************
320 ! * *
321 ! * Compute dry deposition velocity. *
322 ! * *
323 ! * IUSE is the fraction of the grid square occupied by surface ldt in *
324 ! * units of per mil (IUSE=500 -> 50% of the grid square). Add the *
325 ! * contribution of surface type ldt to the deposition velocity; this is *
326 ! * a loop over all surface types in the gridbox. *
327 ! * *
328 ! * Total resistance = Ra + Rb + Rs.
329 ! * *
330 ! ****************************************************************************
332 rttl = ra + rb + rs
333 vd = vd + 1._8/rttl
335 ! ------ Load array DVEL ------
336 if(ipr.eq.1) write(0,*)rs,ra,rb,vd
337 dvel(i,j) = vd !* 1.2
339 ! -- Set a minimum value for DVEL
340 ! MIN(VdSO2) = 2.0e-3 m/s over ice
341 ! = 3.0e-3 m/s over land
342 ! MIN(vd_aerosol) = 1.0e-4 m/s
344 IF (dvel(i,j) < 1.0E-4_8) dvel(i,j) = 1.0E-4_8
345 drydf(i,j) = dvel(i,j) / delz_sfc(i,j)
347 END DO i_loop
348 END DO j_loop
350 END SUBROUTINE depvel_gocart
354 END MODULE module_gocart_drydep