Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cam_molec_diff.F @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 18.5 KB

Line
1	#define WRF_PORT
2
3	module molec_diff
4
5	!------------------------------------------------------------------------------------------------- !
6	! Module to compute molecular diffusivity for various constituents !
7	! !
8	! Public interfaces : !
9	! !
10	! init_molec_diff Initializes time independent coefficients !
11	! init_timestep_molec_diff Time-step initialization for molecular diffusivity !
12	! compute_molec_diff Computes constituent-independent terms for moleculuar diffusivity !
13	! vd_lu_qdecomp Computes constituent-dependent terms for moleculuar diffusivity and !
14	! updates terms in the triadiagonal matrix used for the implicit !
15	! solution of the diffusion equation !
16	! !
17	!---------------------------Code history---------------------------------------------------------- !
18	! Modularized : J. McCaa, September 2004 !
19	! Lastly Arranged : S. Park, January. 2010 !
20	!------------------------------------------------------------------------------------------------- !
21
22	#ifndef WRF_PORT
23	use perf_mod
24	use cam_logfile, only : iulog
25	#else
26	use module_cam_support, only: iulog, t_stopf, t_startf
27	#endif
28
29	implicit none
30	private
31	save
32
33	public init_molec_diff
34	#ifndef WRF_PORT
35	public init_timestep_molec_diff
36	#endif
37	public compute_molec_diff
38	public vd_lu_qdecomp
39
40	! ---------- !
41	! Parameters !
42	! ---------- !
43
44	integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
45
46	real(r8), parameter :: km_fac = 3.55E-7_r8 ! Molecular viscosity constant [ unit ? ]
47	real(r8), parameter :: pr_num = 1._r8 ! Prandtl number [ no unit ]
48	real(r8), parameter :: pwr = 2._r8/3._r8 ! Exponentiation factor [ unit ? ]
49	real(r8), parameter :: d0 = 1.52E20_r8 ! Diffusion factor [ m-1 s-1 ] molec sqrt(kg/kmol/K) [ unit ? ]
50	! Aerononmy, Part B, Banks and Kockarts (1973), p39
51	! Note text cites 1.52E18 cm-1 ...
52
53	real(r8) :: rair ! Gas constant for dry air
54	real(r8) :: mw_dry ! Molecular weight of dry air
55	real(r8) :: n_avog ! Avogadro's number [ molec/kmol ]
56	real(r8) :: gravit
57	real(r8) :: cpair
58	real(r8) :: kbtz ! Boltzman constant
59
60	integer :: ntop_molec ! Top interface level to which molecular vertical diffusion is applied ( = 1 )
61	integer :: nbot_molec ! Bottom interface level to which molecular vertical diffusion is applied ( = pver )
62	real(r8), allocatable :: mw_fac(:) ! sqrt(1/M_q + 1/M_d) in constituent diffusivity [ unit ? ]
63
64	contains
65
66	!============================================================================ !
67	! !
68	!============================================================================ !
69
70	subroutine init_molec_diff( kind, ncnst, rair_in, ntop_molec_in, nbot_molec_in, &
71	mw_dry_in, n_avog_in, gravit_in, cpair_in, kbtz_in )
72
73	use constituents, only : cnst_mw
74	use upper_bc, only : ubc_init
75
76	integer, intent(in) :: kind ! Kind of reals being passed in
77	integer, intent(in) :: ncnst ! Number of constituents
78	integer, intent(in) :: ntop_molec_in ! Top interface level to which molecular vertical diffusion is applied ( = 1 )
79	integer, intent(in) :: nbot_molec_in ! Bottom interface level to which molecular vertical diffusion is applied.
80	real(r8), intent(in) :: rair_in
81	real(r8), intent(in) :: mw_dry_in ! Molecular weight of dry air
82	real(r8), intent(in) :: n_avog_in ! Avogadro's number [ molec/kmol ]
83	real(r8), intent(in) :: gravit_in
84	real(r8), intent(in) :: cpair_in
85	real(r8), intent(in) :: kbtz_in ! Boltzman constant
86	integer :: m ! Constituent index
87
88	if( kind .ne. r8 ) then
89	write(iulog,*) 'KIND of reals passed to init_molec_diff -- exiting.'
90	stop 'init_molec_diff'
91	endif
92
93	rair = rair_in
94	mw_dry = mw_dry_in
95	n_avog = n_avog_in
96	gravit = gravit_in
97	cpair = cpair_in
98	kbtz = kbtz_in
99	ntop_molec = ntop_molec_in
100	nbot_molec = nbot_molec_in
101
102	! Initialize upper boundary condition variables
103
104	call ubc_init()
105
106	! Molecular weight factor in constitutent diffusivity
107	! *** FAKE THIS FOR NOW USING MOLECULAR WEIGHT OF DRY AIR FOR ALL TRACERS **
108
109	allocate(mw_fac(ncnst))
110	do m = 1, ncnst
111	mw_fac(m) = d0 * mw_dry * sqrt(1._r8/mw_dry + 1._r8/cnst_mw(m)) / n_avog
112	end do
113
114	end subroutine init_molec_diff
115
116	!============================================================================ !
117	! !
118	!============================================================================ !
119	#ifndef WRF_PORT
120	subroutine init_timestep_molec_diff(state)
121	!--------------------------- !
122	! Timestep dependent setting !
123	!--------------------------- !
124	use upper_bc, only : ubc_timestep_init
125	use physics_types,only: physics_state
126	use ppgrid, only: begchunk, endchunk
127
128	type(physics_state), intent(in) :: state(begchunk:endchunk)
129
130	call ubc_timestep_init( state)
131
132	end subroutine init_timestep_molec_diff
133	#endif
134	!============================================================================ !
135	! !
136	!============================================================================ !
137
138	integer function compute_molec_diff( lchnk , &
139	pcols , pver , ncnst , ncol , t , pmid , pint , &
140	zi , ztodt , kvh , kvm , tint , rhoi , tmpi2 , &
141	kq_scal , ubc_t , ubc_mmr , dse_top , cc_top , cd_top , cnst_mw_out , &
142	cnst_fixed_ubc_out , mw_fac_out , ntop_molec_out , nbot_molec_out )
143
144	use upper_bc, only : ubc_get_vals
145	use constituents, only : cnst_mw, cnst_fixed_ubc
146
147	! --------------------- !
148	! Input-Output Argument !
149	! --------------------- !
150
151	integer, intent(in) :: pcols
152	integer, intent(in) :: pver
153	integer, intent(in) :: ncnst
154	integer, intent(in) :: ncol ! Number of atmospheric columns
155	integer, intent(in) :: lchnk ! Chunk identifier
156	real(r8), intent(in) :: t(pcols,pver) ! Temperature input
157	real(r8), intent(in) :: pmid(pcols,pver) ! Midpoint pressures
158	real(r8), intent(in) :: pint(pcols,pver+1) ! Interface pressures
159	real(r8), intent(in) :: zi(pcols,pver+1) ! Interface heights
160	real(r8), intent(in) :: ztodt ! 2 delta-t
161
162	real(r8), intent(inout) :: kvh(pcols,pver+1) ! Diffusivity for heat
163	real(r8), intent(inout) :: kvm(pcols,pver+1) ! Viscosity ( diffusivity for momentum )
164	real(r8), intent(inout) :: tint(pcols,pver+1) ! Interface temperature
165	real(r8), intent(inout) :: rhoi(pcols,pver+1) ! Density ( rho ) at interfaces
166	real(r8), intent(inout) :: tmpi2(pcols,pver+1) ! dt(grho)**2/dp at interfaces
167
168	real(r8), intent(out) :: kq_scal(pcols,pver+1) ! kq_facsqrt(T)m_d/rho for molecular diffusivity
169	real(r8), intent(out) :: ubc_mmr(pcols,ncnst) ! Upper boundary mixing ratios [ kg/kg ]
170	real(r8), intent(out) :: cnst_mw_out(ncnst)
171	logical, intent(out) :: cnst_fixed_ubc_out(ncnst)
172	real(r8), intent(out) :: mw_fac_out(ncnst)
173	real(r8), intent(out) :: dse_top(pcols) ! dse on top boundary
174	real(r8), intent(out) :: cc_top(pcols) ! Lower diagonal at top interface
175	real(r8), intent(out) :: cd_top(pcols) ! cc_top * dse ubc value
176	integer, intent(out) :: ntop_molec_out
177	integer, intent(out) :: nbot_molec_out
178
179	! --------------- !
180	! Local variables !
181	! --------------- !
182
183	integer :: m ! Constituent index
184	integer :: i ! Column index
185	integer :: k ! Level index
186	real(r8) :: mkvisc ! Molecular kinematic viscosity ctint*(2/3)/rho
187	real(r8) :: ubc_t(pcols) ! Upper boundary temperature (K)
188
189	! ----------------------- !
190	! Main Computation Begins !
191	! ----------------------- !
192
193	! Get upper boundary values
194
195	call ubc_get_vals( lchnk, ncol, ntop_molec, pint, zi, ubc_t, ubc_mmr )
196
197	! Below are already computed, just need to be copied for output
198
199	cnst_mw_out(:ncnst) = cnst_mw(:ncnst)
200	cnst_fixed_ubc_out(:ncnst) = cnst_fixed_ubc(:ncnst)
201	mw_fac_out(:ncnst) = mw_fac(:ncnst)
202	ntop_molec_out = ntop_molec
203	nbot_molec_out = nbot_molec
204
205	! Density and related factors for moecular diffusion and ubc.
206	! Always have a fixed upper boundary T if molecular diffusion is active. Why ?
207
208	tint (:ncol,ntop_molec) = ubc_t(:ncol)
209	rhoi (:ncol,ntop_molec) = pint(:ncol,ntop_molec) / ( rair * tint(:ncol,ntop_molec) )
210	tmpi2(:ncol,ntop_molec) = ztodt * ( gravit * rhoi(:ncol,ntop_molec))**2 &
211	/ ( pmid(:ncol,ntop_molec) - pint(:ncol,ntop_molec) )
212
213	! Compute molecular kinematic viscosity, heat diffusivity and factor for constituent diffusivity
214	! This is a key part of the code.
215
216	kq_scal(:ncol,1:ntop_molec-1) = 0._r8
217	do k = ntop_molec, nbot_molec
218	do i = 1, ncol
219	mkvisc = km_fac * tint(i,k)**pwr / rhoi(i,k)
220	kvm(i,k) = kvm(i,k) + mkvisc
221	kvh(i,k) = kvh(i,k) + mkvisc * pr_num
222	kq_scal(i,k) = sqrt(tint(i,k)) / rhoi(i,k)
223	end do
224	end do
225	kq_scal(:ncol,nbot_molec+1:pver+1) = 0._r8
226
227	! Top boundary condition for dry static energy
228
229	dse_top(:ncol) = cpair * tint(:ncol,ntop_molec) + gravit * zi(:ncol,ntop_molec)
230
231	! Top value of cc for dry static energy
232
233	do i = 1, ncol
234	cc_top(i) = ztodt * gravit*2 rhoi(i,ntop_molec) * km_fac * ubc_t(i)**pwr / &
235	( ( pint(i,2) - pint(i,1) ) * ( pmid(i,1) - pint(i,1) ) )
236	enddo
237	cd_top(:ncol) = cc_top(:ncol) * dse_top(:ncol)
238
239	compute_molec_diff = 1
240	return
241	end function compute_molec_diff
242
243	!============================================================================ !
244	! !
245	!============================================================================ !
246
247	integer function vd_lu_qdecomp( pcols , pver , ncol , fixed_ubc , mw , ubc_mmr , &
248	kv , kq_scal, mw_facm , tmpi , rpdel , &
249	ca , cc , dnom , ze , rhoi , &
250	tint , ztodt , ntop_molec , nbot_molec , cd_top )
251
252	!------------------------------------------------------------------------------ !
253	! Add the molecular diffusivity to the turbulent diffusivity for a consitutent. !
254	! Update the superdiagonal (ca(k)), diagonal (cb(k)) and subdiagonal (cc(k)) !
255	! coefficients of the tridiagonal diffusion matrix, also ze and denominator. !
256	!------------------------------------------------------------------------------ !
257
258	! ---------------------- !
259	! Input-Output Arguments !
260	! ---------------------- !
261
262	integer, intent(in) :: pcols
263	integer, intent(in) :: pver
264	integer, intent(in) :: ncol ! Number of atmospheric columns
265
266	integer, intent(in) :: ntop_molec
267	integer, intent(in) :: nbot_molec
268
269	logical, intent(in) :: fixed_ubc ! Fixed upper boundary condition flag
270	real(r8), intent(in) :: kv(pcols,pver+1) ! Eddy diffusivity
271	real(r8), intent(in) :: kq_scal(pcols,pver+1) ! Molecular diffusivity ( kq_facsqrt(T)m_d/rho )
272	real(r8), intent(in) :: mw ! Molecular weight for this constituent
273	real(r8), intent(in) :: ubc_mmr(pcols) ! Upper boundary mixing ratios [ kg/kg ]
274	real(r8), intent(in) :: mw_facm ! sqrt(1/M_q + 1/M_d) for this constituent
275	real(r8), intent(in) :: tmpi(pcols,pver+1) ! dt(g/R)2/dppi(k+1)/(.5(tm(k+1)+tm(k))*2
276	real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel ( thickness bet interfaces )
277	real(r8), intent(in) :: rhoi(pcols,pver+1) ! Density at interfaces [ kg/m3 ]
278	real(r8), intent(in) :: tint(pcols,pver+1) ! Interface temperature [ K ]
279	real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
280
281	real(r8), intent(inout) :: ca(pcols,pver) ! -Upper diagonal
282	real(r8), intent(inout) :: cc(pcols,pver) ! -Lower diagonal
283	real(r8), intent(inout) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)ze(k-1)) , 1./(b(k) - c(k)e(k-1))
284	real(r8), intent(inout) :: ze(pcols,pver) ! Term in tri-diag. matrix system
285
286	real(r8), intent(out) :: cd_top(pcols) ! Term for updating top level with ubc
287
288	! --------------- !
289	! Local Variables !
290	! --------------- !
291
292	integer :: i ! Longitude index
293	integer :: k, kp1 ! Vertical indicies
294
295	real(r8) :: rghd(pcols,pver+1) ! (1/H_i - 1/H) * (rho*g)^(-1)
296	real(r8) :: kmq(ncol) ! Molecular diffusivity for constituent
297	real(r8) :: wrk0(ncol) ! Work variable
298	real(r8) :: wrk1(ncol) ! Work variable
299
300	real(r8) :: cb(pcols,pver) ! - Diagonal
301	real(r8) :: kvq(pcols,pver+1) ! Output vertical diffusion coefficient
302
303	! ----------------------- !
304	! Main Computation Begins !
305	! ----------------------- !
306
307	! --------------------------------------------------------------------- !
308	! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the !
309	! tridiagonal diffusion matrix. The diagonal elements (cb=1+ca+cc) are !
310	! a combination of ca and cc; they are not required by the solver. !
311	!---------------------------------------------------------------------- !
312
313	call t_startf('vd_lu_qdecomp')
314
315	kvq(:,:) = 0._r8
316	cd_top(:) = 0._r8
317
318	! Compute difference between scale heights of constituent and dry air
319
320	do k = ntop_molec, nbot_molec
321	do i = 1, ncol
322	rghd(i,k) = gravit / ( kbtz * n_avog * tint(i,k) ) * ( mw - mw_dry )
323	rghd(i,k) = ztodt * gravit * rhoi(i,k) * rghd(i,k)
324	enddo
325	enddo
326
327	!-------------------- !
328	! Molecular diffusion !
329	!-------------------- !
330
331	do k = nbot_molec - 1, ntop_molec, -1
332	kp1 = k + 1
333	kmq(:ncol) = kq_scal(:ncol,kp1) * mw_facm
334	wrk0(:ncol) = ( kv(:ncol,kp1) + kmq(:ncol) ) * tmpi(:ncol,kp1)
335	wrk1(:ncol) = kmq(:ncol) * 0.5_r8 * rghd(:ncol,kp1)
336	! Add species separation term
337	ca(:ncol,k ) = ( wrk0(:ncol) - wrk1(:ncol) ) * rpdel(:ncol,k)
338	cc(:ncol,kp1) = ( wrk0(:ncol) + wrk1(:ncol) ) * rpdel(:ncol,kp1)
339	kvq(:ncol,kp1) = kmq(:ncol)
340	end do
341
342	if( fixed_ubc ) then
343	cc(:ncol,ntop_molec) = kq_scal(:ncol,ntop_molec) * mw_facm &
344	* ( tmpi(:ncol,ntop_molec) + rghd(:ncol,ntop_molec) ) &
345	* rpdel(:ncol,ntop_molec)
346	end if
347
348	! Calculate diagonal elements
349
350	do k = nbot_molec - 1, ntop_molec + 1, -1
351	kp1 = k + 1
352	cb(:ncol,k) = 1._r8 + ca(:ncol,k) + cc(:ncol,k) &
353	+ rpdel(:ncol,k) * ( kvq(:ncol,kp1) * rghd(:ncol,kp1) &
354	- kvq(:ncol,k) * rghd(:ncol,k) )
355	kvq(:ncol,kp1) = kv(:ncol,kp1) + kvq(:ncol,kp1)
356	end do
357
358	k = ntop_molec
359	kp1 = k + 1
360	if( fixed_ubc ) then
361	cb(:ncol,k) = 1._r8 + ca(:ncol,k) &
362	+ rpdel(:ncol,k) * kvq(:ncol,kp1) * rghd(:ncol,kp1) &
363	+ kq_scal(:ncol,ntop_molec) * mw_facm &
364	* ( tmpi(:ncol,ntop_molec) - rghd(:ncol,ntop_molec) ) &
365	* rpdel(:ncol,ntop_molec)
366	else
367	cb(:ncol,k) = 1._r8 + ca(:ncol,k) &
368	+ rpdel(:ncol,k) * kvq(:ncol,kp1) * rghd(:ncol,kp1)
369	end if
370
371	k = nbot_molec
372	cb(:ncol,k) = 1._r8 + cc(:ncol,k) + ca(:ncol,k) &
373	- rpdel(:ncol,k) * kvq(:ncol,k)*rghd(:ncol,k)
374	do k = 1, nbot_molec + 1, -1
375	cb(:ncol,k) = 1._r8 + ca(:ncol,k) + cc(:ncol,k)
376	end do
377
378	! Compute term for updating top level mixing ratio for ubc
379
380	if( fixed_ubc ) then
381	cd_top(:ncol) = cc(:ncol,ntop_molec) * ubc_mmr(:ncol)
382	end if
383
384	!-------------------------------------------------------- !
385	! Calculate e(k). !
386	! This term is required in solution of tridiagonal matrix !
387	! defined by implicit diffusion equation. !
388	!-------------------------------------------------------- !
389
390	do k = nbot_molec, ntop_molec + 1, -1
391	dnom(:ncol,k) = 1._r8 / ( cb(:ncol,k) - ca(:ncol,k) * ze(:ncol,k+1) )
392	ze(:ncol,k) = cc(:ncol,k) * dnom(:ncol,k)
393	end do
394	k = ntop_molec
395	dnom(:ncol,k) = 1._r8 / ( cb(:ncol,k) - ca(:ncol,k) * ze(:ncol,k+1) )
396
397	vd_lu_qdecomp = 1
398	call t_stopf('vd_lu_qdecomp')
399	return
400
401	end function vd_lu_qdecomp
402
403	end module molec_diff

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format