Line data Source code
1 : module molec_diff
2 :
3 : !------------------------------------------------------------------------------------------------- !
4 : ! Module to compute molecular diffusivity for various constituents !
5 : ! !
6 : ! Public interfaces : !
7 : ! !
8 : ! init_molec_diff Initializes time independent coefficients !
9 : ! init_timestep_molec_diff Time-step initialization for molecular diffusivity !
10 : ! compute_molec_diff Computes constituent-independent terms for moleculuar diffusivity !
11 : ! vd_lu_qdecomp Computes constituent-dependent terms for moleculuar diffusivity and !
12 : ! updates terms in the triadiagonal matrix used for the implicit !
13 : ! solution of the diffusion equation !
14 : ! !
15 : !---------------------------Code history---------------------------------------------------------- !
16 : ! Modularized : J. McCaa, September 2004 !
17 : ! Lastly Arranged : S. Park, January. 2010 !
18 : ! M. Mills, November 2011
19 : !------------------------------------------------------------------------------------------------- !
20 :
21 : use perf_mod
22 : use air_composition, only: mbarv
23 : use phys_control, only: waccmx_is !WACCM-X runtime switch
24 :
25 : implicit none
26 : private
27 : save
28 :
29 : public init_molec_diff
30 : public compute_molec_diff
31 : public vd_lu_qdecomp
32 :
33 : ! ---------- !
34 : ! Parameters !
35 : ! ---------- !
36 :
37 : integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
38 :
39 : real(r8), parameter :: km_fac = 3.55E-7_r8 ! Molecular viscosity constant [ unit ? ]
40 : real(r8), parameter :: pwr = 2._r8/3._r8 ! Exponentiation factor [ no unit ]
41 : real(r8), parameter :: d0 = 1.52E20_r8 ! Diffusion factor [ m-1 s-1 ] molec sqrt(kg/kmol/K) [ unit ? ]
42 : ! Aerononmy, Part B, Banks and Kockarts (1973), p39
43 : ! Note text cites 1.52E18 cm-1 ...
44 :
45 : real(r8) :: mw_dry ! Molecular weight of dry air
46 : real(r8) :: n_avog ! Avogadro's number [ molec/kmol ]
47 :
48 : real(r8), allocatable :: mw_fac(:) ! sqrt(1/M_q + 1/M_d) in constituent diffusivity [ unit ? ]
49 : real(r8), allocatable :: alphath(:) ! Thermal diffusion factor, -0.38 for H, 0 for others
50 :
51 : logical :: waccmx_mode = .false.
52 :
53 : contains
54 :
55 : !============================================================================ !
56 : ! !
57 : !============================================================================ !
58 :
59 0 : subroutine init_molec_diff( kind, ncnst, mw_dry_in, n_avog_in, &
60 : errstring)
61 :
62 : use constituents, only: cnst_mw, cnst_get_ind
63 :
64 : integer, intent(in) :: kind ! Kind of reals being passed in
65 : integer, intent(in) :: ncnst ! Number of constituents
66 : real(r8), intent(in) :: mw_dry_in ! Molecular weight of dry air
67 : real(r8), intent(in) :: n_avog_in ! Avogadro's number [ molec/kmol ]
68 :
69 : character(len=*), intent(out) :: errstring
70 :
71 : ! Local
72 :
73 : integer :: m ! Constituent index
74 : integer :: indx_H ! Constituent index for H
75 : integer :: ierr ! Allocate error check
76 :
77 0 : errstring = ' '
78 :
79 0 : mw_dry = mw_dry_in
80 0 : n_avog = n_avog_in
81 :
82 0 : if( kind /= r8 ) then
83 0 : errstring = 'inconsistent KIND of reals passed to init_molec_diff'
84 0 : return
85 : end if
86 :
87 : ! Determine whether WACCM-X is on.
88 0 : waccmx_mode = waccmx_is('ionosphere') .or. waccmx_is('neutral')
89 :
90 : ! Molecular weight factor in constitutent diffusivity
91 : ! ***** FAKE THIS FOR NOW USING MOLECULAR WEIGHT OF DRY AIR FOR ALL TRACERS ****
92 :
93 0 : allocate(mw_fac(ncnst))
94 0 : do m = 1, ncnst
95 0 : mw_fac(m) = d0 * mw_dry * sqrt(1._r8/mw_dry + 1._r8/cnst_mw(m)) / n_avog
96 : end do
97 :
98 : !--------------------------------------------------------------------------------------------
99 : ! For WACCM-X, get H data index and initialize thermal diffusion coefficient
100 : !--------------------------------------------------------------------------------------------
101 0 : if ( waccmx_mode ) then
102 :
103 0 : call cnst_get_ind('H', indx_H)
104 :
105 0 : allocate(alphath(ncnst), stat=ierr)
106 0 : if ( ierr /= 0 ) then
107 0 : errstring = 'allocate failed in init_molec_diff'
108 0 : return
109 : end if
110 0 : alphath(:ncnst) = 0._r8
111 0 : alphath(indx_H) = -0.38_r8
112 :
113 : endif
114 :
115 0 : end subroutine init_molec_diff
116 :
117 : !============================================================================ !
118 : ! !
119 : !============================================================================ !
120 :
121 0 : subroutine compute_molec_diff(lchnk, pcols, pver, ncnst, ncol, &
122 0 : kvm, kvt, tint, rhoi, kq_scal, cnst_mw, &
123 0 : mw_fac_out, nbot_molec)
124 :
125 : use cam_thermo, only: kmvis, kmcnd
126 : use air_composition, only: cpairv
127 :
128 : ! --------------------- !
129 : ! Input-Output Argument !
130 : ! --------------------- !
131 :
132 : integer, intent(in) :: pcols
133 : integer, intent(in) :: pver
134 : integer, intent(in) :: ncnst
135 : integer, intent(in) :: ncol ! Number of atmospheric columns
136 : integer, intent(in) :: lchnk ! Chunk identifier
137 :
138 : real(r8), intent(inout) :: kvm(pcols,pver+1) ! Viscosity ( diffusivity for momentum )
139 : real(r8), intent(out) :: kvt(pcols,pver+1) ! Kinematic molecular conductivity
140 : real(r8), intent(in) :: tint(pcols,pver+1) ! Interface temperature [ K ]
141 : real(r8), intent(in) :: rhoi(pcols,pver+1) ! Density ( rho ) at interfaces
142 : real(r8), intent(in) :: cnst_mw(ncnst) ! Constituent molecular weight
143 :
144 : real(r8), intent(out) :: kq_scal(pcols,pver+1) ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
145 : real(r8), intent(out) :: mw_fac_out(pcols,pver+1,ncnst) ! composition dependent mw_fac on interface level
146 : integer, intent(in) :: nbot_molec
147 :
148 : ! --------------- !
149 : ! Local variables !
150 : ! --------------- !
151 :
152 : integer :: m ! Constituent index
153 : integer :: k ! Level index
154 :
155 0 : real(r8) :: mbarvi(pcols,nbot_molec,ncnst) ! mbarv on interface level
156 :
157 0 : real(r8) :: mkvisc(ncol) ! Molecular kinematic viscosity c*tint**(2/3)/rho
158 :
159 : ! ----------------------- !
160 : ! Main Computation Begins !
161 : ! ----------------------- !
162 :
163 : !
164 : ! Need variable mw_fac for kvt and constant otherwise.
165 : !
166 : ! Then compute molecular kinematic viscosity, heat diffusivity and
167 : ! factor for constituent diffusivity.
168 : ! This is a key part of the code. For WACCM-X, use constituent
169 : ! dependent molecular viscosity and conductivity.
170 : !
171 :
172 0 : kvt = 0._r8
173 0 : kq_scal = 0._r8
174 0 : if (waccmx_mode) then
175 0 : do m = 1, ncnst
176 0 : mbarvi(:ncol,1,m) = .75_r8*mbarv(:ncol,1,lchnk)+0.5_r8*mbarv(:ncol,2,lchnk) &
177 0 : -.25_r8*mbarv(:ncol,3,lchnk)
178 0 : do k = 2, nbot_molec
179 0 : mbarvi(:ncol,k,m) = 0.5_r8 * (mbarv(:ncol,k-1,lchnk)+mbarv(:ncol,k,lchnk))
180 0 : mw_fac_out(:ncol,k,m) = d0 * mbarvi(:ncol,k,m) * &
181 0 : sqrt(1._r8/mbarvi(:ncol,k,m) + 1._r8/cnst_mw(m)) / n_avog
182 : enddo
183 0 : mw_fac_out(:ncol,1,m) = 1.5_r8*mw_fac_out(:ncol,2,m)-.5_r8*mw_fac_out(:ncol,3,m)
184 0 : do k = nbot_molec+1, pver+1
185 0 : mw_fac_out(:ncol,k,m) = mw_fac_out(:ncol,nbot_molec,m)
186 : enddo
187 : end do
188 :
189 0 : do k = 1, nbot_molec
190 0 : mkvisc = kmvis(:ncol,k,lchnk) / rhoi(:ncol,k)
191 0 : kvm(:ncol,k) = kvm(:ncol,k) + mkvisc
192 0 : mkvisc = kmcnd(:ncol,k,lchnk) / rhoi(:ncol,k)
193 0 : kvt(:ncol,k) = mkvisc
194 0 : kq_scal(:ncol,k) = sqrt(tint(:ncol,k)) / rhoi(:ncol,k)
195 : end do
196 :
197 : else
198 0 : do m = 1, ncnst
199 0 : mw_fac_out(:,:,m) = mw_fac(m)
200 : end do
201 :
202 0 : do k = 1, nbot_molec
203 0 : mkvisc = km_fac * tint(:ncol,k)**pwr / rhoi(:ncol,k)
204 0 : kvm(:ncol,k) = kvm(:ncol,k) + mkvisc
205 0 : kvt(:ncol,k) = mkvisc * cpairv(:ncol,k,lchnk)
206 0 : kq_scal(:ncol,k) = sqrt(tint(:ncol,k)) / rhoi(:ncol,k)
207 : end do
208 : endif
209 :
210 0 : end subroutine compute_molec_diff
211 :
212 : !============================================================================ !
213 : ! !
214 : !============================================================================ !
215 :
216 0 : function vd_lu_qdecomp( &
217 : pcols , pver , ncol , fixed_ubc , mw , &
218 0 : kv , kq_scal, mw_facm , dpidz_sq , p , &
219 : interface_boundary, molec_boundary, &
220 0 : tint , ztodt , nbot_molec , &
221 0 : lchnk , t , m , no_molec_decomp) result(decomp)
222 :
223 : use coords_1d, only: Coords1D
224 : use linear_1d_operators, only: BoundaryType, TriDiagDecomp
225 : use vdiff_lu_solver, only: fin_vol_lu_decomp
226 :
227 : !------------------------------------------------------------------------------ !
228 : ! Add the molecular diffusivity to the turbulent diffusivity for a consitutent. !
229 : ! Update the superdiagonal (ca(k)), diagonal (cb(k)) and subdiagonal (cc(k)) !
230 : ! coefficients of the tridiagonal diffusion matrix, also ze and denominator. !
231 : !------------------------------------------------------------------------------ !
232 :
233 : ! ---------------------- !
234 : ! Input-Output Arguments !
235 : ! ---------------------- !
236 :
237 : integer, intent(in) :: pcols
238 : integer, intent(in) :: pver
239 : integer, intent(in) :: ncol ! Number of atmospheric columns
240 :
241 : integer, intent(in) :: nbot_molec
242 :
243 : logical, intent(in) :: fixed_ubc ! Fixed upper boundary condition flag
244 : real(r8), intent(in) :: kv(pcols,pver+1) ! Eddy diffusivity
245 : real(r8), intent(in) :: kq_scal(pcols,pver+1) ! Molecular diffusivity ( kq_fac*sqrt(T)*m_d/rho )
246 : real(r8), intent(in) :: mw ! Molecular weight for this constituent
247 : real(r8), intent(in) :: mw_facm(pcols,pver+1) ! composition dependent sqrt(1/M_q + 1/M_d) for this constituent
248 : real(r8), intent(in) :: dpidz_sq(ncol,pver+1) ! (g*rho)**2 (square of vertical derivative of pint)
249 : type(Coords1D), intent(in) :: p ! Pressure coordinates
250 : type(BoundaryType), intent(in) :: interface_boundary ! Boundary on grid edge.
251 : type(BoundaryType), intent(in) :: molec_boundary ! Boundary at edge of molec_diff region.
252 : real(r8), intent(in) :: tint(pcols,pver+1) ! Interface temperature [ K ]
253 : real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
254 :
255 : integer, intent(in) :: lchnk ! Chunk number
256 : real(r8), intent(in) :: t(pcols,pver) ! temperature
257 : integer, intent(in) :: m ! cnst index
258 :
259 : ! Decomposition covering levels without vertical diffusion.
260 : type(TriDiagDecomp), intent(in) :: no_molec_decomp
261 :
262 : ! LU decomposition information for solver.
263 : type(TriDiagDecomp) :: decomp
264 :
265 : ! --------------- !
266 : ! Local Variables !
267 : ! --------------- !
268 :
269 : ! Level index.
270 : integer :: k
271 :
272 : ! Molecular diffusivity for constituent.
273 0 : real(r8) :: kmq(ncol,nbot_molec+1)
274 :
275 : ! Term for drift due to molecular separation: (m_i/m - 1) / p
276 0 : real(r8) :: mw_term(ncol,nbot_molec+1)
277 :
278 : ! Diffusion coefficient.
279 0 : real(r8) :: diff_coef(ncol,nbot_molec+1)
280 : ! Advection velocity.
281 0 : real(r8) :: advect_v(ncol,nbot_molec+1)
282 :
283 : ! 1/mbar * d(mbar)/dp
284 0 : real(r8) :: gradm(ncol,nbot_molec+1)
285 :
286 : ! alphaTh/T * dT/dp, for now alphaTh is non-zero only for H.
287 0 : real(r8) :: gradt(ncol,nbot_molec+1)
288 :
289 : ! mbarv at interface
290 0 : real(r8) :: mbarvi(ncol)
291 :
292 : ! ----------------------- !
293 : ! Main Computation Begins !
294 : ! ----------------------- !
295 :
296 : ! --------------------------------------------------------------------- !
297 : ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the !
298 : ! tridiagonal diffusion matrix. The diagonal elements (cb=1+ca+cc) are !
299 : ! a combination of ca and cc; they are not required by the solver. !
300 : !---------------------------------------------------------------------- !
301 :
302 0 : call t_startf('vd_lu_qdecomp')
303 :
304 0 : kmq = 0._r8
305 0 : mw_term = 0._r8
306 0 : gradm = 0._r8
307 0 : gradt = 0._r8
308 :
309 : ! Compute difference between scale heights of constituent and dry air
310 :
311 0 : if ( waccmx_mode ) then
312 :
313 : ! Top level first.
314 0 : k = 1
315 0 : mbarvi = .75_r8*mbarv(:ncol,k,lchnk)+0.5_r8*mbarv(:ncol,k+1,lchnk) &
316 0 : -.25_r8*mbarv(:ncol,k+2,lchnk)
317 0 : mw_term(:,k) = (mw/mbarvi - 1._r8) / p%ifc(:,k)
318 0 : gradm(:,k) = (mbarv(:ncol,k,lchnk)-mbarvi)/ &
319 0 : (p%mid(:,k)-p%ifc(:,k))/ &
320 0 : (mbarv(:ncol,k,lchnk)+mbarvi)*2._r8
321 :
322 0 : if (alphath(m) /= 0._r8) then
323 0 : gradt(:,k) = alphath(m)*(t(:ncol,k)-tint(:ncol,k))/ &
324 0 : (p%mid(:ncol,k)-p%ifc(:ncol,k))/ &
325 0 : (t(:ncol,k)+tint(:ncol,k))*2._r8
326 : end if
327 :
328 : ! Interior of molecular diffusion region.
329 0 : do k = 2, nbot_molec
330 0 : mbarvi = 0.5_r8 * (mbarv(:ncol,k-1,lchnk)+mbarv(:ncol,k,lchnk))
331 0 : mw_term(:,k) = (mw/mbarvi - 1._r8) / p%ifc(:,k)
332 0 : gradm(:,k) = (mbarv(:ncol,k,lchnk)-mbarv(:ncol,k-1,lchnk)) * &
333 0 : p%rdst(:,k-1)/mbarvi
334 : enddo
335 :
336 0 : if (alphath(m) /= 0._r8) then
337 0 : do k = 2, nbot_molec
338 0 : gradt(:,k) = alphath(m)*(t(:ncol,k)-t(:ncol,k-1)) &
339 0 : *p%rdst(:,k-1)/tint(:ncol,k)
340 : end do
341 : end if
342 :
343 : ! Leave nbot_molec+1 terms as zero, because molecular diffusion is
344 : ! small at the lower boundary.
345 :
346 : else
347 :
348 0 : do k = 1, nbot_molec
349 0 : mw_term(:,k) = (mw/mw_dry - 1._r8) / p%ifc(:ncol,k)
350 : enddo
351 :
352 : endif
353 :
354 : !-------------------- !
355 : ! Molecular diffusion !
356 : !-------------------- !
357 :
358 : ! Start with non-molecular portion of diffusion.
359 :
360 : ! Molecular diffusion coefficient.
361 0 : do k = 1, nbot_molec
362 0 : kmq(:,k) = kq_scal(:ncol,k) * mw_facm(:ncol,k)
363 : end do
364 :
365 0 : diff_coef = kv(:ncol,:nbot_molec+1) + kmq
366 :
367 : ! "Drift" terms.
368 0 : advect_v = kmq*mw_term
369 0 : if ( waccmx_mode ) then
370 : advect_v = advect_v - kmq*gradt - &
371 0 : (kv(:ncol,:nbot_molec+1) + kmq)*gradm
372 : end if
373 :
374 : ! Convert from z to pressure representation.
375 0 : diff_coef = dpidz_sq(:,:nbot_molec+1) * diff_coef
376 0 : advect_v = dpidz_sq(:,:nbot_molec+1) * advect_v
377 :
378 0 : if( fixed_ubc ) then
379 : decomp = fin_vol_lu_decomp(ztodt, p, &
380 : coef_q_diff=diff_coef, coef_q_adv=advect_v, &
381 : upper_bndry=interface_boundary, &
382 : lower_bndry=molec_boundary, &
383 0 : graft_decomp=no_molec_decomp)
384 : else
385 : decomp = fin_vol_lu_decomp(ztodt, p, &
386 : coef_q_diff=diff_coef, coef_q_adv=advect_v, &
387 : lower_bndry=molec_boundary, &
388 0 : graft_decomp=no_molec_decomp)
389 : end if
390 :
391 0 : call t_stopf('vd_lu_qdecomp')
392 :
393 0 : end function vd_lu_qdecomp
394 :
395 : end module molec_diff
|
