Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cam_bl_diffusion_solver.F @ 1

Last change on this file since 1 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: 47.0 KB

Line
1	#define WRF_PORT
2
3	module diffusion_solver
4
5	!------------------------------------------------------------------------------------ !
6	! Module to solve vertical diffusion equations using a tri-diagonal solver. !
7	! The module will also apply countergradient fluxes, and apply molecular !
8	! diffusion for constituents. !
9	! !
10	! Public interfaces : !
11	! init_vdiff initializes time independent coefficients !
12	! compute_vdiff solves diffusion equations !
13	! vd_lu_solve tridiagonal solver also used by gwd (should be private) !
14	! vd_lu_decomp tridiagonal solver also used by gwd (should be private) !
15	! vdiff_selector type for storing fields selected to be diffused !
16	! vdiff_select selects fields to be diffused !
17	! operator(.not.) extends .not. to operate on type vdiff_selector !
18	! any provides functionality of intrinsic any for type vdiff_selector !
19	! !
20	!------------------------------------ Code History ---------------------------------- !
21	! Initial subroutines : B. Boville and others, 1991-2004 !
22	! Modularization : J. McCaa, September 2004 !
23	! Most Recent Code : Sungsu Park, Aug. 2006, Dec. 2008, Jan. 2010. !
24	!------------------------------------------------------------------------------------ !
25
26
27	#ifndef WRF_PORT
28	use cam_logfile, only : iulog
29	#else
30	use module_cam_support, only: iulog
31	#endif
32
33	implicit none
34	private
35	save
36
37	integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
38
39	! ----------------- !
40	! Public interfaces !
41	! ----------------- !
42
43	public init_vdiff ! Initialization
44	public compute_vdiff ! Full routine
45	public vd_lu_solve ! Tridiagonal solver also used by gwd ( should be private! )
46	public vd_lu_decomp ! Tridiagonal solver also used by gwd ( should be private! )
47	public vdiff_selector ! Type for storing fields selected to be diffused
48	public vdiff_select ! Selects fields to be diffused
49	public operator(.not.) ! Extends .not. to operate on type vdiff_selector
50	public any ! Provides functionality of intrinsic any for type vdiff_selector
51
52	integer, public :: nbot_molec ! Bottom level where molecular diffusivity is applied
53
54	! Below stores logical array of fields to be diffused
55
56	type vdiff_selector
57	private
58	logical, pointer, dimension(:) :: fields
59	end type vdiff_selector
60
61	! Below extends .not. to operate on type vdiff_selector
62
63	interface operator(.not.)
64	module procedure not
65	end interface
66
67	! Below provides functionality of intrinsic any for type vdiff_selector
68
69	interface any
70	module procedure my_any
71	end interface
72
73	! ------------ !
74	! Private data !
75	! ------------ !
76
77	real(r8), private :: cpair ! Specific heat of dry air
78	real(r8), private :: gravit ! Acceleration due to gravity
79	real(r8), private :: rair ! Gas constant for dry air
80	real(r8), private :: zvir ! rh2o/rair - 1
81	real(r8), private :: latvap ! Latent heat of vaporization
82	real(r8), private :: karman ! von Karman constant
83
84	! Parameters used for Turbulent Mountain Stress
85
86	real(r8), parameter :: z0fac = 0.025_r8 ! Factor determining z_0 from orographic standard deviation
87	real(r8), parameter :: z0max = 100._r8 ! Max value of z_0 for orography
88	real(r8), parameter :: horomin = 10._r8 ! Min value of subgrid orographic height for mountain stress
89	real(r8), parameter :: dv2min = 0.01_r8 ! Minimum shear squared
90	real(r8), private :: oroconst ! Converts from standard deviation to height
91
92	contains
93
94	! =============================================================================== !
95	! !
96	! =============================================================================== !
97
98	subroutine init_vdiff( kind, ncnst, rair_in, gravit_in, fieldlist_wet, fieldlist_dry, errstring )
99
100	integer, intent(in) :: kind ! Kind used for reals
101	integer, intent(in) :: ncnst ! Number of constituents
102	real(r8), intent(in) :: rair_in ! Input gas constant for dry air
103	real(r8), intent(in) :: gravit_in ! Input gravititational acceleration
104	type(vdiff_selector), intent(out) :: fieldlist_wet ! List of fields to be diffused using moist mixing ratio
105	type(vdiff_selector), intent(out) :: fieldlist_dry ! List of fields to be diffused using dry mixing ratio
106	character(128), intent(out) :: errstring ! Output status
107
108	errstring = ''
109	if( kind .ne. r8 ) then
110	write(iulog,*) 'KIND of reals passed to init_vdiff -- exiting.'
111	errstring = 'init_vdiff'
112	return
113	endif
114
115	rair = rair_in
116	gravit = gravit_in
117
118	allocate( fieldlist_wet%fields( 3 + ncnst ) )
119	fieldlist_wet%fields(:) = .false.
120
121	allocate( fieldlist_dry%fields( 3 + ncnst ) )
122	fieldlist_dry%fields(:) = .false.
123
124	end subroutine init_vdiff
125
126	! =============================================================================== !
127	! !
128	! =============================================================================== !
129
130	subroutine compute_vdiff( lchnk , &
131	pcols , pver , ncnst , ncol , pmid , &
132	pint , rpdel , t , ztodt , taux , &
133	tauy , shflx , cflx , ntop , nbot , &
134	kvh , kvm , kvq , cgs , cgh , &
135	zi , ksrftms , qmincg , fieldlist , &
136	u , v , q , dse , &
137	tautmsx , tautmsy , dtk , topflx , errstring , &
138	tauresx , tauresy , itaures , &
139	do_molec_diff , compute_molec_diff , vd_lu_qdecomp )
140
141	!-------------------------------------------------------------------------- !
142	! Driver routine to compute vertical diffusion of momentum, moisture, trace !
143	! constituents and dry static energy. The new temperature is computed from !
144	! the diffused dry static energy. !
145	! Turbulent diffusivities and boundary layer nonlocal transport terms are !
146	! obtained from the turbulence module. !
147	!-------------------------------------------------------------------------- !
148	#ifndef WRF_PORT
149	use phys_debug_util, only : phys_debug_col
150	use time_manager, only : is_first_step, get_nstep
151	use phys_control, only : phys_getopts
152	#endif
153
154	! Modification : Ideally, we should diffuse 'liquid-ice static energy' (sl), not the dry static energy.
155	! Also, vertical diffusion of cloud droplet number concentration and aerosol number
156	! concentration should be done very carefully in the future version.
157
158	! --------------- !
159	! Input Arguments !
160	! --------------- !
161
162	integer, intent(in) :: lchnk
163	integer, intent(in) :: pcols
164	integer, intent(in) :: pver
165	integer, intent(in) :: ncnst
166	integer, intent(in) :: ncol ! Number of atmospheric columns
167	integer, intent(in) :: ntop ! Top interface level to which vertical diffusion is applied ( = 1 ).
168	integer, intent(in) :: nbot ! Bottom interface level to which vertical diffusion is applied ( = pver ).
169	integer, intent(in) :: itaures ! Indicator determining whether 'tauresx,tauresy' is updated (1) or non-updated (0) in this subroutine.
170
171	real(r8), intent(in) :: pmid(pcols,pver) ! Mid-point pressures [ Pa ]
172	real(r8), intent(in) :: pint(pcols,pver+1) ! Interface pressures [ Pa ]
173	real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel
174	real(r8), intent(in) :: t(pcols,pver) ! Temperature [ K ]
175	real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
176	real(r8), intent(in) :: taux(pcols) ! Surface zonal stress. Input u-momentum per unit time per unit area into the atmosphere [ N/m2 ]
177	real(r8), intent(in) :: tauy(pcols) ! Surface meridional stress. Input v-momentum per unit time per unit area into the atmosphere [ N/m2 ]
178	real(r8), intent(in) :: shflx(pcols) ! Surface sensible heat flux [ W/m2 ]
179	real(r8), intent(in) :: cflx(pcols,ncnst) ! Surface constituent flux [ kg/m2/s ]
180	real(r8), intent(in) :: zi(pcols,pver+1) ! Interface heights [ m ]
181	real(r8), intent(in) :: ksrftms(pcols) ! Surface drag coefficient for turbulent mountain stress. > 0. [ kg/s/m2 ]
182	real(r8), intent(in) :: qmincg(ncnst) ! Minimum constituent mixing ratios from cg fluxes
183
184	logical, intent(in) :: do_molec_diff ! Flag indicating multiple constituent diffusivities
185	integer, external, optional :: compute_molec_diff ! Constituent-independent moleculuar diffusivity routine
186	integer, external, optional :: vd_lu_qdecomp ! Constituent-dependent moleculuar diffusivity routine
187	type(vdiff_selector), intent(in) :: fieldlist ! Array of flags selecting which fields to diffuse
188
189	! ---------------------- !
190	! Input-Output Arguments !
191	! ---------------------- !
192
193	real(r8), intent(inout) :: kvh(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
194	real(r8), intent(inout) :: kvm(pcols,pver+1) ! Eddy viscosity ( Eddy diffusivity for momentum ) [ m2/s ]
195	real(r8), intent(inout) :: kvq(pcols,pver+1) ! Eddy diffusivity for constituents
196	real(r8), intent(inout) :: cgs(pcols,pver+1) ! Counter-gradient star [ cg/flux ]
197	real(r8), intent(inout) :: cgh(pcols,pver+1) ! Counter-gradient term for heat
198
199	real(r8), intent(inout) :: u(pcols,pver) ! U wind. This input is the 'raw' input wind to PBL scheme without iterative provisional update. [ m/s ]
200	real(r8), intent(inout) :: v(pcols,pver) ! V wind. This input is the 'raw' input wind to PBL scheme without iterative provisional update. [ m/s ]
201	real(r8), intent(inout) :: q(pcols,pver,ncnst) ! Moisture and trace constituents [ kg/kg, #/kg ? ]
202	real(r8), intent(inout) :: dse(pcols,pver) ! Dry static energy [ J/kg ]
203
204	real(r8), intent(inout) :: tauresx(pcols) ! Input : Reserved surface stress at previous time step
205	real(r8), intent(inout) :: tauresy(pcols) ! Output : Reserved surface stress at current time step
206
207	! ---------------- !
208	! Output Arguments !
209	! ---------------- !
210
211	real(r8), intent(out) :: dtk(pcols,pver) ! T tendency from KE dissipation
212	real(r8), intent(out) :: tautmsx(pcols) ! Implicit zonal turbulent mountain surface stress [ N/m2 = kg m/s /s/m2 ]
213	real(r8), intent(out) :: tautmsy(pcols) ! Implicit meridional turbulent mountain surface stress [ N/m2 = kg m/s /s/m2 ]
214	real(r8), intent(out) :: topflx(pcols) ! Molecular heat flux at the top interface
215	character(128), intent(out) :: errstring ! Output status
216
217	! --------------- !
218	! Local Variables !
219	! --------------- !
220
221	integer :: i, k, m, icol ! Longitude, level, constituent indices
222	integer :: status ! Status indicator
223	integer :: ntop_molec ! Top level where molecular diffusivity is applied
224	logical :: lqtst(pcols) ! Adjust vertical profiles
225	logical :: need_decomp ! Whether to compute a new decomposition
226	logical :: cnst_fixed_ubc(ncnst) ! Whether upper boundary condition is fixed
227	logical :: do_iss ! Use implicit turbulent surface stress computation
228
229	real(r8) :: tmpm(pcols,pver) ! Potential temperature, ze term in tri-diag sol'n
230	real(r8) :: ca(pcols,pver) ! - Upper diag of tri-diag matrix
231	real(r8) :: cc(pcols,pver) ! - Lower diag of tri-diag matrix
232	real(r8) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
233
234	real(r8) :: tmp1(pcols) ! Temporary storage
235	real(r8) :: tmpi1(pcols,pver+1) ! Interface KE dissipation
236	real(r8) :: tint(pcols,pver+1) ! Interface temperature
237	real(r8) :: rhoi(pcols,pver+1) ! rho at interfaces
238	real(r8) :: tmpi2(pcols,pver+1) ! dt(grho)**2/dp at interfaces
239	real(r8) :: rrho(pcols) ! 1./bottom level density
240
241	real(r8) :: zero(pcols) ! Zero array for surface heat exchange coefficients
242	real(r8) :: tautotx(pcols) ! Total surface stress ( zonal )
243	real(r8) :: tautoty(pcols) ! Total surface stress ( meridional )
244
245	real(r8) :: dinp_u(pcols,pver+1) ! Vertical difference at interfaces, input u
246	real(r8) :: dinp_v(pcols,pver+1) ! Vertical difference at interfaces, input v
247	real(r8) :: dout_u ! Vertical difference at interfaces, output u
248	real(r8) :: dout_v ! Vertical difference at interfaces, output v
249	real(r8) :: dse_top(pcols) ! dse on top boundary
250	real(r8) :: cc_top(pcols) ! Lower diagonal at top interface
251	real(r8) :: cd_top(pcols) !
252	real(r8) :: rghd(pcols,pver+1) ! (1/H_i - 1/H) (grho)^(-1)
253
254	real(r8) :: qtm(pcols,pver) ! Temporary copy of q
255	real(r8) :: kq_scal(pcols,pver+1) ! kq_facsqrt(T)m_d/rho for molecular diffusivity
256	real(r8) :: mw_fac(ncnst) ! sqrt(1/M_q + 1/M_d) for this constituent
257	real(r8) :: cnst_mw(ncnst) ! Molecular weight [ kg/kmole ]
258	real(r8) :: ubc_mmr(pcols,ncnst) ! Upper boundary mixing ratios [ kg/kg ]
259	real(r8) :: ubc_t(pcols) ! Upper boundary temperature [ K ]
260
261	real(r8) :: ws(pcols) ! Lowest-level wind speed [ m/s ]
262	real(r8) :: tau(pcols) ! Turbulent surface stress ( not including mountain stress )
263	real(r8) :: ksrfturb(pcols) ! Surface drag coefficient of 'normal' stress. > 0. Virtual mass input per unit time per unit area [ kg/s/m2 ]
264	real(r8) :: ksrf(pcols) ! Surface drag coefficient of 'normal' stress + Surface drag coefficient of 'tms' stress. > 0. [ kg/s/m2 ]
265	real(r8) :: usum_in(pcols) ! Vertical integral of input u-momentum. Total zonal momentum per unit area in column [ sum of u*dp/g = kg m/s m-2 ]
266	real(r8) :: vsum_in(pcols) ! Vertical integral of input v-momentum. Total meridional momentum per unit area in column [ sum of v*dp/g = kg m/s m-2 ]
267	real(r8) :: usum_mid(pcols) ! Vertical integral of u-momentum after adding explicit residual stress
268	real(r8) :: vsum_mid(pcols) ! Vertical integral of v-momentum after adding explicit residual stress
269	real(r8) :: usum_out(pcols) ! Vertical integral of u-momentum after doing implicit diffusion
270	real(r8) :: vsum_out(pcols) ! Vertical integral of v-momentum after doing implicit diffusion
271	real(r8) :: tauimpx(pcols) ! Actual net stress added at the current step other than mountain stress
272	real(r8) :: tauimpy(pcols) ! Actual net stress added at the current step other than mountain stress
273	real(r8) :: wsmin ! Minimum sfc wind speed for estimating frictional transfer velocity ksrf. [ m/s ]
274	real(r8) :: ksrfmin ! Minimum surface drag coefficient [ kg/s/m^2 ]
275	real(r8) :: timeres ! Relaxation time scale of residual stress ( >= dt ) [ s ]
276	real(r8) :: ramda ! dt/timeres [ no unit ]
277	real(r8) :: psum
278	real(r8) :: u_in, u_res, tauresx_in
279	real(r8) :: v_in, v_res, tauresy_in
280
281	! ------------------------------------------------ !
282	! Parameters for implicit surface stress treatment !
283	! ------------------------------------------------ !
284
285	wsmin = 1._r8 ! Minimum wind speed for ksrfturb computation [ m/s ]
286	ksrfmin = 1.e-4_r8 ! Minimum surface drag coefficient [ kg/s/m^2 ]
287	timeres = 7200._r8 ! Relaxation time scale of residual stress ( >= dt ) [ s ]
288	#ifndef WRF_PORT
289	call phys_getopts( do_iss_out = do_iss )
290	#endif
291	do_iss = .true.
292
293	! ----------------------- !
294	! Main Computation Begins !
295	! ----------------------- !
296
297	errstring = ''
298	if( ( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) .and. .not. diffuse(fieldlist,'s') ) then
299	errstring = 'diffusion_solver.compute_vdiff: must diffuse s if diffusing u or v'
300	return
301	end if
302	zero(:) = 0._r8
303
304	! Compute 'rho' and 'dt(grho)^2/dp' at interfaces
305
306	tint(:ncol,1) = t(:ncol,1)
307	rhoi(:ncol,1) = pint(:ncol,1) / (rair*tint(:ncol,1))
308	do k = 2, pver
309	do i = 1, ncol
310	tint(i,k) = 0.5_r8 * ( t(i,k) + t(i,k-1) )
311	rhoi(i,k) = pint(i,k) / (rair*tint(i,k))
312	tmpi2(i,k) = ztodt * ( gravitrhoi(i,k) )*2 / ( pmid(i,k) - pmid(i,k-1) )
313	end do
314	end do
315	tint(:ncol,pver+1) = t(:ncol,pver)
316	rhoi(:ncol,pver+1) = pint(:ncol,pver+1) / ( rair*tint(:ncol,pver+1) )
317
318	rrho(:ncol) = rair * t(:ncol,pver) / pmid(:ncol,pver)
319	tmp1(:ncol) = ztodt * gravit * rpdel(:ncol,pver)
320
321	!--------------------------------------- !
322	! Computation of Molecular Diffusivities !
323	!--------------------------------------- !
324
325	! Modification : Why 'kvq' is not changed by molecular diffusion ?
326
327	if( do_molec_diff ) then
328
329	if( (.not.present(compute_molec_diff)) .or. (.not.present(vd_lu_qdecomp)) ) then
330	errstring = 'compute_vdiff: do_molec_diff true but compute_molec_diff or vd_lu_qdecomp missing'
331	return
332	endif
333
334	! The next subroutine 'compute_molec_diff' :
335	! Modifies : kvh, kvm, tint, rhoi, and tmpi2
336	! Returns : kq_scal, ubc_t, ubc_mmr, dse_top, cc_top, cd_top, cnst_mw,
337	! cnst_fixed_ubc , mw_fac , ntop_molec
338
339	status = compute_molec_diff( lchnk , &
340	pcols , pver , ncnst , ncol , t , pmid , pint , &
341	zi , ztodt , kvh , kvm , tint , rhoi , tmpi2 , &
342	kq_scal , ubc_t , ubc_mmr , dse_top , cc_top , cd_top , cnst_mw , &
343	cnst_fixed_ubc , mw_fac , ntop_molec , nbot_molec )
344
345	else
346
347	kq_scal(:,:) = 0._r8
348	cd_top(:) = 0._r8
349	cc_top(:) = 0._r8
350
351	endif
352
353	!---------------------------- !
354	! Diffuse Horizontal Momentum !
355	!---------------------------- !
356
357	if( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) then
358
359	! Compute the vertical upward differences of the input u,v for KE dissipation
360	! at each interface.
361	! Velocity = 0 at surface, so difference at the bottom interface is -u,v(pver)
362	! These 'dinp_u, dinp_v' are computed using the non-diffused input wind.
363
364	do i = 1, ncol
365	dinp_u(i,1) = 0._r8
366	dinp_v(i,1) = 0._r8
367	dinp_u(i,pver+1) = -u(i,pver)
368	dinp_v(i,pver+1) = -v(i,pver)
369	end do
370	do k = 2, pver
371	do i = 1, ncol
372	dinp_u(i,k) = u(i,k) - u(i,k-1)
373	dinp_v(i,k) = v(i,k) - v(i,k-1)
374	end do
375	end do
376
377	! -------------------------------------------------------------- !
378	! Do 'Implicit Surface Stress' treatment for numerical stability !
379	! in the lowest model layer. !
380	! -------------------------------------------------------------- !
381
382	if( do_iss ) then
383
384	! Compute surface drag coefficient for implicit diffusion
385	! including turbulent turbulent mountain stress.
386
387	do i = 1, ncol
388	ws(i) = max( sqrt( u(i,pver)2._r8 + v(i,pver)2._r8 ), wsmin )
389	tau(i) = sqrt( taux(i)2._r8 + tauy(i)2._r8 )
390	ksrfturb(i) = max( tau(i) / ws(i), ksrfmin )
391	end do
392	ksrf(:ncol) = ksrfturb(:ncol) + ksrftms(:ncol) ! Do all surface stress ( normal + tms ) implicitly
393
394	! Vertical integration of input momentum.
395	! This is total horizontal momentum per unit area [ kg*m/s/m2 ] in each column.
396	! Note (u,v) are the raw input to the PBL scheme, not the
397	! provisionally-marched ones within the iteration loop of the PBL scheme.
398
399	do i = 1, ncol
400	usum_in(i) = 0._r8
401	vsum_in(i) = 0._r8
402	do k = 1, pver
403	usum_in(i) = usum_in(i) + (1._r8/gravit)*u(i,k)/rpdel(i,k)
404	vsum_in(i) = vsum_in(i) + (1._r8/gravit)*v(i,k)/rpdel(i,k)
405	end do
406	end do
407
408	! Add residual stress of previous time step explicitly into the lowest
409	! model layer with a relaxation time scale of 'timeres'.
410
411	ramda = ztodt / timeres
412	u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)tauresx(:ncol)ramda
413	v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)tauresy(:ncol)ramda
414
415	! Vertical integration of momentum after adding explicit residual stress
416	! into the lowest model layer.
417
418	do i = 1, ncol
419	usum_mid(i) = 0._r8
420	vsum_mid(i) = 0._r8
421	do k = 1, pver
422	usum_mid(i) = usum_mid(i) + (1._r8/gravit)*u(i,k)/rpdel(i,k)
423	vsum_mid(i) = vsum_mid(i) + (1._r8/gravit)*v(i,k)/rpdel(i,k)
424	end do
425	end do
426
427	! Debug
428	! icol = phys_debug_col(lchnk)
429	! if ( icol > 0 .and. get_nstep() .ge. 1 ) then
430	! tauresx_in = tauresx(icol)
431	! tauresy_in = tauresy(icol)
432	! u_in = u(icol,pver) - tmp1(icol) * tauresx(icol) * ramda
433	! v_in = v(icol,pver) - tmp1(icol) * tauresy(icol) * ramda
434	! u_res = u(icol,pver)
435	! v_res = v(icol,pver)
436	! endif
437	! Debug
438
439	else
440
441	! In this case, do 'turbulent mountain stress' implicitly,
442	! but do 'normal turbulent stress' explicitly.
443	! In this case, there is no 'redisual stress' as long as 'tms' is
444	! treated in a fully implicit wway, which is true.
445
446	! 1. Do 'tms' implicitly
447
448	ksrf(:ncol) = ksrftms(:ncol)
449
450	! 2. Do 'normal stress' explicitly
451
452	u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*taux(:ncol)
453	v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauy(:ncol)
454
455	end if ! End of 'do iss' ( implicit surface stress )
456
457	! --------------------------------------------------------------------------------------- !
458	! Diffuse horizontal momentum implicitly using tri-diagnonal matrix. !
459	! The 'u,v' are input-output: the output 'u,v' are implicitly diffused winds. !
460	! For implicit 'normal' stress : ksrf = ksrftms + ksrfturb, !
461	! u(pver) : explicitly include 'redisual normal' stress !
462	! For explicit 'normal' stress : ksrf = ksrftms !
463	! u(pver) : explicitly include 'normal' stress !
464	! Note that in all the two cases above, 'tms' is fully implicitly treated. !
465	! --------------------------------------------------------------------------------------- !
466
467	call vd_lu_decomp( pcols , pver , ncol , &
468	ksrf , kvm , tmpi2 , rpdel , ztodt , zero , &
469	ca , cc , dnom , tmpm , ntop , nbot )
470
471	call vd_lu_solve( pcols , pver , ncol , &
472	u , ca , tmpm , dnom , ntop , nbot , zero )
473
474	call vd_lu_solve( pcols , pver , ncol , &
475	v , ca , tmpm , dnom , ntop , nbot , zero )
476
477	! ---------------------------------------------------------------------- !
478	! Calculate 'total' ( tautotx ) and 'tms' ( tautmsx ) stresses that !
479	! have been actually added into the atmosphere at the current time step. !
480	! Also, update residual stress, if required. !
481	! ---------------------------------------------------------------------- !
482
483	do i = 1, ncol
484
485	! Compute the implicit 'tms' using the updated winds.
486	! Below 'tautmsx(i),tautmsy(i)' are pure implicit mountain stresses
487	! that has been actually added into the atmosphere both for explicit
488	! and implicit approach.
489
490	tautmsx(i) = -ksrftms(i)*u(i,pver)
491	tautmsy(i) = -ksrftms(i)*v(i,pver)
492
493	if( do_iss ) then
494
495	! Compute vertical integration of final horizontal momentum
496
497	usum_out(i) = 0._r8
498	vsum_out(i) = 0._r8
499	do k = 1, pver
500	usum_out(i) = usum_out(i) + (1._r8/gravit)*u(i,k)/rpdel(i,k)
501	vsum_out(i) = vsum_out(i) + (1._r8/gravit)*v(i,k)/rpdel(i,k)
502	end do
503
504	! Compute net stress added into the atmosphere at the current time step.
505	! Note that the difference between 'usum_in' and 'usum_out' are induced
506	! by 'explicit residual stress + implicit total stress' for implicit case, while
507	! by 'explicit normal stress + implicit tms stress' for explicit case.
508	! Here, 'tautotx(i)' is net stress added into the air at the current time step.
509
510	tauimpx(i) = ( usum_out(i) - usum_in(i) ) / ztodt
511	tauimpy(i) = ( vsum_out(i) - vsum_in(i) ) / ztodt
512
513	tautotx(i) = tauimpx(i)
514	tautoty(i) = tauimpy(i)
515
516	! Compute redisual stress and update if required.
517	! Note that the total stress we should have added at the current step is
518	! the sum of 'taux(i) - ksrftms(i)*u(i,pver) + tauresx(i)'.
519
520	if( itaures .eq. 1 ) then
521	tauresx(i) = taux(i) + tautmsx(i) + tauresx(i) - tauimpx(i)
522	tauresy(i) = tauy(i) + tautmsy(i) + tauresy(i) - tauimpy(i)
523	endif
524
525	else
526
527	tautotx(i) = tautmsx(i) + taux(i)
528	tautoty(i) = tautmsy(i) + tauy(i)
529	tauresx(i) = 0._r8
530	tauresy(i) = 0._r8
531
532	end if ! End of 'do_iss' routine
533
534	end do ! End of 'do i = 1, ncol' routine
535
536	! Debug
537	! icol = phys_debug_col(lchnk)
538	! if ( icol > 0 .and. get_nstep() .ge. 1 ) then
539	! write(iulog,*)
540	! write(iulog,*) 'diffusion_solver debug'
541	! write(iulog,*)
542	! write(iulog,*) 'u_in, u_res, u_out'
543	! write(iulog,*) u_in, u_res, u(icol,pver)
544	! write(iulog,*) 'tauresx_in, tautmsx, tauimpx(actual), tauimpx(derived), tauresx_out, taux'
545	! write(iulog,) tauresx_in, tautmsx(icol), tauimpx(icol), -ksrf(icol)u(icol,pver), tauresx(icol), taux(icol)
546	! write(iulog,*)
547	! write(iulog,*) 'v_in, v_res, v_out'
548	! write(iulog,*) v_in, v_res, v(icol,pver)
549	! write(iulog,*) 'tauresy_in, tautmsy, tauimpy(actual), tauimpy(derived), tauresy_out, tauy'
550	! write(iulog,) tauresy_in, tautmsy(icol), tauimpy(icol), -ksrf(icol)v(icol,pver), tauresy(icol), tauy(icol)
551	! write(iulog,*)
552	! write(iulog,*) 'itaures, ksrf, ksrfturb, ksrftms'
553	! write(iulog,*) itaures, ksrf(icol), ksrfturb(icol), ksrftms(icol)
554	! write(iulog,*)
555	! endif
556	! Debug
557
558	! ------------------------------------ !
559	! Calculate kinetic energy dissipation !
560	! ------------------------------------ !
561
562	! Modification : In future, this should be set exactly same as
563	! the ones in the convection schemes
564
565	! 1. Compute dissipation term at interfaces
566	! Note that 'u,v' are already diffused wind, and 'tautotx,tautoty' are
567	! implicit stress that has been actually added. On the other hand,
568	! 'dinp_u, dinp_v' were computed using non-diffused input wind.
569
570	! Modification : I should check whether non-consistency between 'u' and 'dinp_u'
571	! is correctly intended approach. I think so.
572
573	k = pver + 1
574	do i = 1, ncol
575	tmpi1(i,1) = 0._r8
576	tmpi1(i,k) = 0.5_r8 * ztodt * gravit * &
577	( (-u(i,k-1) + dinp_u(i,k))tautotx(i) + (-v(i,k-1) + dinp_v(i,k))tautoty(i) )
578	end do
579
580	do k = 2, pver
581	do i = 1, ncol
582	dout_u = u(i,k) - u(i,k-1)
583	dout_v = v(i,k) - v(i,k-1)
584	tmpi1(i,k) = 0.25_r8 * tmpi2(i,k) * kvm(i,k) * &
585	( dout_u2 + dout_v2 + dout_udinp_u(i,k) + dout_vdinp_v(i,k) )
586	end do
587	end do
588
589	! 2. Compute dissipation term at midpoints, add to dry static energy
590
591	do k = 1, pver
592	do i = 1, ncol
593	dtk(i,k) = ( tmpi1(i,k+1) + tmpi1(i,k) ) * rpdel(i,k)
594	dse(i,k) = dse(i,k) + dtk(i,k)
595	end do
596	end do
597
598	end if ! End of diffuse horizontal momentum, diffuse(fieldlist,'u') routine
599
600	!-------------------------- !
601	! Diffuse Dry Static Energy !
602	!-------------------------- !
603
604	! Modification : In future, we should diffuse the fully conservative
605	! moist static energy,not the dry static energy.
606
607	if( diffuse(fieldlist,'s') ) then
608
609	! Add counter-gradient to input static energy profiles
610
611	do k = 1, pver
612	dse(:ncol,k) = dse(:ncol,k) + ztodt * rpdel(:ncol,k) * gravit * &
613	( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgh(:ncol,k+1) &
614	- rhoi(:ncol,k ) * kvh(:ncol,k ) * cgh(:ncol,k ) )
615	end do
616
617	! Add the explicit surface fluxes to the lowest layer
618
619	dse(:ncol,pver) = dse(:ncol,pver) + tmp1(:ncol) * shflx(:ncol)
620
621	! Diffuse dry static energy
622
623	call vd_lu_decomp( pcols , pver , ncol , &
624	zero , kvh , tmpi2 , rpdel , ztodt , cc_top, &
625	ca , cc , dnom , tmpm , ntop , nbot )
626
627	call vd_lu_solve( pcols , pver , ncol , &
628	dse , ca , tmpm , dnom , ntop , nbot , cd_top )
629
630	! Calculate flux at top interface
631
632	! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
633
634	if( do_molec_diff ) then
635	topflx(:ncol) = - kvh(:ncol,ntop_molec) * tmpi2(:ncol,ntop_molec) / (ztodtgravit) &
636	( dse(:ncol,ntop_molec) - dse_top(:ncol) )
637	end if
638
639	endif
640
641	!---------------------------- !
642	! Diffuse Water Vapor Tracers !
643	!---------------------------- !
644
645	! Modification : For aerosols, I need to use separate treatment
646	! for aerosol mass and aerosol number.
647
648	! Loop through constituents
649
650	need_decomp = .true.
651
652	do m = 1, ncnst
653
654	if( diffuse(fieldlist,'q',m) ) then
655
656	! Add the nonlocal transport terms to constituents in the PBL.
657	! Check for neg q's in each constituent and put the original vertical
658	! profile back if a neg value is found. A neg value implies that the
659	! quasi-equilibrium conditions assumed for the countergradient term are
660	! strongly violated.
661
662	qtm(:ncol,:pver) = q(:ncol,:pver,m)
663
664	do k = 1, pver
665	q(:ncol,k,m) = q(:ncol,k,m) + &
666	ztodt * rpdel(:ncol,k) * gravit * ( cflx(:ncol,m) * rrho(:ncol) ) * &
667	( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgs(:ncol,k+1) &
668	- rhoi(:ncol,k ) * kvh(:ncol,k ) * cgs(:ncol,k ) )
669	end do
670	lqtst(:ncol) = all(q(:ncol,1:pver,m) >= qmincg(m), 2)
671	do k = 1, pver
672	q(:ncol,k,m) = merge( q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol) )
673	end do
674
675	! Add the explicit surface fluxes to the lowest layer
676
677	q(:ncol,pver,m) = q(:ncol,pver,m) + tmp1(:ncol) * cflx(:ncol,m)
678
679	! Diffuse constituents.
680
681	if( need_decomp ) then
682
683	call vd_lu_decomp( pcols , pver , ncol , &
684	zero , kvq , tmpi2 , rpdel , ztodt , zero , &
685	ca , cc , dnom , tmpm , ntop , nbot )
686
687	if( do_molec_diff ) then
688
689	! Update decomposition in molecular diffusion range, include separation velocity term
690
691	status = vd_lu_qdecomp( pcols , pver , ncol , cnst_fixed_ubc(m), cnst_mw(m), ubc_mmr(:,m), &
692	kvq , kq_scal, mw_fac(m) , tmpi2 , rpdel , &
693	ca , cc , dnom , tmpm , rhoi , &
694	tint , ztodt , ntop_molec, nbot_molec , cd_top )
695	else
696	need_decomp = .false.
697	endif
698	end if
699
700	call vd_lu_solve( pcols , pver , ncol , &
701	q(1,1,m) , ca, tmpm , dnom , ntop , nbot , cd_top )
702	end if
703	end do
704
705	return
706	end subroutine compute_vdiff
707
708	! =============================================================================== !
709	! !
710	! =============================================================================== !
711
712	subroutine vd_lu_decomp( pcols, pver, ncol , &
713	ksrf , kv , tmpi , rpdel, ztodt , cc_top, &
714	ca , cc , dnom , ze , ntop , nbot )
715	!---------------------------------------------------------------------- !
716	! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the !
717	! tridiagonal diffusion matrix. !
718	! The diagonal elements (1+ca(k)+cc(k)) are not required by the solver. !
719	! Also determine ze factor and denominator for ze and zf (see solver). !
720	!---------------------------------------------------------------------- !
721
722	! --------------------- !
723	! Input-Output Argument !
724	! --------------------- !
725
726	integer, intent(in) :: pcols ! Number of allocated atmospheric columns
727	integer, intent(in) :: pver ! Number of allocated atmospheric levels
728	integer, intent(in) :: ncol ! Number of computed atmospheric columns
729	integer, intent(in) :: ntop ! Top level to operate on
730	integer, intent(in) :: nbot ! Bottom level to operate on
731	real(r8), intent(in) :: ksrf(pcols) ! Surface "drag" coefficient [ kg/s/m2 ]
732	real(r8), intent(in) :: kv(pcols,pver+1) ! Vertical diffusion coefficients [ m2/s ]
733	real(r8), intent(in) :: tmpi(pcols,pver+1) ! dt(g/R)2/dppi(k+1)/(.5(tm(k+1)+tm(k))*2
734	real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel (thickness bet interfaces)
735	real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
736	real(r8), intent(in) :: cc_top(pcols) ! Lower diagonal on top interface (for fixed ubc only)
737
738	real(r8), intent(out) :: ca(pcols,pver) ! Upper diagonal
739	real(r8), intent(out) :: cc(pcols,pver) ! Lower diagonal
740	real(r8), intent(out) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
741	real(r8), intent(out) :: ze(pcols,pver) ! Term in tri-diag. matrix system
742
743	! --------------- !
744	! Local Variables !
745	! --------------- !
746
747	integer :: i ! Longitude index
748	integer :: k ! Vertical index
749
750	! ----------------------- !
751	! Main Computation Begins !
752	! ----------------------- !
753
754	! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the
755	! tridiagonal diffusion matrix. The diagonal elements (cb=1+ca+cc) are
756	! a combination of ca and cc; they are not required by the solver.
757
758	do k = nbot - 1, ntop, -1
759	do i = 1, ncol
760	ca(i,k ) = kv(i,k+1) * tmpi(i,k+1) * rpdel(i,k )
761	cc(i,k+1) = kv(i,k+1) * tmpi(i,k+1) * rpdel(i,k+1)
762	end do
763	end do
764
765	! The bottom element of the upper diagonal (ca) is zero (element not used).
766	! The subdiagonal (cc) is not needed in the solver.
767
768	do i = 1, ncol
769	ca(i,nbot) = 0._r8
770	end do
771
772	! Calculate e(k). This term is
773	! required in solution of tridiagonal matrix defined by implicit diffusion eqn.
774
775	do i = 1, ncol
776	dnom(i,nbot) = 1._r8/(1._r8 + cc(i,nbot) + ksrf(i)ztodtgravit*rpdel(i,nbot))
777	ze(i,nbot) = cc(i,nbot)*dnom(i,nbot)
778	end do
779
780	do k = nbot - 1, ntop + 1, -1
781	do i = 1, ncol
782	dnom(i,k) = 1._r8/(1._r8 + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1))
783	ze(i,k) = cc(i,k)*dnom(i,k)
784	end do
785	end do
786
787	do i = 1, ncol
788	dnom(i,ntop) = 1._r8/(1._r8 + ca(i,ntop) + cc_top(i) - ca(i,ntop)*ze(i,ntop+1))
789	end do
790
791	return
792	end subroutine vd_lu_decomp
793
794	! =============================================================================== !
795	! !
796	! =============================================================================== !
797
798	subroutine vd_lu_solve( pcols , pver , ncol , &
799	q , ca , ze , dnom , ntop , nbot , cd_top )
800	!----------------------------------------------------------------------------------- !
801	! Solve the implicit vertical diffusion equation with zero flux boundary conditions. !
802	! Procedure for solution of the implicit equation follows Richtmyer and !
803	! Morton (1967,pp 198-200). !
804	! !
805	! The equation solved is !
806	! !
807	! -ca(k)q(k+1) + cb(k)q(k) - cc(k)*q(k-1) = d(k), !
808	! !
809	! where d(k) is the input profile and q(k) is the output profile !
810	! !
811	! The solution has the form !
812	! !
813	! q(k) = ze(k)*q(k-1) + zf(k) !
814	! !
815	! ze(k) = cc(k) * dnom(k) !
816	! !
817	! zf(k) = [d(k) + ca(k)zf(k+1)] dnom(k) !
818	! !
819	! dnom(k) = 1/[cb(k) - ca(k)ze(k+1)] = 1/[1 + ca(k) + cc(k) - ca(k)ze(k+1)] !
820	! !
821	! Note that the same routine is used for temperature, momentum and tracers, !
822	! and that input variables are replaced. !
823	! ---------------------------------------------------------------------------------- !
824
825	! --------------------- !
826	! Input-Output Argument !
827	! --------------------- !
828
829	integer, intent(in) :: pcols ! Number of allocated atmospheric columns
830	integer, intent(in) :: pver ! Number of allocated atmospheric levels
831	integer, intent(in) :: ncol ! Number of computed atmospheric columns
832	integer, intent(in) :: ntop ! Top level to operate on
833	integer, intent(in) :: nbot ! Bottom level to operate on
834	real(r8), intent(in) :: ca(pcols,pver) ! -Upper diag coeff.of tri-diag matrix
835	real(r8), intent(in) :: ze(pcols,pver) ! Term in tri-diag solution
836	real(r8), intent(in) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - ca(k)*ze(k+1))
837	real(r8), intent(in) :: cd_top(pcols) ! cc_top * ubc value
838
839	real(r8), intent(inout) :: q(pcols,pver) ! Constituent field
840
841	! --------------- !
842	! Local Variables !
843	! --------------- !
844
845	real(r8) :: zf(pcols,pver) ! Term in tri-diag solution
846	integer i, k ! Longitude, vertical indices
847
848	! ----------------------- !
849	! Main Computation Begins !
850	! ----------------------- !
851
852	! Calculate zf(k). Terms zf(k) and ze(k) are required in solution of
853	! tridiagonal matrix defined by implicit diffusion equation.
854	! Note that only levels ntop through nbot need be solved for.
855
856	do i = 1, ncol
857	zf(i,nbot) = q(i,nbot)*dnom(i,nbot)
858	end do
859
860	do k = nbot - 1, ntop + 1, -1
861	do i = 1, ncol
862	zf(i,k) = (q(i,k) + ca(i,k)zf(i,k+1))dnom(i,k)
863	end do
864	end do
865
866	! Include boundary condition on top element
867
868	k = ntop
869	do i = 1, ncol
870	zf(i,k) = (q(i,k) + cd_top(i) + ca(i,k)zf(i,k+1))dnom(i,k)
871	end do
872
873	! Perform back substitution
874
875	do i = 1, ncol
876	q(i,ntop) = zf(i,ntop)
877	end do
878
879	do k = ntop + 1, nbot, +1
880	do i = 1, ncol
881	q(i,k) = zf(i,k) + ze(i,k)*q(i,k-1)
882	end do
883	end do
884
885	return
886	end subroutine vd_lu_solve
887
888	! =============================================================================== !
889	! !
890	! =============================================================================== !
891
892	character(128) function vdiff_select( fieldlist, name, qindex )
893	! --------------------------------------------------------------------- !
894	! This function sets the field with incoming name as one to be diffused !
895	! --------------------------------------------------------------------- !
896	type(vdiff_selector), intent(inout) :: fieldlist
897	character(*), intent(in) :: name
898	integer, intent(in), optional :: qindex
899
900	vdiff_select = ''
901	select case (name)
902	case ('u','U')
903	fieldlist%fields(1) = .true.
904	case ('v','V')
905	fieldlist%fields(2) = .true.
906	case ('s','S')
907	fieldlist%fields(3) = .true.
908	case ('q','Q')
909	if( present(qindex) ) then
910	fieldlist%fields(3 + qindex) = .true.
911	else
912	fieldlist%fields(4) = .true.
913	endif
914	case default
915	write(vdiff_select,*) 'Bad argument to vdiff_index: ', name
916	end select
917	return
918
919	end function vdiff_select
920
921	type(vdiff_selector) function not(a)
922	! ------------------------------------------------------------- !
923	! This function extends .not. to operate on type vdiff_selector !
924	! ------------------------------------------------------------- !
925	type(vdiff_selector), intent(in) :: a
926	allocate(not%fields(size(a%fields)))
927	not%fields(:) = .not. a%fields(:)
928	end function not
929
930	logical function my_any(a)
931	! -------------------------------------------------- !
932	! This function extends the intrinsic function 'any' !
933	! to operate on type vdiff_selector !
934	! -------------------------------------------------- !
935	type(vdiff_selector), intent(in) :: a
936	my_any = any(a%fields)
937	end function my_any
938
939	logical function diffuse(fieldlist,name,qindex)
940	! ---------------------------------------------------------------------------- !
941	! This function reports whether the field with incoming name is to be diffused !
942	! ---------------------------------------------------------------------------- !
943	type(vdiff_selector), intent(in) :: fieldlist
944	character(*), intent(in) :: name
945	integer, intent(in), optional :: qindex
946
947	select case (name)
948	case ('u','U')
949	diffuse = fieldlist%fields(1)
950	case ('v','V')
951	diffuse = fieldlist%fields(2)
952	case ('s','S')
953	diffuse = fieldlist%fields(3)
954	case ('q','Q')
955	if( present(qindex) ) then
956	diffuse = fieldlist%fields(3 + qindex)
957	else
958	diffuse = fieldlist%fields(4)
959	endif
960	case default
961	diffuse = .false.
962	end select
963	return
964	end function diffuse
965
966	end module diffusion_solver

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

Download in other formats:

Original Format