Line data Source code
1 : module dust_sediment_mod
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : !
6 : ! Contains routines to compute tendencies from sedimentation of dust
7 : !
8 : ! Author: Phil Rasch
9 : !
10 : !---------------------------------------------------------------------------------
11 :
12 : use shr_kind_mod, only: r8=>shr_kind_r8
13 : use ppgrid, only: pcols, pver, pverp
14 : use physconst, only: gravit, rair
15 : use cam_logfile, only: iulog
16 : use cam_abortutils, only: endrun
17 :
18 : private
19 : public :: dust_sediment_vel, dust_sediment_tend
20 :
21 :
22 : real (r8), parameter :: vland = 2.8_r8 ! dust fall velocity over land (cm/s)
23 : real (r8), parameter :: vocean = 1.5_r8 ! dust fall velocity over ocean (cm/s)
24 : real (r8), parameter :: mxsedfac = 0.99_r8 ! maximum sedimentation flux factor
25 :
26 : contains
27 :
28 : !===============================================================================
29 0 : subroutine dust_sediment_vel (ncol, &
30 : icefrac , landfrac, ocnfrac , pmid , pdel , t , &
31 : dustmr , pvdust )
32 :
33 : !----------------------------------------------------------------------
34 :
35 : ! Compute gravitational sedimentation velocities for dust
36 :
37 : implicit none
38 :
39 : ! Arguments
40 : integer, intent(in) :: ncol ! number of colums to process
41 :
42 : real(r8), intent(in) :: icefrac (pcols) ! sea ice fraction (fraction)
43 : real(r8), intent(in) :: landfrac(pcols) ! land fraction (fraction)
44 : real(r8), intent(in) :: ocnfrac (pcols) ! ocean fraction (fraction)
45 : real(r8), intent(in) :: pmid (pcols,pver) ! pressure of midpoint levels (Pa)
46 : real(r8), intent(in) :: pdel (pcols,pver) ! pressure diff across layer (Pa)
47 : real(r8), intent(in) :: t (pcols,pver) ! temperature (K)
48 : real(r8), intent(in) :: dustmr(pcols,pver) ! dust (kg/kg)
49 :
50 : real(r8), intent(out) :: pvdust (pcols,pverp) ! vertical velocity of dust (Pa/s)
51 : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
52 :
53 : ! Local variables
54 : real (r8) :: rho(pcols,pver) ! air density in kg/m3
55 : real (r8) :: vfall(pcols) ! settling velocity of dust particles (m/s)
56 :
57 : integer i,k
58 :
59 : real (r8) :: lbound, ac, bc, cc
60 :
61 : !-----------------------------------------------------------------------
62 : !--------------------- dust fall velocity ----------------------------
63 : !-----------------------------------------------------------------------
64 :
65 0 : do k = 1,pver
66 0 : do i = 1,ncol
67 :
68 : ! merge the dust fall velocities for land and ocean (cm/s)
69 : ! SHOULD ALSO ACCOUNT FOR ICEFRAC
70 0 : vfall(i) = vland*landfrac(i) + vocean*(1._r8-landfrac(i))
71 : !! vfall(i) = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)
72 :
73 : ! fall velocity (assume positive downward)
74 0 : pvdust(i,k+1) = vfall(i)
75 : end do
76 : end do
77 :
78 0 : return
79 : end subroutine dust_sediment_vel
80 :
81 :
82 : !===============================================================================
83 0 : subroutine dust_sediment_tend ( &
84 : ncol, dtime, pint, pmid, pdel, t, &
85 : dustmr ,pvdust, dusttend, sfdust )
86 :
87 : !----------------------------------------------------------------------
88 : ! Apply Particle Gravitational Sedimentation
89 : !----------------------------------------------------------------------
90 :
91 : implicit none
92 :
93 : ! Arguments
94 : integer, intent(in) :: ncol ! number of colums to process
95 :
96 : real(r8), intent(in) :: dtime ! time step
97 : real(r8), intent(in) :: pint (pcols,pverp) ! interfaces pressure (Pa)
98 : real(r8), intent(in) :: pmid (pcols,pver) ! midpoint pressures (Pa)
99 : real(r8), intent(in) :: pdel (pcols,pver) ! pressure diff across layer (Pa)
100 : real(r8), intent(in) :: t (pcols,pver) ! temperature (K)
101 : real(r8), intent(in) :: dustmr(pcols,pver) ! dust (kg/kg)
102 : real(r8), intent(in) :: pvdust (pcols,pverp) ! vertical velocity of dust drops (Pa/s)
103 : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
104 :
105 : real(r8), intent(out) :: dusttend(pcols,pver) ! dust tend
106 : real(r8), intent(out) :: sfdust (pcols) ! surface flux of dust (rain, kg/m/s)
107 :
108 : ! Local variables
109 : real(r8) :: fxdust(pcols,pverp) ! fluxes at the interfaces, dust (positive = down)
110 :
111 : integer :: i,k
112 : !----------------------------------------------------------------------
113 :
114 : ! initialize variables
115 0 : fxdust (:ncol,:) = 0._r8 ! flux at interfaces (dust)
116 0 : dusttend(:ncol,:) = 0._r8 ! tend (dust)
117 0 : sfdust(:ncol) = 0._r8 ! sedimentation flux out bot of column (dust)
118 :
119 : ! fluxes at interior points
120 0 : call getflx(ncol, pint, dustmr, pvdust, dtime, fxdust)
121 :
122 : ! calculate fluxes at boundaries
123 0 : do i = 1,ncol
124 0 : fxdust(i,1) = 0
125 : ! surface flux by upstream scheme
126 0 : fxdust(i,pverp) = dustmr(i,pver) * pvdust(i,pverp) * dtime
127 : end do
128 :
129 : ! filter out any negative fluxes from the getflx routine
130 0 : do k = 2,pver
131 0 : fxdust(:ncol,k) = max(0._r8, fxdust(:ncol,k))
132 : end do
133 :
134 : ! Limit the flux out of the bottom of each cell to the water content in each phase.
135 : ! Apply mxsedfac to prevent generating very small negative cloud water/ice
136 : ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
137 : ! ***Should we include the flux in the top, to allow for thin surface layers?
138 : ! ***Requires simple treatment of cloud overlap, already included below.
139 0 : do k = 1,pver
140 0 : do i = 1,ncol
141 0 : fxdust(i,k+1) = min( fxdust(i,k+1), mxsedfac * dustmr(i,k) * pdel(i,k) )
142 : !!$ fxdust(i,k+1) = min( fxdust(i,k+1), dustmr(i,k) * pdel(i,k) + fxdust(i,k))
143 : end do
144 : end do
145 :
146 : ! Now calculate the tendencies
147 0 : do k = 1,pver
148 0 : do i = 1,ncol
149 : ! net flux into cloud changes cloud dust/ice (all flux is out of cloud)
150 0 : dusttend(i,k) = (fxdust(i,k) - fxdust(i,k+1)) / (dtime * pdel(i,k))
151 : end do
152 : end do
153 :
154 : ! convert flux out the bottom to mass units Pa -> kg/m2/s
155 0 : sfdust(:ncol) = fxdust(:ncol,pverp) / (dtime*gravit)
156 :
157 0 : return
158 : end subroutine dust_sediment_tend
159 :
160 : !===============================================================================
161 0 : subroutine getflx(ncol, xw, phi, vel, deltat, flux)
162 :
163 : !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
164 : !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
165 : !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
166 : !.........phi1......phi2.......phi3.....phi4.......phi5.......
167 :
168 :
169 : implicit none
170 :
171 : integer ncol ! number of colums to process
172 :
173 : integer i
174 : integer k
175 :
176 : real (r8) vel(pcols,pverp)
177 : real (r8) flux(pcols,pverp)
178 : real (r8) xw(pcols,pverp)
179 : real (r8) psi(pcols,pverp)
180 : real (r8) phi(pcols,pverp-1)
181 : real (r8) fdot(pcols,pverp)
182 : real (r8) xx(pcols)
183 : real (r8) fxdot(pcols)
184 : real (r8) fxdd(pcols)
185 :
186 : real (r8) psistar(pcols)
187 : real (r8) deltat
188 :
189 : real (r8) xxk(pcols,pver)
190 :
191 0 : do i = 1,ncol
192 : ! integral of phi
193 0 : psi(i,1) = 0._r8
194 : ! fluxes at boundaries
195 0 : flux(i,1) = 0
196 0 : flux(i,pverp) = 0._r8
197 : end do
198 :
199 : ! integral function
200 0 : do k = 2,pverp
201 0 : do i = 1,ncol
202 0 : psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
203 : end do
204 : end do
205 :
206 :
207 : ! calculate the derivatives for the interpolating polynomial
208 0 : call cfdotmc_pro (ncol, xw, psi, fdot)
209 :
210 : ! NEW WAY
211 : ! calculate fluxes at interior pts
212 0 : do k = 2,pver
213 0 : do i = 1,ncol
214 0 : xxk(i,k) = xw(i,k)-vel(i,k)*deltat
215 : end do
216 : end do
217 0 : do k = 2,pver
218 0 : call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
219 0 : do i = 1,ncol
220 0 : flux(i,k) = (psi(i,k)-psistar(i))
221 : end do
222 : end do
223 :
224 :
225 0 : return
226 : end subroutine getflx
227 :
228 :
229 :
230 : !##############################################################################
231 :
232 0 : subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar)
233 :
234 :
235 : implicit none
236 :
237 : ! input
238 : integer ncol ! number of colums to process
239 :
240 : real (r8) x(pcols, pverp)
241 : real (r8) f(pcols, pverp)
242 : real (r8) fdot(pcols, pverp)
243 : real (r8) xin(pcols)
244 :
245 : ! output
246 : real (r8) fxdot(pcols)
247 : real (r8) fxdd(pcols)
248 : real (r8) psistar(pcols)
249 :
250 : integer i
251 : integer k
252 : integer intz(pcols)
253 : real (r8) dx
254 : real (r8) s
255 : real (r8) c2
256 : real (r8) c3
257 : real (r8) xx
258 : real (r8) xinf
259 : real (r8) psi1, psi2, psi3, psim
260 : real (r8) cfint
261 : real (r8) cfnew
262 : real (r8) xins(pcols)
263 :
264 : ! the minmod function
265 : real (r8) a, b, c
266 : real (r8) minmod
267 : real (r8) medan
268 : minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
269 : medan(a,b,c) = a + minmod(b-a,c-a)
270 :
271 0 : do i = 1,ncol
272 0 : xins(i) = medan(x(i,1), xin(i), x(i,pverp))
273 0 : intz(i) = 0
274 : end do
275 :
276 : ! first find the interval
277 0 : do k = 1,pverp-1
278 0 : do i = 1,ncol
279 0 : if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0._r8) then
280 0 : intz(i) = k
281 : endif
282 : end do
283 : end do
284 :
285 0 : do i = 1,ncol
286 0 : if (intz(i).eq.0) then
287 0 : write(iulog,*) ' interval was not found for col i ', i
288 0 : call endrun('DUST_SEDIMENT_MOD:cfint2 -- interval was not found ')
289 : endif
290 : end do
291 :
292 : ! now interpolate
293 0 : do i = 1,ncol
294 0 : k = intz(i)
295 0 : dx = (x(i,k+1)-x(i,k))
296 0 : s = (f(i,k+1)-f(i,k))/dx
297 0 : c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
298 0 : c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
299 0 : xx = (xins(i)-x(i,k))
300 0 : fxdot(i) = (3*c3*xx + 2*c2)*xx + fdot(i,k)
301 0 : fxdd(i) = 6*c3*xx + 2*c2
302 0 : cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)
303 :
304 : ! limit the interpolant
305 0 : psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
306 0 : if (k.eq.1) then
307 0 : psi2 = f(i,1)
308 : else
309 0 : psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
310 : endif
311 0 : if (k+1.eq.pverp) then
312 0 : psi3 = f(i,pverp)
313 : else
314 0 : psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
315 : endif
316 0 : psim = medan(psi1, psi2, psi3)
317 0 : cfnew = medan(cfint, psi1, psim)
318 : if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36_r8) .gt..03_r8) then
319 : ! CHANGE THIS BACK LATER!!!
320 : ! $ .gt..1) then
321 :
322 :
323 : ! UNCOMMENT THIS LATER!!!
324 : ! write(iulog,*) ' cfint2 limiting important ', cfint, cfnew
325 :
326 :
327 : endif
328 0 : psistar(i) = cfnew
329 : end do
330 :
331 0 : return
332 : end subroutine cfint2
333 :
334 :
335 :
336 : !##############################################################################
337 :
338 0 : subroutine cfdotmc_pro (ncol, x, f, fdot)
339 :
340 : ! prototype version; eventually replace with final SPITFIRE scheme
341 :
342 : ! calculate the derivative for the interpolating polynomial
343 : ! multi column version
344 :
345 :
346 : implicit none
347 :
348 : ! input
349 : integer ncol ! number of colums to process
350 :
351 : real (r8) x(pcols, pverp)
352 : real (r8) f(pcols, pverp)
353 : ! output
354 : real (r8) fdot(pcols, pverp) ! derivative at nodes
355 :
356 : ! assumed variable distribution
357 : ! x1.......x2.......x3.......x4.......x5.......x6 1,pverp points
358 : ! f1.......f2.......f3.......f4.......f5.......f6 1,pverp points
359 : ! ...sh1.......sh2......sh3......sh4......sh5.... 1,pver points
360 : ! .........d2.......d3.......d4.......d5......... 2,pver points
361 : ! .........s2.......s3.......s4.......s5......... 2,pver points
362 : ! .............dh2......dh3......dh4............. 2,pver-1 points
363 : ! .............eh2......eh3......eh4............. 2,pver-1 points
364 : ! ..................e3.......e4.................. 3,pver-1 points
365 : ! .................ppl3......ppl4................ 3,pver-1 points
366 : ! .................ppr3......ppr4................ 3,pver-1 points
367 : ! .................t3........t4.................. 3,pver-1 points
368 : ! ................fdot3.....fdot4................ 3,pver-1 points
369 :
370 :
371 : ! work variables
372 :
373 :
374 : integer i
375 : integer k
376 :
377 : real (r8) a ! work var
378 : real (r8) b ! work var
379 : real (r8) c ! work var
380 : real (r8) s(pcols,pverp) ! first divided differences at nodes
381 : real (r8) sh(pcols,pverp) ! first divided differences between nodes
382 : real (r8) d(pcols,pverp) ! second divided differences at nodes
383 : real (r8) dh(pcols,pverp) ! second divided differences between nodes
384 : real (r8) e(pcols,pverp) ! third divided differences at nodes
385 : real (r8) eh(pcols,pverp) ! third divided differences between nodes
386 : real (r8) pp ! p prime
387 : real (r8) ppl(pcols,pverp) ! p prime on left
388 : real (r8) ppr(pcols,pverp) ! p prime on right
389 : real (r8) qpl
390 : real (r8) qpr
391 : real (r8) ttt
392 : real (r8) t
393 : real (r8) tmin
394 : real (r8) tmax
395 : real (r8) delxh(pcols,pverp)
396 :
397 :
398 : ! the minmod function
399 : real (r8) minmod
400 : real (r8) medan
401 : minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
402 : medan(a,b,c) = a + minmod(b-a,c-a)
403 :
404 0 : do k = 1,pver
405 :
406 :
407 : ! first divided differences between nodes
408 0 : do i = 1, ncol
409 0 : delxh(i,k) = (x(i,k+1)-x(i,k))
410 0 : sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
411 : end do
412 :
413 : ! first and second divided differences at nodes
414 0 : if (k.ge.2) then
415 0 : do i = 1,ncol
416 0 : d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
417 0 : s(i,k) = minmod(sh(i,k),sh(i,k-1))
418 : end do
419 : endif
420 : end do
421 :
422 : ! second and third divided diffs between nodes
423 0 : do k = 2,pver-1
424 0 : do i = 1, ncol
425 0 : eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
426 0 : dh(i,k) = minmod(d(i,k),d(i,k+1))
427 : end do
428 : end do
429 :
430 : ! treat the boundaries
431 0 : do i = 1,ncol
432 0 : e(i,2) = eh(i,2)
433 0 : e(i,pver) = eh(i,pver-1)
434 : ! outside level
435 : fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1) &
436 0 : - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
437 0 : fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
438 : fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver) &
439 0 : + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
440 0 : fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
441 : ! one in from boundary
442 0 : fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
443 0 : fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
444 : fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver) &
445 0 : - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
446 0 : fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
447 : end do
448 :
449 :
450 0 : do k = 3,pver-1
451 0 : do i = 1,ncol
452 0 : e(i,k) = minmod(eh(i,k),eh(i,k-1))
453 : end do
454 : end do
455 :
456 :
457 :
458 0 : do k = 3,pver-1
459 :
460 0 : do i = 1,ncol
461 :
462 : ! p prime at k-0.5
463 0 : ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)
464 : ! p prime at k+0.5
465 0 : ppr(i,k)=sh(i,k) - dh(i,k) *delxh(i,k)
466 :
467 0 : t = minmod(ppl(i,k),ppr(i,k))
468 :
469 : ! derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
470 0 : pp = sh(i,k-1) + d(i,k)*delxh(i,k-1)
471 :
472 : ! quartic estimate of fdot
473 : fdot(i,k) = pp &
474 : - delxh(i,k-1)*delxh(i,k) &
475 0 : *( eh(i,k-1)*(x(i,k+2)-x(i,k )) &
476 0 : + eh(i,k )*(x(i,k )-x(i,k-2)) &
477 0 : )/(x(i,k+2)-x(i,k-2))
478 :
479 : ! now limit it
480 : qpl = sh(i,k-1) &
481 : + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
482 0 : d(i,k) -e(i,k)*delxh(i,k))
483 : qpr = sh(i,k) &
484 : + delxh(i,k )*minmod(d(i,k) +e(i,k)*delxh(i,k-1), &
485 0 : d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))
486 :
487 0 : fdot(i,k) = medan(fdot(i,k), qpl, qpr)
488 :
489 0 : ttt = minmod(qpl, qpr)
490 0 : tmin = min(0._r8,3*s(i,k),1.5_r8*t,ttt)
491 0 : tmax = max(0._r8,3*s(i,k),1.5_r8*t,ttt)
492 :
493 0 : fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))
494 :
495 : end do
496 :
497 : end do
498 :
499 0 : return
500 : end subroutine cfdotmc_pro
501 : end module dust_sediment_mod
|
