Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_boulac.F @ 1562

Last change on this file since 1562 was 1, checked in by lfita, 11 years ago

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

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 37.2 KB

Line
1	MODULE module_bl_boulac
2
3	!USE module_model_constants
4
5
6	!------------------------------------------------------------------------
7	! Calculation of the tendency due to momentum, heat
8	! and moisture turbulent fluxes follwing the approach
9	! of Bougeault and Lacarrere, 1989 (MWR, 117, 1872-1890).
10	! The scheme computes a prognostic ecuation for TKE and derives
11	! dissipation and turbulent coefficients using length scales.
12	!
13	! Subroutine written by Alberto Martilli, CIEMAT, Spain,
14	! e-mail:alberto_martilli@ciemat.es
15	! August 2006.
16	!------------------------------------------------------------------------
17	! IN THIS VERSION TKE IS NOT ADVECTED!!!!
18	! TO BE CHANGED IN THE FUTURE
19	!
20	! -----------------------------------------------------------------------
21	! Constant used in the module
22	! ck_b=constant used in the compuation of diffusion coefficients
23	! ceps_b=constant used inthe computation of dissipation
24	! temin= minimum value allowed for TKE
25	! vk=von karman constant
26	! -----------------------------------------------------------------------
27
28	real ck_b,ceps_b,vk,temin ! constant for Bougeault and Lacarrere
29	parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4) ! impose minimum values for tke similar to those of MYJ
30	! -----------------------------------------------------------------------
31
32
33	CONTAINS
34
35	subroutine boulac(frc_urb2d,idiff,flag_bep,dz8w,dt,u_phy,v_phy &
36	,th_phy,rho,qv_curr,hfx &
37	,qfx,ustar,cp,g &
38	,rublten,rvblten,rthblten &
39	,rqvblten,rqcblten &
40	,tke,dlk,wu,wv,wt,wq,exch_h,exch_m,pblh &
41	,a_u_bep,a_v_bep,a_t_bep,a_q_bep &
42	,a_e_bep,b_u_bep,b_v_bep &
43	,b_t_bep,b_q_bep,b_e_bep,dlg_bep &
44	,dl_u_bep,sf_bep,vl_bep &
45	,ids,ide, jds,jde, kds,kde &
46	,ims,ime, jms,jme, kms,kme &
47	,its,ite, jts,jte, kts,kte)
48
49	implicit none
50
51
52
53	!-----------------------------------------------------------------------
54	! Input
55	!------------------------------------------------------------------------
56	INTEGER:: ids,ide, jds,jde, kds,kde, &
57	ims,ime, jms,jme, kms,kme, &
58	its,ite, jts,jte, kts,kte
59
60	integer, INTENT(IN) :: idiff
61	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: DZ8W !vertical resolution
62	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: qv_curr !moisture
63	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: RHO !air density
64	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: TH_PHY !potential temperature
65	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: U_PHY !x-component of wind
66	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: V_PHY !y-component of wind
67	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: ustar !friction velocity
68	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: hfx !sensible heat flux (W/m2) at surface
69	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: qfx !moisture flux at surface
70	real, INTENT(IN ) :: g,cp !gravity and Cp
71	REAL, INTENT(IN ):: DT ! Time step
72
73	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: FRC_URB2D !fraction cover urban
74	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH !PBL height
75	!
76	! variable added for urban
77	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_u_bep ! Implicit component for the momemtum in X-direction
78	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_v_bep ! Implicit component for the momemtum in Y-direction
79	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_t_bep ! Implicit component for the Pot. Temp.
80	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_q_bep ! Implicit component for Moisture
81	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_e_bep ! Implicit component for the TKE
82	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_u_bep ! Explicit component for the momemtum in X-direction
83	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_v_bep ! Explicit component for the momemtum in Y-direction
84	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_t_bep ! Explicit component for the Pot. Temp.
85	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_q_bep ! Explicit component for Moisture
86	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_e_bep ! Explicit component for the TKE
87
88	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep ! Height above ground (L_ground in formula (24) of the BLM paper).
89	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::dl_u_bep ! Length scale (lb in formula (22) ofthe BLM paper).
90	! urban surface and volumes
91	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::sf_bep ! surface of the urban grid cells
92	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::vl_bep ! volume of the urban grid cells
93	LOGICAL, INTENT(IN) :: flag_bep !flag for BEP
94
95	!
96	!-----------------------------------------------------------------------
97	! Local, carried on from one timestep to the other
98	!------------------------------------------------------------------------
99	! real, save, allocatable, dimension (:,:,:)::TKE ! Turbulent kinetic energy
100	real, dimension (ims:ime, kms:kme, jms:jme) ::th_0 ! reference state for potential temperature
101
102	!------------------------------------------------------------------------
103	! Output
104	!------------------------------------------------------------------------
105	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_h ! exchange coefficient for heat
106	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_m ! exchange coefficient for momentum
107	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: tke ! Turbulence Kinetic Energy
108	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wu ! Turbulent flux of momentum (x)
109	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wv ! Turbulent flux of momentum (y)
110	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wt ! Turbulent flux of temperature
111	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wq ! Turbulent flux of water vapor
112	real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: dlk ! Turbulent flux of water vapor
113	! only if idiff not equal 1:
114	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RUBLTEN !tendency for U_phy
115	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RVBLTEN !tendency for V_phy
116	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RTHBLTEN !tendency for TH_phy
117	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVBLTEN !tendency for QV_curr
118	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQCBLTEN !tendency for QV_curr
119
120	!--------------------------------------------------------------
121	! Local
122	!--------------------------------------------------------------
123	! 1D array used for the input and output of the routine boulac1D
124
125	real z1D(kms:kme) ! vertical coordinates (faces of the grid)
126	real dz1D(kms:kme) ! vertical resolution
127	real u1D(kms:kme) ! wind speed in the x directions
128	real v1D(kms:kme) ! wind speed in the y directions
129	real th1D(kms:kme) ! potential temperature
130	real q1D(kms:kme) ! potential temperature
131	real rho1D(kms:kme) ! air density
132	real rhoz1D(kms:kme) ! air density at the faces
133	real tke1D(kms:kme) ! air pressure
134	real th01D(kms:kme) ! reference potential temperature
135	real dlk1D(kms:kme) ! dlk
136	real dls1D(kms:kme) ! dls
137	real exch1D(kms:kme) ! exch
138	real sf1D(kms:kme) ! surface of the grid cells
139	real vl1D(kms:kme) ! volume of the grid cells
140	real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
141	real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
142	real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
143	real a_q1D(kms:kme) ! Implicit component of the moisture sources or sinks
144	real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
145	real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
146	real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
147	real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
148	real b_q1D(kms:kme) ! Explicit component of the moisture sources or sinks
149	real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
150	real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
151	real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
152	real sh1D(kms:kme) ! shear
153	real bu1D(kms:kme) ! buoyancy
154	real wu1D(kms:kme) ! turbulent flux of momentum (x component)
155	real wv1D(kms:kme) ! turbulent flux of momentum (y component)
156	real wt1D(kms:kme) ! turbulent flux of temperature
157	real wq1D(kms:kme) ! turbulent flux of water vapor
158	! local added only for diagnostic output
159	real a_e(ims:ime,kms:kme,jms:jme) ! implicit term in TKE
160	real b_e(ims:ime,kms:kme,jms:jme) ! explicit term in TKE
161	real bu(ims:ime,kms:kme,jms:jme) ! buoyancy term in TKE
162	real sh(ims:ime,kms:kme,jms:jme) ! shear term in TKE
163	real wrk(ims:ime) ! working array
164	integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1
165	real ufrac_int ! urban fraction
166	real vect,time_tke,hour,zzz
167	real ustarf
168	real summ1,summ2,summ3
169	save time_tke,hour
170	!
171	! store reference state for potential temperature
172	! and initialize tke
173	!
174
175	!here I fix the value of the reference state equal to the value of the potnetial temperature
176	! the only use of this variable in the code is to compute the paramter BETA = g/th0
177	! I fix it to 300K.
178
179	do ix=its,ite
180	do iy=jts,jte
181	do iz=kts,kte
182	! th_0(ix,iz,iy)=th_phy(ix,iz,iy)
183	th_0(ix,iz,iy)=300.
184	enddo
185	enddo
186	enddo
187
188	bu1d=0.
189	sh1d=0.
190	b_e1d=0.
191	b_u1d=0.
192	b_v1d=0.
193	b_t1d=0.
194	b_q1d=0.
195	a_e1d=0.
196	a_u1d=0.
197	a_v1d=0.
198	a_t1d=0.
199	a_q1d=0.
200	z1D=0. ! vertical coordinates (faces of the grid)
201	dz1D=0. ! vertical resolution
202	u1D =0. ! wind speed in the x directions
203	v1D =0. ! wind speed in the y directions
204	th1D=0. ! potential temperature
205	q1D=0. ! potential temperature
206	rho1D=0. ! air density
207	rhoz1D=0. ! air density at the faces
208	tke1D =0. ! air pressure
209	th01D =0. ! reference potential temperature
210	dlk1D =0. ! dlk
211	dls1D =0. ! dls
212	exch1D=0. ! exch
213	sf1D =1. ! surface of the grid cells
214	vl1D =1. ! volume of the grid cells
215	a_u1D =0. ! Implicit component of the momentum sources or sinks in the X-direction
216	a_v1D =0. ! Implicit component of the momentum sources or sinks in the Y-direction
217	a_t1D =0. ! Implicit component of the heat sources or sinks
218	a_q1D =0. ! Implicit component of the moisture sources or sinks
219	a_e1D =0. ! Implicit component of the TKE sources or sinks
220	b_u1D =0. ! Explicit component of the momentum sources or sinks in the X-direction
221	b_v1D =0. ! Explicit component of the momentum sources or sinks in the Y-direction
222	b_t1D =0. ! Explicit component of the heat sources or sinks
223	b_q1D =0. ! Explicit component of the moisture sources or sinks
224	b_e1D =0. ! Explicit component of the TKE sources or sinks
225	dlg1D =0. ! Height above ground (L_ground in formula (24) of the BLM paper).
226	dl_u1D=0. ! Length scale (lb in formula (22) ofthe BLM paper)
227	sh1D =0. ! shear
228	bu1D =0. ! buoyancy
229	wu1D =0. ! turbulent flux of momentum (x component)
230	wv1D =0. ! turbulent flux of momentum (y component)
231	wt1D =0. ! turbulent flux of temperature
232	wq1D =0. ! turbulent flux of water vapor
233	! local added only for diagnostic output
234
235	! loop over the columns.
236	! put variables in 1D temporary arrays
237	!
238
239	do ix=its,ite
240	do iy=jts,jte
241	z1d(kts)=0.
242	do iz= kts,kte
243	u1D(iz)=u_phy(ix,iz,iy)
244	v1D(iz)=v_phy(ix,iz,iy)
245	th1D(iz)=th_phy(ix,iz,iy)
246	q1D(iz)=qv_curr(ix,iz,iy)
247	tke1D(iz)=tke(ix,iz,iy)
248	rho1D(iz)=rho(ix,iz,iy)
249	th01D(iz)=th_0(ix,iz,iy)
250	dz1D(iz)=dz8w(ix,iz,iy)
251	z1D(iz+1)=z1D(iz)+dz1D(iz)
252	enddo
253	rhoz1D(kts)=rho1D(kts)
254	do iz=kts+1,kte
255	rhoz1D(iz)=(rho1D(iz)dz1D(iz-1)+rho1D(iz-1)dz1D(iz))/(dz1D(iz-1)+dz1D(iz))
256	enddo
257	rhoz1D(kte+1)=rho1D(kte)
258	if(flag_bep)then
259	do iz=kts,kte
260	a_e1D(iz)=a_e_bep(ix,iz,iy)
261	b_e1D(iz)=b_e_bep(ix,iz,iy)
262	dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.(1.-frc_urb2d(ix,iy))+dlg_bep(ix,iz,iy)frc_urb2d(ix,iy)
263	dl_u1D(iz)=dl_u_bep(ix,iz,iy)
264	if((1.-frc_urb2d(ix,iy)).lt.1.)dl_u1D(iz)=dl_u1D(iz)/frc_urb2d(ix,iy)
265	vl1D(iz)=vl_bep(ix,iz,iy)
266	sf1D(iz)=sf_bep(ix,iz,iy)
267	enddo
268	ufrac_int=frc_urb2d(ix,iy)
269	sf1D(kte+1)=sf_bep(ix,1,iy)
270	else
271	do iz=kts,kte
272	a_e1D(iz)=0.
273	b_e1D(iz)=0.
274	dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.
275	dl_u1D(iz)=0.
276	vl1D(iz)=1.
277	sf1D(iz)=1.
278	enddo
279	ufrac_int=0.
280	sf1D(kte+1)=1.
281	endif
282	! call the routine that will solve the turbulence in 1D for tke
283
284	call boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz1d,z1D,dt,u1D,v1D,th1D,rho1D,rhoz1D,q1D,th01D,&
285	tke1D,ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),cp,g, &
286	a_e1D,b_e1D, &
287	dlg1D,dl_u1D,sf1D,vl1D,dlk1D,dls1D,exch1D,sh1D,bu1D)
288
289	call pbl_height(kms,kme,kts,kte,dz1d,z1d,th1D,q1D,pblh(ix,iy))
290
291
292	! store turbulent exchange coefficients, TKE, and other variables
293
294	do iz= kts,kte
295	a_e(ix,iz,iy)=a_e1D(iz)
296	b_e(ix,iz,iy)=b_e1D(iz)
297	if(flag_bep)then
298	dlg_bep(ix,iz,iy)=dlg1D(iz)
299	endif
300	tke(ix,iz,iy)=tke1D(iz)
301	dlk(ix,iz,iy)=dlk1D(iz)
302	sh(ix,iz,iy)=sh1D(iz)
303	bu(ix,iz,iy)=bu1D(iz)
304	exch_h(ix,iz,iy)=exch1D(iz)
305	exch_m(ix,iz,iy)=exch1D(iz)
306	enddo
307
308	if(idiff.ne.1)then
309
310	! estimate the tendencies
311
312	if(flag_bep)then
313	do iz=kts,kte
314	a_t1D(iz)=a_t_bep(ix,iz,iy)
315	b_t1D(iz)=b_t_bep(ix,iz,iy)
316	a_u1D(iz)=a_u_bep(ix,iz,iy)
317	b_u1D(iz)=b_u_bep(ix,iz,iy)
318	a_v1D(iz)=a_v_bep(ix,iz,iy)
319	b_v1D(iz)=b_v_bep(ix,iz,iy)
320	a_q1D(iz)=a_q_bep(ix,iz,iy)
321	b_q1D(iz)=b_q_bep(ix,iz,iy)
322	enddo
323	else
324	do iz=kts,kte
325	a_t1D(iz)=0.
326	b_t1D(iz)=0.
327	a_u1D(iz)=0.
328	b_u1D(iz)=0.
329	a_v1D(iz)=0.
330	b_v1D(iz)=0.
331	a_q1D(iz)=0.
332	b_q1D(iz)=0.
333	enddo
334	b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp
335	b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1)
336	a_u1D(1)=(-ustar(ix,iy)ustar(ix,iy)/dz1D(1)/((u1D(1)2.+v1D(1)2.)*.5))
337	a_v1D(1)=(-ustar(ix,iy)ustar(ix,iy)/dz1D(1)/((u1D(1)2.+v1D(1)2.)*.5))
338	endif
339
340	! solve diffusion equation for momentum x component
341	call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D)
342
343	! solve diffusion equation for momentum y component
344	call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D)
345
346	! solve diffusion equation for potential temperature
347	call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D)
348
349	! solve diffusion equation for water vapor mixing ratio
350	call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D)
351
352	! compute the tendencies
353
354	do iz= kts,kte
355	rthblten(ix,iz,iy)=rthblten(ix,iz,iy)+(th1D(iz)-th_phy(ix,iz,iy))/dt
356	rqvblten(ix,iz,iy)=rqvblten(ix,iz,iy)+(q1D(iz)-qv_curr(ix,iz,iy))/dt
357	rublten(ix,iz,iy)=rublten(ix,iz,iy)+(u1D(iz)-u_phy(ix,iz,iy))/dt
358	rvblten(ix,iz,iy)=rvblten(ix,iz,iy)+(v1D(iz)-v_phy(ix,iz,iy))/dt
359	wu(ix,iz,iy)=wu1D(iz)
360	wv(ix,iz,iy)=wv1D(iz)
361	wt(ix,iz,iy)=wt1D(iz)
362	wq(ix,iz,iy)=wq1D(iz)
363	enddo
364	endif
365
366	enddo ! iy
367	enddo ! ix
368
369
370	time_tke=time_tke+dt
371
372
373	return
374	end subroutine boulac
375
376
377	! ===6=8===============================================================72
378
379	! ===6=8===============================================================72
380
381	subroutine boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz,z,dt,u,v,th,rho,rhoz,qa,th0,te, &
382	ustar,hfx,qfx,cp,g, &
383	a_e,b_e, &
384	dlg,dl_u,sf,vl,dlk,dls,exch,sh,bu)
385
386	! ----------------------------------------------------------------------
387	! 1D resolution of TKE following Bougeault and Lacarrere
388	! ----------------------------------------------------------------------
389
390	implicit none
391
392	integer iz,ix,iy
393
394	! ----------------------------------------------------------------------
395	! INPUT:
396	! ----------------------------------------------------------------------
397
398
399	integer kms,kme,kts,kte
400	real z(kms:kme) ! Altitude above the ground of the cell interfaces.
401	real dz(kms:kme) ! vertical resolution
402	real u(kms:kme) ! Wind speed in the x direction
403	real v(kms:kme) ! Wind speed in the y direction
404	real th(kms:kme) ! Potential temperature
405	real rho(kms:kme) ! Air density
406	real g ! gravity
407	real cp !
408	real te(kms:kme) ! turbulent kinetic energy
409	real qa(kms:kme) ! air humidity
410	real th0(kms:kme) ! Reference potential temperature
411	real dt ! Time step
412	real ustar ! ustar
413	real hfx ! sensbile heat flux
414	real qfx ! kinematic latent heat flux
415	real sf(kms:kme) ! surface of the urban grid cells
416	real vl(kms:kme) ! volume of the urban grid cells
417	real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
418	real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
419	real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
420	real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
421	real ufrac_int ! urban fraction
422	! local variables not needed in principle, but that can be used to estimate the budget and turbulent fluxes
423
424	real we(kms:kme),dwe(kms:kme)
425
426	! local variables
427	real sh(kms:kme) ! shear term in TKE eqn.
428	real bu(kms:kme) ! buoyancy term in TKE eqn.
429	real td(kms:kme) ! dissipation term in TKE eqn.
430	real exch(kms:kme) ! turbulent diffusion coefficients (defined at the faces)
431	real dls(kms:kme) ! dissipation length scale
432	real dlk(kms:kme) ! length scale used to estimate exch
433	real dlu(kms:kme) ! l_up
434	real dld(kms:kme) ! l_down
435	real rhoz(kms:kme) !air density at the faces of the cell
436	real tstar ! derived from hfx and ustar
437	real beta
438	real summ1,summ2,summ3,summ4
439	! interpolate air density at the faces
440
441
442
443	! estimation of tstar
444
445	tstar=-hfx/rho(1)/cp/ustar
446
447	! first compute values of dlu and dld (length scales up and down).
448
449	call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
450
451	!then average them to obtain dls and dlk (length scales for dissipation and eddy coefficients)
452
453	call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
454
455	! compute the turbulent diffusion coefficients exch
456
457	call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
458
459	! compute source and sink terms in the TKE equation (shear, buoyancy and dissipation)
460
461	call tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,dls,td,sh,bu,b_e,a_e,sf,ufrac_int)
462
463	! solve for tke
464
465	call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
466
467	! avoid negative values for tke
468
469	do iz=kts,kte
470	if(te(iz).lt.temin) te(iz)=temin
471	enddo
472
473	return
474	end subroutine boulac1d
475	!
476	! ===6=8===============================================================72
477
478	! ===6=8===============================================================72
479	subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
480	! compute the length scales up and down
481	implicit none
482	integer kms,kme,kts,kte,iz,izz,ix,iy
483	real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g
484	real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme)
485	real z(kms:kme),th(kms:kme),th0(kms:kme)
486
487	do iz=kts,kte
488	zup=0.
489	dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2.
490	zzz=0.
491	zup_inf=0.
492	beta=g/th0(iz) !Buoyancy coefficient
493	do izz=iz,kte-1
494	dzt=(dz(izz+1)+dz(izz))/2.
495	zup=zup-betath(iz)dzt
496	zup=zup+beta(th(izz+1)+th(izz))dzt/2.
497	zzz=zzz+dzt
498	if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then
499	bbb=(th(izz+1)-th(izz))/dzt
500	if(bbb.ne.0)then
501	tl=(-beta(th(izz)-th(iz))+sqrt( max(0.,(beta(th(izz)-th(iz)))*2.+2.bbbbeta(te(iz)-zup_inf))))/bbb/beta
502	else
503	if(th(izz).ne.th(iz))then
504	tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz)))
505	else
506	tl=0.
507	endif
508	endif
509	dlu(iz)=zzz-dzt+tl
510	endif
511	zup_inf=zup
512	enddo
513
514	zdo=0.
515	zdo_sup=0.
516	dld(iz)=z(iz)+dz(iz)/2.
517	zzz=0.
518	do izz=iz,kts+1,-1
519	dzt=(dz(izz-1)+dz(izz))/2.
520	zdo=zdo+betath(iz)dzt
521	zdo=zdo-beta(th(izz-1)+th(izz))dzt/2.
522	zzz=zzz+dzt
523	if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then
524	bbb=(th(izz)-th(izz-1))/dzt
525	if(bbb.ne.0.)then
526	tl=(beta(th(izz)-th(iz))+sqrt( max(0.,(beta(th(izz)-th(iz)))*2.+2.bbbbeta(te(iz)-zdo_sup))))/bbb/beta
527	else
528	if(th(izz).ne.th(iz))then
529	tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz)))
530	else
531	tl=0.
532	endif
533	endif
534
535	dld(iz)=zzz-dzt+tl
536	endif
537	zdo_sup=zdo
538	enddo
539	enddo
540
541
542	end subroutine dissip_bougeault
543	!
544	! ===6=8===============================================================72
545	! ===6=8===============================================================72
546	subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
547	! compute the length scales for dissipation and turbulent coefficients
548	implicit none
549	integer kms,kme,kts,kte,iz,ix,iy
550	real dlu(kms:kme),dld(kms:kme),dl_u(kms:kme)
551	real dls(kms:kme),dlk(kms:kme),dlg(kms:kme)
552
553	do iz=kts,kte
554	dld(iz)=min(dld(iz),dlg(iz))
555	dls(iz)=sqrt(dlu(iz)*dld(iz))
556	dlk(iz)=min(dlu(iz),dld(iz))
557
558	if(dl_u(iz).gt.0.)then
559	dls(iz)=1./(1./dls(iz)+1./dl_u(iz))
560	dlk(iz)=1./(1./dlk(iz)+1./dl_u(iz))
561	endif
562	enddo
563
564	return
565	end subroutine length_bougeault
566	!
567
568	! ===6=8===============================================================72
569	! ===6=8===============================================================72
570
571	subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
572	! compute turbulent coefficients
573	implicit none
574	integer iz,kms,kme,kts,kte,ix,iy
575	real te_m,dlk_m
576	real te(kms:kme),exch(kms:kme)
577	real dz(kms:kme),z(kms:kme)
578	real dlk(kms:kme)
579	real fact
580
581	exch(kts)=0.
582
583	! do iz=2,nz-1
584	do iz=kts+1,kte
585	te_m=(te(iz-1)dz(iz)+te(iz)dz(iz-1))/(dz(iz)+dz(iz-1))
586	dlk_m=(dlk(iz-1)dz(iz)+dlk(iz)dz(iz-1))/(dz(iz)+dz(iz-1))
587	exch(iz)=ck_bdlk_msqrt(te_m)
588	! exch(iz)=max(exch(iz),0.0001)
589	exch(iz)=max(exch(iz),0.1)
590	enddo
591
592	exch(kte+1)=0.1
593
594	return
595	end subroutine cdtur_bougeault
596
597
598	! ===6=8===============================================================72
599	! ===6=8===============================================================72
600
601	subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)
602
603
604	!------------------------------------------------------------------------
605	! Calculation of the diffusion in 1D
606	!------------------------------------------------------------------------
607	! - Input:
608	! nz : number of points
609	! iz1 : first calculated point
610	! co : concentration of the variable of interest
611	! dz : vertical levels
612	! cd : diffusion coefficients
613	! dtext : external time step
614	! rho : density of the air at the center
615	! rhoz : density of the air at the face
616	! itest : if itest eq 1 then update co, else store in a flux array
617	! - Output:
618	! co :concentration of the variable of interest
619
620	! - Internal:
621	! cddz : constant terms in the equations
622	! dt : diffusion time step
623	! nt : number of the diffusion time steps
624	! cstab : ratio of the stability condition for the time step
625	!---------------------------------------------------------------------
626
627	implicit none
628	integer iz,iz1,izf
629	integer kms,kme,kts,kte
630	real dt,dzv
631	real co(kms:kme),cd(kms:kme),dz(kms:kme)
632	real rho(kms:kme),rhoz(kms:kme)
633	real cddz(kms:kme+1),fc(kms:kme),df(kms:kme)
634	real a(kms:kme,3),c(kms:kme)
635	real sf(kms:kme),vl(kms:kme)
636	real aa(kms:kme),bb(kms:kme)
637
638
639	! Compute cddz=2*cd/dz
640
641	cddz(kts)=sf(kts)rhoz(kts)cd(kts)/dz(kts)
642	do iz=kts+1,kte
643	cddz(iz)=2.sf(iz)rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
644	enddo
645	cddz(kte+1)=sf(kte+1)rhoz(kte+1)cd(kte+1)/dz(kte)
646
647	do iz=kts,iz1-1
648	a(iz,1)=0.
649	a(iz,2)=1.
650	a(iz,3)=0.
651	c(iz)=co(iz)
652	enddo
653
654	do iz=iz1,kte-izf
655	dzv=vl(iz)*dz(iz)
656	a(iz,1)=-cddz(iz)*dt/dzv/rho(iz)
657	a(iz,2)=1+dt(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)dt
658	a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz)
659	c(iz)=co(iz)+bb(iz)*dt
660	enddo
661
662	do iz=kte-(izf-1),kte
663	a(iz,1)=0.
664	a(iz,2)=1
665	a(iz,3)=0.
666	c(iz)=co(iz)
667	enddo
668
669	call invert (kms,kme,kts,kte,a,c,co)
670
671	do iz=kts,iz1
672	fc(iz)=0.
673	enddo
674
675	do iz=iz1+1,kte
676	fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz)
677	enddo
678
679	! do iz=1,iz1
680	! df(iz)=0.
681	! enddo
682	!
683	! do iz=iz1+1,nz-izf
684	! dzv=vl(iz)*dz(iz)
685	! df(iz)=+(co(iz-1)cddz(iz)-co(iz)(cddz(iz)+cddz(iz+1))+co(iz+1)*cddz(iz+1))/dzv/rho(iz)
686	! enddo
687	!
688	! do iz=nz-izf,nz
689	! df(iz)=0.
690	! enddo
691
692	return
693	end subroutine diff
694
695	! ===6=8===============================================================72
696	! ===6=8===============================================================72
697
698	subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
699	! compute buoyancy term
700	implicit none
701	integer kms,kme,kts,kte,iz,ix,iy
702	real dtdz1,dtdz2,cdm,dtmdz,g
703	real th(kms:kme),exch(kms:kme),dz(kms:kme),bu(kms:kme)
704	real th0(kms:kme),ustar,tstar,ufrac_int
705
706	! bu(1)=-ustartstarg/th0(1)*(1.-ufrac_int)
707	bu(kts)=0.
708
709
710	do iz=kts+1,kte-1
711	dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
712	dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))
713	dtmdz=0.5*(dtdz1+dtdz2)
714	cdm=0.5*(exch(iz+1)+exch(iz))
715	bu(iz)=-cdmdtmdzg/th0(iz)
716	enddo
717	!
718
719	bu(kte)=0.
720
721	return
722	end subroutine buoy
723
724	! ===6=8===============================================================72
725	! ===6=8===============================================================72
726
727	subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int)
728	! compute shear term
729	implicit none
730	integer kms,kme,kts,kte,iz,ix,iy
731	real dudz1,dudz2,dvdz1,dvdz2,cdm,dumdz,ustar
732	real tstar,th,al,phim,g
733	real u(kms:kme),v(kms:kme),cdua(kms:kme),dz(kms:kme),sh(kms:kme)
734	real u1,u2,v1,v2,ufrac_int
735
736	! al=vkgtstar/(th(ustar*2.))
737	! if(al.ge.0.)phim=1.+4.7dz(1)/2.al
738	! if(al.lt.0.)phim=(1.-15dz(1)/2.al)**(-0.25)
739	!
740	! sh(1)=(ustar*3.)/vk/(dz(1)/2.)(1.-ufrac_int)
741	sh(kts)=0.
742	do iz=kts+1,kte-1
743	u2=(dz(iz+1)u(iz)+dz(iz)u(iz+1))/(dz(iz)+dz(iz+1))
744	u1=(dz(iz)u(iz-1)+dz(iz-1)u(iz))/(dz(iz-1)+dz(iz))
745	v2=(dz(iz+1)v(iz)+dz(iz)v(iz+1))/(dz(iz)+dz(iz+1))
746	v1=(dz(iz)v(iz-1)+dz(iz-1)v(iz))/(dz(iz-1)+dz(iz))
747	cdm=0.5*(cdua(iz)+cdua(iz+1))
748	dumdz=((u2-u1)/dz(iz))2.+((v2-v1)/dz(iz))2.
749	sh(iz)=cdm*dumdz
750	enddo
751
752	!!!!!!!
753	sh(kte)=0.
754
755	return
756	end subroutine shear
757
758	! ===6=8===============================================================72
759	! ===6=8===============================================================72
760
761	subroutine invert(kms,kme,kts,kte,a,c,x)
762
763	!ccccccccccccccccccccccccccccccc
764	! Aim: Inversion and resolution of a tridiagonal matrix
765	! A X = C
766	! Input:
767	! a(*,1) lower diagonal (Ai,i-1)
768	! a(*,2) principal diagonal (Ai,i)
769	! a(*,3) upper diagonal (Ai,i+1)
770	! c
771	! Output
772	! x results
773	!ccccccccccccccccccccccccccccccc
774
775	implicit none
776	integer in
777	integer kts,kte,kms,kme
778	real a(kms:kme,3),c(kms:kme),x(kms:kme)
779
780	do in=kte-1,kts,-1
781	c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
782	a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
783	enddo
784
785	do in=kts+1,kte
786	c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
787	enddo
788
789	do in=kts,kte
790
791	x(in)=c(in)/a(in,2)
792
793	enddo
794
795	return
796	end subroutine invert
797
798	! ===6=8===============================================================72
799	! ===6=8===============================================================72
800
801	subroutine tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch, &
802	dls,td,sh,bu,b_e,a_e,sf,ufrac_int)
803	! in this routine the shear, buoyancy and part of the dissipation terms
804	! of the TKE equation are computed
805
806	implicit none
807	integer kms,kme,kts,kte,iz,ix,iy
808	real g,ustar,tstar,ufrac_int
809	real z(kms:kme),dz(kms:kme),u(kms:kme),v(kms:kme),th(kms:kme),th0(kms:kme),te(kms:kme)
810	real exch(kms:kme),dls(kms:kme),td(kms:kme),sh(kms:kme),bu(kms:kme)
811	real a_e(kms:kme),b_e(kms:kme)
812	real vl(kms:kme),sf(kms:kme)
813	real te1,dl1
814
815	call shear(ix,iy,g,kms,kme,kts,kte,u,v,exch,dz,sh,ustar,tstar,th(kts),ufrac_int)
816
817	call buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
818
819	do iz=kts,kte
820	te1=max(te(iz),temin)
821	dl1=max(dls(iz),0.1)
822	td(iz)=-ceps_b*sqrt(te1)/dl1
823	sh(iz)=sh(iz)*sf(iz)
824	bu(iz)=bu(iz)*sf(iz)
825	a_e(iz)=a_e(iz)+td(iz)
826	b_e(iz)=b_e(iz)+sh(iz)+bu(iz)
827	enddo
828
829
830	return
831	end subroutine tke_bougeault
832	!######################################################################
833	subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh)
834
835	! this routine computes the PBL height
836	! with an approach similar to MYNN
837	implicit none
838	integer kms,kme,kts,kte,iz
839	real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
840	real pblh
841	!Local
842	real thv(kms:kme),zc(kms:kme)
843	real thsfc
844
845	! compute the height of the center of the grid cells
846	do iz=kts,kte
847	zc(iz)=z(iz)+dz(iz)/2.
848	enddo
849
850	! compute the virtual potential temperature
851
852	do iz=kts,kte
853	thv(iz)=th(iz)(1.+0.61q(iz))
854	enddo
855	! now compute the PBL height
856
857	pblh=0.
858	thsfc=thv(kts)+0.5
859	do iz=kts+1,kte
860	if(pblh.eq.0.and.thv(iz).gt.thsfc)then
861	pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1))
862	! pblh=z(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(z(iz)-z(iz-1))
863	endif
864	enddo
865
866	return
867	end subroutine pbl_height
868
869
870	! ===6=8===============================================================72
871
872
873	! ===6=8===============================================================72
874	SUBROUTINE BOULACINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
875	& TKE_PBL,EXCH_H,RESTART,ALLOWED_TO_READ, &
876	& IDS,IDE,JDS,JDE,KDS,KDE, &
877	& IMS,IME,JMS,JME,KMS,KME, &
878	& ITS,ITE,JTS,JTE,KTS,KTE )
879	!-----------------------------------------------------------------------
880	IMPLICIT NONE
881	!-----------------------------------------------------------------------
882	LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
883	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
884	& IMS,IME,JMS,JME,KMS,KME, &
885	& ITS,ITE,JTS,JTE,KTS,KTE
886
887	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
888	& RUBLTEN, &
889	& RVBLTEN, &
890	& RTHBLTEN, &
891	& RQVBLTEN, &
892	& TKE_PBL
893	INTEGER :: I,J,K,ITF,JTF,KTF
894	!-----------------------------------------------------------------------
895	!-----------------------------------------------------------------------
896
897	JTF=MIN0(JTE,JDE-1)
898	KTF=MIN0(KTE,KDE-1)
899	ITF=MIN0(ITE,IDE-1)
900
901	IF(.NOT.RESTART)THEN
902	DO J=JTS,JTF
903	DO K=KTS,KTF
904	DO I=ITS,ITF
905	TKE_PBL(I,K,J)=0.0001
906	RUBLTEN(I,K,J)=0.
907	RVBLTEN(I,K,J)=0.
908	RTHBLTEN(I,K,J)=0.
909	RQVBLTEN(I,K,J)=0.
910	EXCH_H(I,K,J)=0.
911	ENDDO
912	ENDDO
913	ENDDO
914	ENDIF
915
916	END SUBROUTINE BOULACINIT
917
918
919
920	END MODULE module_bl_boulac
921

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

Download in other formats:

Original Format