Context Navigation

gettherm.pro @ 205

Last change on this file since 205 was 205, checked in by acolaitis, 13 years ago
Small commit to add A.C. thermals diagnostics IDL routines. gettherm is for extracting data and saving it in thermiques.dat, and thermiques is the monstruous 4000 lines plotting routine.
File size: 118.2 KB

Line
1	PRO gettherm
2
3	spawn, 'clear'
4	print, ''
5	print, ' Thermals Analysis '
6	print, ' (usine à gaz) '
7	print, ''
8
9	full='true'
10	f_offset='false'
11	overplot_convadj='false'
12	plot_3d = 'false'
13
14	lctu_gcm = 8.
15
16	; *********** Best values for LES thermals with tau =0.5
17	betalpha = 1.3
18	afact = 1.8
19	fact_epsilon = 0.0008
20	detr_min = 0.0007
21
22
23	;betalpha = 1.
24	;afact = 2.4
25	;fact_epsilon = 0.0007
26	;detr_min = 0.0007
27
28	; --------
29
30	; *********** Best values for LES thermals with tau =1.
31	;betalpha = 1.
32	;afact = 1.1
33	;fact_epsilon = 0.00025
34	;detr_min = 0.0007
35	; --------
36
37	;datname='thermiques.dat.scale1.2'
38	;datname='thermiques.dat.scale1.4'
39	;datname='thermiques.dat.scale0.6'
40	datname='thermiques.dat'
41	;datname='thermiques.dat.scale0.8'
42	;datname='thermiques.dat' ;scale =1.0, sigmao =1.0
43
44	ns = 0 ; number of points for time-smoothing of LES data : 2*ns+1 points, ns = 9 eq to 30mn (-15mn//+15mn)
45	nstot = float(2.*ns+1.)
46
47	GcmSubCase = ''
48	LayerCase=''
49	s_trac1 = 'qtrac1'
50	s_trac2 = 'qtrac2'
51	got_pdt = 'true'
52	TestCase = 'Case_A'
53	SubCase = '_11_shorter'
54	Histo = 'true'
55	newtest = ''
56	visualization_mode = 'false'
57	label_init:
58	spawn, 'clear'
59	print, ' Available simulations :'
60	print, '----'
61	print, ' 1/ Case_A_11 : 45x45x71, ztop=10km,dx=100m,dz=140m,Ls=47.1°,(21.8N;205.0E),55tiu,A=0.275,Tau=0.5'
62	print, ' 2/ Case_A_4 : 101x101x201 ztop=15km,dx=100m,dz= 75m,Ls=47.1°,(21.8N;205.0E),55tiu,A=0.275,Tau=0.5'
63	print, ' 3/ Case_A_4_shorter : Case A4 with Dtrac1 = 5 mn and Dtrac2 = 10 mn (compared to 20 and 100)'
64	print, ' 4/ Case_A_4_shorter_winds : Case A4_shorter with bckgrnd wind u=10 m/s'
65	print, ' 5/ Case_A_4_shorter_winds_30 : Case A4_shorter with bckgrnd wind u=30 m/s'
66	print, ' 6/ Case_A_4_shorter_winds_tau1 : Case A4_shorter with bckgrnd wind u=10 m/s and tau=1'
67	print, ' 7/ Case_A_4_shorter_winds_tau2 : Case A4_shorter with bckgrnd wind u=10 m/s and tau=2'
68	print, ' 8/ Case_ExtremeCase : 101x101x201 ztop=15km,dx=100m,dz=75m,Ls=0°,(0.N;0.E),50tiu,A=0.1,Tau=0.05'
69	print, '----'
70	print, ' 9/ Case_C_4_shorter_winds : '
71	print, ' 10/ Case_I_4_shorter_winds : '
72	print, ' 11/ Case_Z_4_shorter_winds : '
73	print, ' 12/ 1D : 124 layers'
74	print, ' 13/ 1D : 32 layers'
75	print, ' 14/ 1D : low dt'
76	print, ''
77	print, ' 0/ PLUME VISUALISATION : '+visualization_mode
78	print, ' 999/ CLEAR thermiques.dat for considered case'
79	print, ' 100/ OVERPLOT CONVADJ ONLY RESULTS : '+overplot_convadj
80	print, ''
81	print, ' SIMULATION NUMBER : '
82	print, ' '+TestCase+SubCase+LayerCase+' '
83	print, ''
84	print, 'Any change ? (number of new case to change, or any other key to continue)'
85	read, newtest
86	if (newtest eq '1') then begin
87	TestCase = 'Case_A'
88	SubCase = '_11'
89	pGround = 867.5594
90	goto,label_init
91	endif
92	if (newtest eq '2') then begin
93	TestCase = 'Case_A'
94	SubCase = '_4'
95	got_pdt = 'false'
96	pGround = 867.5594
97	goto,label_init
98	endif
99	if (newtest eq '3') then begin
100	TestCase = 'Case_A'
101	SubCase = '_4_shorter'
102	;s_trac1 = 'qtrac2'
103	;s_trac2 = 'qtrac1'
104	full = 'true'
105	f_offset='false'
106	pGround = 867.5594
107	goto,label_init
108	endif
109	if (newtest eq '4') then begin
110	TestCase = 'Case_A'
111	SubCase = '_4_shorter_wind'
112	;s_trac1 = 'qtrac2'
113	;s_trac2 = 'qtrac1'
114	full = 'true'
115	f_offset='false'
116	pGround = 867.5594
117	goto,label_init
118	endif
119	if (newtest eq '5') then begin
120	TestCase = 'Case_A'
121	SubCase = '_4_shorter_wind_30'
122	GcmSubCase = '_wind_30'
123	;s_trac1 = 'qtrac2'
124	;s_trac2 = 'qtrac1'
125	full = 'true'
126	pGround = 867.5594
127	f_offset='false'
128	goto,label_init
129	endif
130	if (newtest eq '6') then begin
131	TestCase = 'Case_A'
132	SubCase = '_4_shorter_wind_tau1'
133	GcmSubCase = '_tau1'
134	;s_trac1 = 'qtrac2'
135	;s_trac2 = 'qtrac1'
136	full = 'true'
137	f_offset='false'
138	pGround = 867.5594
139	goto,label_init
140	endif
141	if (newtest eq '7') then begin
142	TestCase = 'Case_A'
143	SubCase = '_4_shorter_wind_tau2'
144	GcmSubCase = '_tau2'
145	;s_trac1 = 'qtrac2'
146	;s_trac2 = 'qtrac1'
147	full = 'true'
148	f_offset='false'
149	pGround = 867.5594
150	goto,label_init
151	endif
152	if (newtest eq '8') then begin
153	TestCase = 'ExtremeCase'
154	SubCase = ''
155	GcmSubCase = ''
156	;s_trac1 = 'qtrac2'
157	;s_trac2 = 'qtrac1'
158	full = 'true'
159	f_offset='false'
160	pGround = 677.722
161	;lctu_gcm = 6.
162	goto,label_init
163	endif
164	if (newtest eq '9') then begin
165	TestCase = 'Case_C'
166	SubCase = '_4_shorter_wind'
167	;s_trac1 = 'qtrac2'
168	;s_trac2 = 'qtrac1'
169	full = 'true'
170	f_offset='false'
171	pGround = 483.
172	goto,label_init
173	endif
174	if (newtest eq '10') then begin
175	TestCase = 'Case_I'
176	SubCase = '_4_shorter_wind'
177	;s_trac1 = 'qtrac2'
178	;s_trac2 = 'qtrac1'
179	full = 'true'
180	f_offset='false'
181	pGround = 630.
182	goto,label_init
183	endif
184	if (newtest eq '11') then begin
185	TestCase = 'Case_Z'
186	SubCase = '_4_shorter_wind'
187	;s_trac1 = 'qtrac2'
188	;s_trac2 = 'qtrac1'
189	full = 'true'
190	f_offset='false'
191	pGround = 266.
192	goto,label_init
193	endif
194	if (newtest eq '12') then begin
195	LayerCase=''
196	goto,label_init
197	endif
198	if (newtest eq '13') then begin
199	LayerCase='_32lev'
200	goto,label_init
201	endif
202	if (newtest eq '14') then begin
203	LayerCase='_low_dt'
204	goto,label_init
205	endif
206	if (newtest eq '100') then begin
207	if (overplot_convadj eq 'true') then overplot_convadj = 'false' else overplot_convadj = 'true'
208	goto,label_init
209	endif
210	if (newtest eq '0') then begin
211	visualization_mode = 'true'
212	spawn, 'clear'
213	print, 'The first timestep of the considered file will be shown.'
214	print, 'Defaut file is : file 6 (lt ~ 13h10)'
215	print, 'Which file do you want ? (lt ~= file_number + 7)'
216	loop_special = '6'
217	read, loop_special
218	print, 'Do you wish to plot histograms or manipulate volumic data ?'
219	print, '1/ Histogram'
220	print, '2/ Volumic data'
221	x=''
222	read, x
223	if (x eq '1') then Histo='true' else Histo='false'
224	goto,label_init
225	endif
226
227	;les_path='/san0/acolmd/SIMUS/LES_'+TestCase+SubCase
228	;gcm_path='/san0/acolmd/SIMUS/GCM_'+TestCase+LayerCase+'_2'
229	;gcm_path='/san0/acolmd/SIMUS/GCM_'+TestCase+GcmSubCase+LayerCase
230	;gcm_convadj_path=gcm_path+'_convadj'
231
232	;les_path='/data/acolmd/Thermiques/LES_Case_A_4_2trac_wind10_tau05/'
233	;les_path='/data/acolmd/Thermiques/LES_Case_A_257x257x301_wind10_tau05'
234	;les_path='/data/acolmd/Thermiques/LES_Case_E_257x257x301_wind10_tau005'
235	;les_path='/data/acolmd/Thermiques/LES_Case_A_4_2trac_wind30_tau05/'
236	;les_path='/data/acolmd/Thermiques/LES_Case_A_4_2trac_wind10_tau1/'
237	;les_path='/data/acolmd/Thermiques/LES_Case_A_4_2trac_wind10_tau2/'
238	;les_path='/data/acolmd/Thermiques/LES_Case_C_4_2trac_wind10_tau05/'
239	;les_path='/data/acolmd/Thermiques/LES_Case_I_4_2trac_wind10_tau05/'
240	;les_path='/data/acolmd/Thermiques/LES_Case_Z_4_2trac_wind10_tau05/'
241	;les_path='/data/acolmd/Thermiques/ExtremeCase/'
242
243	;les_path='/data/acolmd/Thermiques/LES_Case_A_101x101x201_tracup_wind10_tau05_gcmsoil/'
244	;les_path='/data/acolmd/Thermiques/LES_Case_E_101x101x201_tracup_wind10_tau005_gcmsoil/'
245	;les_path='/data/acolmd/Thermiques/LES_Case_C_101x101x201_tracup_wind10_tau05_gcmsoil'
246	;les_path='/data/acolmd/Thermiques/LES_Case_I_101x101x201_tracup_wind10_tau05_gcmsoil'
247	;les_path='/data/acolmd/Thermiques/LES_Case_Z_101x101x201_tracup_wind10_tau05_gcmsoil'
248	;les_path='/data/acolmd/Thermiques/LES_Case_A_101x101x201_tracup_wind30_tau05_gcmsoil'
249	les_path='/data/acolmd/Thermiques/LES_Case_A_101x101x201_tracup_wind10_tau2_gcmsoil'
250
251	gcm_path='/data/acolmd/Thermiques/THgcm/'
252	gcm_convadj_path=gcm_path+'_convadj'
253
254
255	if (newtest eq '999') then spawn, 'rm -f '+les_path+'/'+datname
256	print, ''
257	print, ' -- Loading LES data -- '
258
259	print, 'LES DATA IN : '
260	print, les_path
261	print, 'GCM DATA IN : '
262	print, gcm_path
263
264	p0=610. & t0=220. & r_cp=1./3.89419 & grav=3.72 & R=191.182
265	history_interval_s = 100.
266	;lctu_gcm = 8. ; Initial local time of gcm 1d simu
267	scale = 1. ; Scaling factor for conditional sampling
268	decimate = 10. ; Coeff for subsampling the data for sigma integral
269	sigmao= 0.3 ; multiplicative coeff for the computation of Sigma0 in the CS
270	sigmao_ude = 0.2 ; number of standard deviation away from mean for the selection of downdraft in UDE
271	NBINS=100.
272
273	openr,unit,les_path+'/'+datname,/get_lun,error=err
274	IF (err ne 0) THEN BEGIN
275
276	OPENR, 22, les_path+'/input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22
277	OPENR, 23, les_path+'/input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
278
279	domain='d01' & filesWRF = FindFile(les_path+'/wrfout_'+domain+'_????-??-??_??:??:??') & nf=n_elements(filesWRF)
280	;print, filesWRF
281	;domain='d01' & filesWRF = les_path+'/'+['wrfout_d01_9999-01-01_03:05:00','wrfout_d01_9999-01-01_04:06:40'] & nf=n_elements(filesWRF)
282	;ce fichier utilise aussi offset_localtime = 3.05
283
284	; WARNING WARNING : FOR THE CASE 4_SHORTER, THE TRAC2 HAS 10 MN LIFETIME, WE WANT TO USE IT MORE EXTENSIVELY THAN THE TRAC1 (5mn) SO
285	; we switch the names of trac1 and trac2 in the initialization of this routine, in the "case".
286
287	; WARNING WARNING : in this version (thermiques2), one of the tracers in the data is just a neutral tracer which will be used to compute the turbulent flux of tracer MMR. This is a priori done by the tracer 2 (after inversion)
288
289	id=ncdf_open(filesWRF(0))
290	NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east' ), dummy, nx & NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north' ), dummy, ny
291	NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top' ), dummy, nz & NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), dummy, nt
292	NCDF_CLOSE, id
293	id=ncdf_open(filesWRF(nf-1)) ;; for interrupted runs
294	NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), dummy, ntlast
295	NCDF_CLOSE, id
296	nttot = (nf-1)*nt + ntlast
297	wt = fltarr(nz,nttot) & wq = fltarr(nz,nttot) & wq_updraft = fltarr(nz,nttot) & wq_downdraft = fltarr(nz,nttot) & wq_env_ude=fltarr(nz,nttot)
298	q_mean_up = fltarr(nz,nttot) & q_mean_down = fltarr(nz,nttot) & q_mean_env_ude = fltarr(nz,nttot) & q_mean = fltarr(nz,nttot)
299	tke_les = fltarr(nz,nttot) & ztke = fltarr(nz,nttot) & t = fltarr(nz,nttot)
300	p = fltarr(nz) & ph = fltarr(nz) & pht = fltarr(nz,nttot) & pt = fltarr(nz,nttot)
301	xtke = fltarr(nz,nttot) & ytke = fltarr(nz,nttot) & temp_les = fltarr(nz,nttot)
302	wmax = fltarr(nttot)
303	alpha1 = fltarr(nz) & alpha2 = fltarr(nz)
304	alpha1out = fltarr(nz,nttot) & alpha2out = fltarr(nz,nttot)
305	zqtrac1 = dblarr(nx,ny,nz) & zqtrac2 = dblarr(nx,ny,nz)
306	sigmazqtrac1 = fltarr(nz) & sigmazqtrac2 = fltarr(nz)
307	sigmazminqtrac1 = fltarr(nz) & sigmazminqtrac2 = fltarr(nz)
308	fm_trac1_les = fltarr(nz,nttot) & fm_trac2_les = fltarr(nz,nttot)
309	anomalqtrac1 = fltarr(nx,ny,nz) & anomalqtrac2 = fltarr(nx,ny,nz)
310	e_trac1_les = fltarr(nz,nttot) & e_trac2_les = fltarr(nz,nttot)
311	dtempdztmp = fltarr(nx,ny,nz)
312	localtime = lctu + history_interval_s*findgen(nttot)/3700.
313	w_mean1 = fltarr(nz,nttot) & w_mean2 = fltarr(nz,nttot)
314	w_mean1_env = fltarr(nz,nttot) & w_mean1_down = fltarr(nz,nttot)
315	w_mean1_env_ude = fltarr(nz,nttot) & w_mean1_full = fltarr(nz,nttot)
316	buoyancy1_les = fltarr(nz,nttot) & buoyancy2_les = fltarr(nz,nttot)
317	e1_term2 = fltarr(nz,nttot) & e1_term3 = fltarr(nz,nttot)
318	dtetadttmp = fltarr(nx,ny,nz)
319	rhomoy1 = fltarr(nz,nttot)
320	plumeIndex1out = make_array(nxny,nz,VALUE=-1.) & envIndex1out = make_array(nxny,nz,VALUE=-1.)
321	hf1tmp = fltarr(nz,nttot) & hf1tmpenv = fltarr(nz,nttot)
322	tplume1moy = fltarr(nz,nttot) & tenv1moy = fltarr(nz,nttot) & tenv1moy_ude = fltarr(nz,nttot)
323	tmoy_full = fltarr(nz,nttot) & tdown1moy = fltarr(nz,nttot)
324	dteta1moydt_entr = fltarr(nz,nttot) & dteta1moydt_detr = fltarr(nz,nttot)
325	d1_term1 = fltarr(nz,nttot) & d1_term2 = fltarr(nz,nttot) & d1_term3=fltarr(nz,nttot)
326	hf1_term1 = fltarr(nz,nttot) & hf1_term2 = fltarr(nz,nttot) & hf1_term3 = fltarr(nz,nttot)
327	d1_term1_ude = fltarr(nz,nttot) & d1_term2_ude = fltarr(nz,nttot) & d1_term3_ude=fltarr(nz,nttot)
328	e1_term1_ude = fltarr(nz,nttot) & e1_term2_ude = fltarr(nz,nttot) & e1_term3_ude=fltarr(nz,nttot)
329	downward_flux1 = fltarr(nz,nttot) & beta1out = fltarr(nz,nttot)
330	hf1_ude_term1 = fltarr(nz,nttot) & hf1_ude_term2 = fltarr(nz,nttot) & hf1_ude_term3 = fltarr(nz,nttot) & hf1_ude_term4 = fltarr(nz,nttot)
331	hf1tmpenv_ude = fltarr(nz,nttot) & hf1tmp_down = fltarr(nz,nttot)
332	dTeta_phys = make_array(nz,nttot)
333	exner = fltarr(nz,nttot)
334	uv_moy = fltarr(nz,nttot)
335	tsurf = fltarr(nttot)
336	Gamma_1 = fltarr(nz,nttot) & Gamma_2 = fltarr(nz,nttot) & Gamma_3 = fltarr(nz,nttot)
337	Gamma_1_tmp = fltarr(nz,nttot) & dgamma1tmp = fltarr(nz,nttot)
338	ptotprime = fltarr(nx,ny,nz) & anomalptot = fltarr(nx,ny,nz) & dptotprimedztmp = fltarr(nx,ny,nz)
339	hfx = fltarr(nttot) & flxrad = fltarr(nttot) & flxgrd = fltarr(nttot) & lwdownz = fltarr(nttot) & swdownz = fltarr(nttot)
340
341	;wBin = fltarr(NBINS,nttot) & wBinEnv_ude = fltarr(NBINS,nttot) & wBinUp = fltarr(NBINS,nttot) & wBinDown = fltarr(NBINS,nttot)
342	l=0
343
344	FOR loop = 0, nf-1 DO BEGIN
345	timetime = SYSTIME(1)
346	if (loop ne nf-1) then nloop2=nt else nloop2=ntlast
347	if (loop ne 0) then loop2_init=0 else loop2_init=1 ;le tout premier pas de temps est l'initialisation, certains champs sont à 0 => bug
348	FOR loop2 = loop2_init, nloop2-1 DO BEGIN
349	pht(*,l) = TOTAL(TOTAL(getget(filesWRF(loop), 'PHTOT', count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny) / 1000. / 3.72
350	ph = TEMPORARY(ph) + pht(*,l) / (nttot-1)
351	ENDFOR
352	print, 'computing altitudes, file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s'
353	ENDFOR
354	altitudes_LES = 1000.*(TEMPORARY(ph) - hgtu/1000.) ;; altitude above ground
355	pht = fltarr(nz,nttot)
356	ph = fltarr(nz)
357
358	FOR loop = 0, nf-1 DO BEGIN
359	timetime = SYSTIME(1)
360	if (loop ne nf-1) then nloop2=nt else nloop2=ntlast
361	if (loop ne 0) then loop2_init=0 else loop2_init=1 ;le tout premier pas de temps est l'initialisation, certains champs sont à 0 => bug
362	FOR loop2 = loop2_init, nloop2-1 DO BEGIN
363
364	anomalt = 1. & anomalu = 1. & anomalv = 1. & anomalw = 1.
365	; --------------------------------------------------------
366	; u' = u and v' = v (car PAS de background wind !)
367	; tke = 0.5 ( <u'^2> + <v'^2> + <w'^2> ) ; u' = u ; v' = v
368	; --------------------------------------------------------
369
370	tprime = getget(filesWRF(loop), 'T', anomaly=anomalt, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; t' = t - <t>
371	t(*,l) = t0 + temporary(anomalt)
372	ztke(,l) = 0.5 TOTAL(TOTAL(getget(filesWRF(loop), 'W', anomaly=anomalw, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny)
373	xtke(,l) = 0.5 TOTAL(TOTAL(getget(filesWRF(loop), 'U', anomaly=anomalu, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny)
374	ytke(,l) = 0.5 TOTAL(TOTAL(getget(filesWRF(loop), 'V', anomaly=anomalv, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny)
375	uv_moy(*,l) = TOTAL(TOTAL(sqrt(getget(filesWRF(loop), 'U', count=[0,0,0,1], offset=[0,0,0,loop2])^2 + getget(filesWRF(loop), 'V', count=[0,0,0,1], offset=[0,0,0,loop2])^2),1),1) / float(nx) / float(ny)
376	tke_les(,l) = xtke(,l) + ytke(,l) + ztke(,l)
377	wprime = getget(filesWRF(loop), 'W', anomaly=anomalw, count=[0,0,0,1], offset=[0,0,0,loop2])
378	pht(*,l) = TOTAL(TOTAL(getget(filesWRF(loop), 'PHTOT', count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny) / 1000. / 3.72
379	pt(*,l) = TOTAL(TOTAL(getget(filesWRF(loop), 'PTOT' , count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny)
380	tsurf(l) = TOTAL(TOTAL(getget(filesWRF(loop), 'TSURF' , count=[0,0,1], offset=[0,0,loop2]),1),1) / float(nx) / float(ny)
381	hfx(l) = TOTAL(TOTAL(getget(filesWRF(loop), 'HFX' , count=[0,0,1], offset=[0,0,loop2]),1),1) / float(nx) / float(ny)
382	flxrad(l) = TOTAL(TOTAL(getget(filesWRF(loop), 'FLXRAD' , count=[0,0,1], offset=[0,0,loop2]),1),1) / float(nx) / float(ny)
383	flxgrd(l) = TOTAL(TOTAL(getget(filesWRF(loop), 'FLXGRD' , count=[0,0,1], offset=[0,0,loop2]),1),1) / float(nx) / float(ny)
384	lwdownz(l) = TOTAL(TOTAL(getget(filesWRF(loop), 'LWDOWNZ' , count=[0,0,1], offset=[0,0,loop2]),1),1) / float(nx) / float(ny)
385	swdownz(l) = TOTAL(TOTAL(getget(filesWRF(loop), 'SWDOWNZ' , count=[0,0,1], offset=[0,0,loop2]),1),1) / float(nx) / float(ny)
386
387	temp_les(,l) = t(,l)(pt(,l)/p0)^r_cp
388	IF (got_pdt eq 'true') then begin
389	exner(,l) = (pt(,l)/p0)^r_cp
390	dTeta_phys(,l) = (TOTAL(TOTAL(getget(filesWRF(loop), 'PDT', count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny))/exner(,l)
391	ENDIF
392	ph = TEMPORARY(ph) + pht(*,l) / (nttot-1)
393	p = TEMPORARY(p ) + pt(*,l) / nttot
394	ptotprime = getget(filesWRF(loop), 'PTOT', count=[0,0,0,1], offset=[0,0,0,loop2])
395	FOR k=0, nz-1 DO BEGIN
396	rhomoy1(,l) = TOTAL(TOTAL(reform(((ptotprime(,,k)/(R(t(k,l)+tprime(,,k))))(p0/ptotprime(,,k))^r_cp)),1),1)/(float(nx)float(ny))
397	anomalptot(,,k) = ptotprime(,,k) - total(total(ptotprime(,,k),1),1)/ float(nx) / float(ny) ; prime par rapport a la moyenne du niveau (plus continu)
398	ENDFOR
399	zqtrac1 = getget(filesWRF(loop), s_trac1, count=[0,0,0,1], offset=[0,0,0,loop2])
400	; zqtrac2 = getget(filesWRF(loop), s_trac2, count=[0,0,0,1], offset=[0,0,0,loop2])
401	FOR i=0,nx-1 DO BEGIN
402	FOR j=0, ny-1 DO BEGIN
403	dtempdztmp(i,j,) = deriv(altitudes_LES, tprime(i,j,) + t(*,l))
404	dptotprimedztmp(i,j,) = deriv(altitudes_LES, anomalptot(i,j,))
405	ENDFOR
406	ENDFOR
407
408	FOR k=0,nz-1 DO BEGIN
409	anomalqtrac1(,,k) = zqtrac1(,,k) - TOTAL(TOTAL(REFORM(zqtrac1(,,k)),1),1)/ float(nx) / float(ny)
410	; anomalqtrac2(,,k) = zqtrac2(,,k) - TOTAL(TOTAL(REFORM(zqtrac2(,,k)),1),1)/ float(nx) / float(ny)
411	sigmazqtrac1(k) = STDDEV(REFORM(zqtrac1(,,k)))
412	IF (k ne 0) THEN BEGIN
413	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
414	sigmazminqtrac1(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac1(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
415	ENDIF ELSE BEGIN
416	sigmazminqtrac1(k) = sigmazqtrac1(k)
417	ENDELSE
418	plumeIndex1 = WHERE((anomalqtrac1(,,k) GT scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) AND ((anomalw(k)+wprime(,*,k)) GT 0.))
419	envIndex1 = WHERE((anomalqtrac1(,,k) LE scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) OR ((anomalw(k)+wprime(,*,k)) LE 0.))
420	; plumeIndex1 = WHERE(anomalqtrac1(,,k) GT scale*MAX([sigmazqtrac1(k),sigmazminqtrac1(k)]))
421	; envIndex1 = WHERE(anomalqtrac1(,,k) LE scale*MAX([sigmazqtrac1(k),sigmazminqtrac1(k)]))
422	IF(plumeIndex1(0) EQ -1) THEN BEGIN
423	fm_trac1_les(k,l)=0.
424	e_trac1_les(k,l)=0.
425	alpha1out(k,l)=0.
426	buoyancy1_les(k,l)=0.
427	w_mean1(k,l)=0.
428	w_mean1_env(k,l)=0.
429	w_mean1_down(k,l)=0.
430	w_mean1_full(k,l)=0.
431	w_mean1_env_ude(k,l)=0.
432	e1_term2(k,l)=0.
433	e1_term3(k,l)=0.
434	e1_term1_ude(k,l)=0.
435	e1_term2_ude(k,l)=0.
436	e1_term3_ude(k,l)=0.
437	hf1tmp(k,l)=0.
438	hf1tmpenv(k,l)=0.
439	Gamma_1_tmp(k,l)=0.
440	plumeIndex1out(*,k)=-1.
441	envIndex1out(*,k)=-1.
442	d1_term1(k,l)=0.
443	d1_term2(k,l)=0.
444	d1_term3(k,l)=0.
445	d1_term1_ude(k,l)=0.
446	d1_term2_ude(k,l)=0.
447	d1_term3_ude(k,l)=0.
448	downward_flux1(k,l)=0.
449	beta1out(k,l)=0.
450	tmoy_full(k,l)=0.
451	tdown1moy(k,l)=0.
452	wq_updraft(k,l)=0.
453	wq_downdraft(k,l)=0.
454	q_mean_up(k,l)=0.
455	q_mean_down(k,l)=0.
456	Gamma_2(k,l)=0.
457	Gamma_3(k,l)=0.
458	ENDIF ELSE BEGIN
459	FOR n=0L,n_elements(plumeIndex1)-1 DO BEGIN
460	plumeIndex1out(n,k)=plumeIndex1(n)
461	ENDFOR
462	FOR n=0L,n_elements(envIndex1)-1 DO BEGIN
463	envIndex1out(n,k)=envIndex1(n)
464	ENDFOR
465	alpha1(k) = n_elements(plumeIndex1) / float(nx) / float(ny)
466	wprimetmp = reform(reform((anomalw(k)+wprime(,,k))),[nx*ny,1])
467	w_mean1_full(k,l) = mean(wprimetmp)
468	w_mean1(k,l) = mean(wprimetmp(plumeIndex1))
469	w_mean1_env(k,l) = mean(wprimetmp(envIndex1))
470	downdraft_index1 = WHERE((abs(anomalw(k)+wprime(,,k)) gt sigmao_udeSTDDEV(wprimetmp(envIndex1))) and (anomalw(k)+wprime(,*,k) lt 0.))
471
472	envIndex1_ude = WHERE(((abs(anomalw(k)+wprime(,,k)) le sigmao_udeSTDDEV(wprimetmp(envIndex1))) or (anomalw(k)+wprime(,,k) ge 0.)) AND ((anomalqtrac1(,,k) LE scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) OR ((anomalw(k)+wprime(,,k)) LE 0.)))
473	IF (envIndex1_ude(0) ne -1) THEN w_mean1_env_ude(k,l) = mean(wprimetmp(envIndex1_ude)) ELSE w_mean1_env_ude(k,l) =0.
474	if (downdraft_index1(0) ne -1) then begin
475	w_mean1_down(k,l)=mean(wprimetmp(downdraft_index1))
476	wprimetmp=0.
477	beta1 = n_elements(downdraft_index1) / float(nx) / float(ny)
478	beta1out(k,l)=beta1
479	downward_flux1(k,l) = beta1rhomoy1(k,l)w_mean1_down(k,l)
480	endif else begin
481	downward_flux1(k,l)=0.
482	beta1out(k,l)=0.
483	w_mean1_down(k,l)=0.
484	tdown1moy(k,l)=0.
485	endelse
486	fm_trac1_les(k,l) = alpha1(k)rhomoy1(k,l)w_mean1(k,l)
487	dtempdztmplin = reform(reform(dtempdztmp(,,k)),[nx*ny,1])
488	alpha1out(k,l)=alpha1(k)
489	tfull=reform(tprime(,,k)+t(k,l),[nx*ny,1])
490	if (downdraft_index1(0) ne -1) then tdown1moy(k,l)=mean(tfull(downdraft_index1))
491	tplume1moy(k,l)=mean(tfull(plumeIndex1))
492	tenv1moy(k,l)=mean(tfull(envIndex1))
493	if (envIndex1_ude(0) ne -1) then tenv1moy_ude(k,l) = mean(tfull(envIndex1_ude)) else tenv1moy_ude(k,l)=0.
494	tmoy_full(k,l) = mean(tfull)
495	buoyancy1_les(k,l)=grav*(tplume1moy(k,l)/tenv1moy(k,l)-1.)
496	e_trac1_les(k,l) = fm_trac1_les(k,l)TOTAL((1./(tenv1moy(k,l)-tplume1moy(k,l)))(dtempdztmplin(plumeIndex1)),1)/float(n_elements(plumeIndex1))
497	d1_term1(k,l) = fm_trac1_les(k,l)TOTAL((1./(tenv1moy(k,l)-tplume1moy(k,l)))(temporary(dtempdztmplin(envIndex1))),1)/float(n_elements(envIndex1))
498	if (envIndex1_ude(0) ne -1) then begin
499	e1_term1_ude(k,l) = fm_trac1_les(k,l)TOTAL((1./(tenv1moy_ude(k,l)-tplume1moy(k,l)))(dtempdztmplin(plumeIndex1)),1)/float(n_elements(plumeIndex1))
500	d1_term1_ude(k,l) = fm_trac1_les(k,l)TOTAL((1./(tenv1moy_ude(k,l)-tplume1moy(k,l)))(temporary(dtempdztmplin(envIndex1_ude))),1)/float(n_elements(envIndex1_ude))
501	endif else begin
502	e1_term1_ude(k,l) = 0.
503	d1_term1_ude(k,l) = 0.
504	endelse
505	wtmp=reform(wprime(,,k)+anomalw(k),[nx*ny,1])
506	ttmp=reform(tprime(,,k)+t(k,l),[nx*ny,1])
507	qtmp=reform(zqtrac1(,,k),[nx*ny,1])
508	pttmp=reform(ptotprime(,,k),[nx*ny,1])
509	pt_mean1=mean(pttmp(plumeIndex1))
510	pttmp=0.
511	; anomalptot(,,k) = ptotprime(,,k) - temporary(pt_mean1) ;prime par rapport a la moyenne dans la plume (oscille fortement : probleme de bonne definition verticale de la plume)
512	Gamma_1_tmp(k,l) = alpha1out(k,l)rhomoy1(k,l)TOTAL((wtmp(plumeIndex1)-w_mean1(k,l))^2,1) / float(n_elements(plumeIndex1))
513	anomalptotlin = reform(anomalptot(,,k),[nx*ny,1])
514	dptotprimedztmplin = reform(dptotprimedztmp(,,k),[nx*ny,1])
515	Gamma_2(k,l) = -(TOTAL(dptotprimedztmplin(plumeIndex1),1)/float(n_elements(plumeIndex1)))/rhomoy1(k,l)
516	Gamma_3(k,l) = -grav*(TOTAL(anomalptotlin(plumeIndex1),1)/float(n_elements(plumeIndex1)))/pt(k,l)
517	hf1tmp(k,l) = TOTAL((wtmp(plumeIndex1)-w_mean1(k,l))*(ttmp(plumeIndex1)-tplume1moy(k,l)),1) / float(n_elements(plumeIndex1))
518	hf1tmpenv(k,l) = TOTAL((wtmp(envIndex1)-w_mean1_env(k,l))*(ttmp(envIndex1)-tenv1moy(k,l)),1) / float(n_elements(envIndex1))
519	q_mean_up(k,l)=mean(qtmp(plumeIndex1))
520	wq_updraft(k,l) = TOTAL((wtmp(plumeIndex1)-w_mean1(k,l))*(qtmp(plumeIndex1)-mean(qtmp(plumeIndex1))),1) / float(n_elements(plumeIndex1))
521
522	if (envIndex1_ude(0) ne -1) then begin
523	hf1tmpenv_ude(k,l) = TOTAL((wtmp(envIndex1_ude)-w_mean1_env_ude(k,l))*(ttmp(envIndex1_ude)-tenv1moy_ude(k,l)),1) / float(n_elements(envIndex1_ude))
524	q_mean_env_ude(k,l)=mean(qtmp(envIndex1_ude))
525	wq_env_ude(k,l)= TOTAL((wtmp(envIndex1_ude)-mean(wtmp(envIndex1_ude)))*(qtmp(envIndex1_ude)-mean(qtmp(envIndex1_ude))),1) / float(n_elements(envIndex1_ude))
526	endif else begin
527	hf1tmpenv_ude(k,l) =0.
528	wq_env_ude(k,l)=0.
529	endelse
530	if (downdraft_index1(0) ne -1) then begin
531	hf1tmp_down(k,l) = TOTAL((wtmp(downdraft_index1)-w_mean1_down(k,l))*(ttmp(downdraft_index1)-tdown1moy(k,l)),1) / float(n_elements(downdraft_index1))
532	q_mean_down(k,l)=mean(qtmp(downdraft_index1))
533	wq_downdraft(k,l) = TOTAL((wtmp(downdraft_index1)-mean(wtmp(downdraft_index1)))*(qtmp(downdraft_index1)-mean(qtmp(downdraft_index1))),1) / float(n_elements(downdraft_index1))
534	endif else begin
535	hf1tmp_down(k,l)=0.
536	wq_downdraft(k,l)=0.
537	endelse
538	q_mean(k,l)=mean(qtmp)
539	IF((n_elements(plumeIndex1) + n_elements(envIndex1)) ne float(nx*ny)) then print, 'WARNING : INDEX PROBLEM : plume / env : ', n_elements(plumeIndex1), n_elements(envIndex1)
540	; IF((n_elements(plumeIndex1) + n_elements(envIndex1_ude) + n_elements(downdraft_index1)) ne float(nx*ny)) then print, 'WARNING : INDEX PROBLEM : plume / env / downdraft : ', n_elements(plumeIndex1), n_elements(envIndex1_ude), n_elements(downdraft_index1)
541	ENDELSE
542	ENDFOR
543
544	; FOR i=0, nx-1 DO BEGIN
545	; FOR j=0, ny-1 DO BEGIN
546	; dptotprimedztmp(i,j,) = deriv(altitudes_LES, anomalptot(i,j,))
547	; ENDFOR
548	; ENDFOR
549	;
550	; FOR k=0, nz-1 DO BEGIN
551	; IF(plumeIndex1out(0,k) eq -1) THEN BEGIN
552	; Gamma_2(k,l)=0.
553	; ENDIF ELSE BEGIN
554	; wpi=where(plumeIndex1out ne -1)
555	; plumeIndex1=reform(plumeIndex1out(temporary(wpi),k))
556	; dptotprimedztmplin = reform(dptotprimedztmp(,,k),[nx*ny,1])
557	; Gamma_2(k,l) = -(TOTAL(dptotprimedztmplin(plumeIndex1),1)/float(n_elements(plumeIndex1)))/rhomoy1(k,l)
558	; ENDELSE
559	; ENDFOR
560
561	dgamma1tmp(,l) = deriv(altitudes_LES,Gamma_1_tmp(,l))
562	FOR k=0, nz-1 DO BEGIN
563	IF (alpha1out(k,l) ne 0.) THEN BEGIN
564	Gamma_1(k,l) = -(1./(alpha1out(k,l)rhomoy1(k,l)))dgamma1tmp(k,l)
565	ENDIF ELSE BEGIN
566	Gamma_1(k,l)=0.
567	ENDELSE
568	ENDFOR
569	drhoahfdztmp = deriv(altitudes_LES,rhomoy1(,l)alpha1out(,l)hf1tmp(*,l))
570	drhoahfdztmpDetr = deriv(altitudes_LES,rhomoy1(,l)(1.-alpha1out(,l))hf1tmpenv(*,l))
571	drhoahfdztmpDetr_ude = deriv(altitudes_LES,rhomoy1(,l)(1.-alpha1out(,l)-beta1out(,l))hf1tmpenv_ude(,l))
572
573	; wBin(*,l)=HISTOGRAM(wtmp,nbins=NBINS)
574	; wBinEnv_ude(*,l)=HISTOGRAM(wtmp(envIndex1_ude),nbins=NBINS)
575	; wBinUp(*,l)=HISTOGRAM(wtmp(plumeIndex1),nbins=NBINS)
576	; wBinDown(*,l)=HISTOGRAM(wtmp(downdraft_index1),nbins=NBINS)
577	wtmp=0.
578	ttmp=0.
579
580	FOR k=0,nz-1 DO BEGIN
581	IF(plumeIndex1out(0,k) eq -1) THEN BEGIN
582	e1_term2(k,l)=0.
583	e1_term2_ude(k,l)=0.
584	ENDIF ELSE BEGIN
585	e1_term2(k,l)=drhoahfdztmp(k)/(tenv1moy(k,l)-tplume1moy(k,l))
586	e1_term2_ude(k,l)=drhoahfdztmp(k)/(tenv1moy_ude(k,l)-tplume1moy(k,l))
587	ENDELSE
588
589	IF(envIndex1out(0,k) eq -1) THEN BEGIN
590	d1_term2(k,l)=0.
591	d1_term2_ude(k,l)=0.
592	ENDIF ELSE BEGIN
593	d1_term2(k,l)=-drhoahfdztmpDetr(k)/(tenv1moy(k,l)-tplume1moy(k,l))
594	d1_term2_ude(k,l)=-drhoahfdztmpDetr_ude(k)/(tenv1moy_ude(k,l)-tplume1moy(k,l))
595	ENDELSE
596	ENDFOR
597
598
599	tfull1=0.
600
601	zqtrac1=0.
602	; zqtrac2 = getget(filesWRF(loop), s_trac2, count=[0,0,0,1], offset=[0,0,0,loop2])
603	; FOR k=0,nz-1 DO BEGIN
604	; anomalqtrac2(,,k) = zqtrac2(,,k) - TOTAL(TOTAL(REFORM(zqtrac2(,,k)),1),1)/ float(nx) / float(ny)
605	; sigmazqtrac2(k) = STDDEV(zqtrac2(,,k))
606	; IF (k ne 0) THEN BEGIN
607	; subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
608	; sigmazminqtrac2(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac2(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
609	; ENDIF ELSE BEGIN
610	; sigmazminqtrac2(k) = sigmazqtrac2(k)
611	; ENDELSE
612	; plumeIndex2 = WHERE((anomalqtrac2(,,k) GT scaleMAX([sigmazqtrac2(k),sigmazminqtrac2(k)])) AND ((wprime(,*,k)+anomalw(k)) GT 0.))
613	; envIndex2 = WHERE((anomalqtrac2(,,k) LE scaleMAX([sigmazqtrac2(k),sigmazminqtrac2(k)])) OR ((wprime(,*,k)+anomalw(k)) LE 0.))
614	; IF(plumeIndex2(0) EQ -1) THEN BEGIN
615	; fm_trac2_les(k,l)=0.
616	; e_trac2_les(k,l)=0.
617	; alpha2out(k,l)=0.
618	; buoyancy2_les(k,l)=0.
619	; w_mean2(k,l)=0.
620	; ENDIF ELSE BEGIN
621	; alpha2(k) = n_elements(plumeIndex2) / float(nx) / float(ny)
622	; wprimetmp = reform(reform((anomalw(k)+wprime(,,k))),[nx*ny,1])
623	; w_mean2(k,l) = mean(wprimetmp(plumeIndex2))
624	; wprimetmp=0.
625	; fm_trac2_les(k,l) = alpha2(k)rhomoy1(k,l)w_mean2(k,l)
626	; tprimetmp = reform(reform(-tprime(,,k)),[nx*ny,1])
627	; dtempdztmplin = reform(reform(dtempdztmp(,,k)),[nx*ny,1])
628	; e_trac2_les(k,l) = TOTAL((1./(temporary(tprimetmp(plumeIndex2))))*(temporary(dtempdztmplin(plumeIndex2))),1)/n_elements(plumeIndex2)
629	; alpha2out(k,l)=alpha2(k)
630	; tfull=reform(tprime(,,k)+t(k,l),[nx*ny,1])
631	; tplume2moy=mean(tfull(plumeIndex2))
632	; tenv2moy=mean(tfull(envIndex2))
633	; buoyancy2_les(k,l)=grav*(tplume2moy/tenv2moy-1.)
634	; ENDELSE
635	; ENDFOR
636	; zqtrac2=0.
637
638	wt(,l) = TOTAL(TOTAL(TEMPORARY(tprime)wprime,1),1) / float(nx) / float(ny)
639	wq(,l) = TOTAL(TOTAL(TEMPORARY(wprime)anomalqtrac1,1),1) / float(nx) / float(ny)
640
641	wmax(l) = max(w_mean1(*,l))
642	l=l+1
643	ENDFOR
644	print, 'file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s'
645
646	ENDFOR
647
648	hf1_term1 = hf1tmp*alpha1out
649	hf1_term2 = temporary(hf1tmpenv)*(1.-alpha1out)
650	hf1_term3 = alpha1out(1.-alpha1out)(w_mean1 - w_mean1_env)*(tplume1moy - tenv1moy)
651
652	hf1_ude_term1 = temporary(hf1tmp)*alpha1out
653	hf1_ude_term2 = temporary(hf1tmp_down)*beta1out
654	hf1_ude_term3 = temporary(hf1tmpenv_ude)*(1.-(alpha1out+beta1out))
655	hf1_ude_term4 = alpha1out(w_mean1 - w_mean1_full)(tplume1moy - tmoy_full) + beta1out(w_mean1_down - w_mean1_full)(tdown1moy - tmoy_full) + (1.- (alpha1out+beta1out))(w_mean1_env_ude - w_mean1_full)(tenv1moy_ude - tmoy_full)
656
657	FOR k=0, nz-1 DO BEGIN
658	; dteta1moydt_entr(k,) = deriv(localtime,tplume1moy(k,))/3700. - dTeta_phys(k,*)
659	; dteta1moydt_detr(k,) = deriv(localtime,tplume1moy(k,))/3700. + dTeta_phys(k,*)
660	dteta1moydt_entr(k,) = deriv(localtime,tplume1moy(k,))/3700. - dTeta_phys(k,*)
661	dteta1moydt_detr(k,) = dTeta_phys(k,) - deriv(localtime,tplume1moy(k,*))/3700.
662	ENDFOR
663
664	FOR k=0, nz-1 DO BEGIN
665	FOR l=0, nttot-1 DO BEGIN
666	IF (tenv1moy(k,l) ne tplume1moy(k,l)) THEN e1_term3(k,l) = rhomoy1(k,l)alpha1out(k,l)dteta1moydt_entr(k,l)/(tenv1moy(k,l)-tplume1moy(k,l)) ELSE e1_term3(k,l)=0.
667	IF (tenv1moy_ude(k,l) ne tplume1moy(k,l)) THEN e1_term3_ude(k,l) = rhomoy1(k,l)alpha1out(k,l)dteta1moydt_entr(k,l)/(tenv1moy_ude(k,l)-tplume1moy(k,l)) ELSE e1_term3_ude(k,l)=0
668	; IF (tenv1moy(k,l) ne tplume1moy(k,l)) THEN d1_term3(k,l) = rhomoy1(k,l)(1.-alpha1out(k,l))dteta1moydt_detr(k,l)/(tenv1moy(k,l)-tplume1moy(k,l)) ELSE d1_term3(k,l)=0.
669	IF (tenv1moy(k,l) ne tplume1moy(k,l)) THEN d1_term3(k,l) = rhomoy1(k,l)(1.-alpha1out(k,l))dteta1moydt_detr(k,l)/(tenv1moy(k,l)-tplume1moy(k,l)) ELSE d1_term3(k,l)=0.
670	IF (tenv1moy_ude(k,l) ne tplume1moy(k,l)) THEN d1_term3_ude(k,l) = rhomoy1(k,l)(1.-alpha1out(k,l)-beta1out(k,l))dteta1moydt_detr(k,l)/(tenv1moy_ude(k,l)-tplume1moy(k,l)) ELSE d1_term3_ude(k,l)=0.
671	ENDFOR
672	ENDFOR
673
674
675	ht = TEMPORARY(pht) - hgtu/1000.
676	save, hfx, flxrad, flxgrd, lwdownz, swdownz, q_mean_up, q_mean_down, q_mean_env_ude, q_mean, Gamma_1, Gamma_2, Gamma_3, w_mean1_env, tsurf, wq, wq_updraft, wq_downdraft, wq_env_ude, d1_term1_ude, d1_term2_ude, d1_term3_ude, e1_term1_ude, e1_term2_ude, e1_term3_ude, tplume1moy, tdown1moy, w_mean1_full, tmoy_full, tenv1moy_ude, w_mean1_env_ude, uv_moy, hf1_ude_term1, hf1_ude_term2, hf1_ude_term3, hf1_ude_term4, w_mean1_down, downward_flux1, beta1out, hf1_term1, hf1_term2, hf1_term3, d1_term1, d1_term2, d1_term3, e1_term2, e1_term3, buoyancy1_les, buoyancy2_les, w_mean1, w_mean2, nx, ny, alpha1out, alpha2out, e_trac1_les, e_trac2_les, tke_les, ztke, altitudes_LES, ht, t, p, pt, localtime, xtke, ytke, wt, temp_les, wmax, fm_trac1_les, fm_trac2_les,filename=les_path+'/'+datname
677
678	nz = n_elements(altitudes_LES)
679
680	ENDIF ELSE BEGIN
681
682	print, 'OK, file is here'
683	restore, filename=les_path+'/'+datname
684	nz = n_elements(altitudes_LES)
685	nttot = n_elements(tmoy_full(0,*))
686
687	OPENR, 23, les_path+'/input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
688
689	ENDELSE
690
691	tenv1moy = tplume1moy/((buoyancy1_les/grav)+1.)
692
693	taverage = string((localtime(nstot)-localtime(1))*3700./60.)
694	print, ''
695	print, ' -- Loading testphys1d data -- '
696
697	file1=gcm_path+'/diagfi.nc'
698	file2=gcm_convadj_path+'/diagfi.nc'
699	file3='/san0/acolmd/SIMUS/GCM3D_TestBed/diagfi.nc'
700
701	getcdf, file=file1, charvar='q2', invar=tke_gcm
702	getcdf, file=file1, charvar='aps', invar=aps
703	getcdf, file=file1, charvar='bps', invar=bps
704	getcdf, file=file1, charvar='co2col', invar=co2_col
705	;getcdf, file=file1, charvar='arcol', invar=ar_col
706	getcdf, file=file1, charvar='tkecol', invar=tke_col
707	;getcdf, file=file1, charvar='ar', invar=ar
708	getcdf, file=file1, charvar='heatFlux_up', invar=heatFlux_up
709	getcdf, file=file1, charvar='heatFlux_down', invar=heatFlux_down
710	getcdf, file=file1, charvar='pplay', invar=pplay
711	getcdf, file=file1, charvar='pplev', invar=pplev
712	getcdf, file=file1, charvar='temp', invar=temp_gcm
713	getcdf, file=file1, charvar='zw2', invar=zw2_lev
714	getcdf, file=file1, charvar='fm_therm', invar=fm_therm_gcm_lev
715	getcdf, file=file1, charvar='entr_therm', invar=zdz_entr_therm_gcm
716	getcdf, file=file1, charvar='detr_therm', invar=zdz_detr_therm_gcm
717	getcdf, file=file1, charvar='fraca', invar=alpha_gcm_lev
718	getcdf, file=file1, charvar='buoyancyOut', invar=buoyancy_gcm
719	getcdf, file=file1, charvar='buoyancyEst', invar=buoyancy_est_gcm
720	getcdf, file=file1, charvar='zkh', invar=zkh
721	getcdf, file=file1, charvar='zh', invar=zh
722	getcdf, file=file1, charvar='tsurf', invar=tsurf_gcm
723
724	if (overplot_convadj eq 'true') then begin
725	getcdf, file=file2, charvar='temp', invar=temp_gcm_convadj
726	getcdf, file=file2, charvar='pplay', invar=pplay_convadj
727	endif
728	if (plot_3d eq 'true') then begin
729	getcdf, file=file3, charvar='temp', invar=temp_gcm_3d
730	getcdf, file=file3, charvar='pplay', invar=pplay_3d
731	getcdf, file=file3, charvar='latitude', invar=latitude_3d
732	getcdf, file=file3, charvar='longitude', invar=longitude_3d
733	nWEmx_3d = n_elements(reform(temp_gcm_3d(*,0,0,0)))
734	nNSmx_3d = n_elements(reform(temp_gcm_3d(0,*,0,0)))
735	nZmx_3d = n_elements(reform(temp_gcm_3d(0,0,*,0)))
736	nTmx_3d = n_elements(reform(temp_gcm_3d(0,0,0,*)))
737	ndays_3d = 1.
738	lctu_gcm_3d = 0.
739	history_interval_s_gcm_3d = ndays_3d*88800./float(nTmx_3d) ; Timestep interval of gcm 1d simu in sec
740	localtime_lon0 = lctu_gcm_3d + history_interval_s_gcm_3d*findgen(nTmx_3d)/3700.
741	Xc = 205.
742	Yc = 21.8
743	plot_index_x = (Xc-longitude_3d(0))/(longitude_3d(1)-longitude_3d(0))
744	plot_index_y = (Yc-latitude_3d(0))/(latitude_3d(1)-latitude_3d(0))
745	localtime_true = localtime_lon0 -(12./180.)*Xc
746	endif
747
748
749	nTmx = n_elements(reform(temp_gcm(0,*)))
750	if (overplot_convadj eq 'true') then begin
751	nTmx_convadj = n_elements(reform(temp_gcm_convadj(0,*)))
752	endif else begin
753	nTmx_convadj = 10000.
754	endelse
755	ndays = 1.
756	print, ''
757	print, 'WARNING ----------------------- '
758	print, 'CONFIGURATION : '+string(ndays,format='(I0)')+' days simulation'
759	print, ''
760	history_interval_s_gcm = ndays*88800./float(nTmx) ; Timestep interval of gcm 1d simu in sec
761	history_interval_s_gcm_convadj = ndays*88800./float(nTmx_convadj)
762	localtime_gcm = lctu_gcm + history_interval_s_gcm*findgen(nTmx)/3700.
763	localtime_gcm_convadj = lctu_gcm + history_interval_s_gcm_convadj*findgen(nTmx_convadj)/3700.
764	; **********************************
765	; ****** PLOTS ****************
766	if (f_offset eq 'true') then begin
767	offset_localtime = 3.108100
768	endif else begin
769	offset_localtime = 0.
770	endelse
771	localtime=localtime+offset_localtime
772	print, '****************************************************************************************************'
773	print, 'local time LES'
774	print, localtime
775	print, '****************************************************************************************************'
776	print, 'local time GCM'
777	print, localtime_gcm
778	print, '****************************************************************************************************'
779
780	time_offset = (ndays-1.)*24.
781
782	lt_plot_ini = 8.
783	lt_plotindex_les_ini = where(localtime eq lt_plot_ini)
784	lt_plotindex_gcm_ini = where(localtime_gcm eq (lt_plot_ini+time_offset))
785
786	lt_plot0 = 10.
787	lt_plotindex_les0 = where(localtime eq lt_plot0)
788	lt_plotindex_gcm0 = where(localtime_gcm eq (lt_plot0+time_offset))
789	lt_plotindex_gcm_convadj0 = where(localtime_gcm_convadj eq (lt_plot0+time_offset))
790
791	;lt_plot = 8.25
792	lt_plot = 12.
793	lt_plotindex_les = where((localtime lt lt_plot+0.01) and (localtime gt lt_plot-0.01))
794	lt_plotindex_gcm = where(localtime_gcm eq (lt_plot+time_offset))
795	lt_plotindex_gcm_convadj = where(localtime_gcm_convadj eq (lt_plot+time_offset))
796	print, 'lt plotindex les 12h'
797	print, lt_plotindex_les
798
799	lt_plot2 = 14.
800	lt_plotindex_les2 = where(localtime eq lt_plot2)
801	lt_plotindex_gcm2 = where(localtime_gcm eq (lt_plot2+time_offset))
802	lt_plotindex_gcm_convadj2 = where(localtime_gcm_convadj eq (lt_plot2+time_offset))
803
804	lt_plot3 = 16.
805	lt_plotindex_les3 = where(localtime eq lt_plot3)
806	lt_plotindex_gcm3 = where(localtime_gcm eq (lt_plot3+time_offset))
807	lt_plotindex_gcm_convadj3 = where(localtime_gcm_convadj eq (lt_plot3+time_offset))
808
809	lt_plot4 = 18.
810	lt_plotindex_les4 = where(localtime eq lt_plot4)
811	lt_plotindex_gcm4 = where(localtime_gcm eq (lt_plot4+time_offset))
812	lt_plotindex_gcm_convadj4 = where(localtime_gcm_convadj eq (lt_plot4+time_offset))
813	;--------------------------------------------------------------------------------
814	;---------------------------------------------------------------------------------
815
816	nTmx_les=n_elements(reform(wt(0,*)))
817	nZmx=n_elements(aps) ; number of vertical levels
818	H_low=9650. ; scale height at low altitudes
819	H_high=15000. ; scale height at high altitudes
820	trans_window=10. ; # of levels over which H(:) goes from H_low to H_high
821	lev_trans=32.+trans_window/2. ; level at which H(lev_trans)=(H_low+H_high)/2
822	P_ref=p0 ; reference surface pressure used to build zsurface -610 Pa-
823	Hgcm = make_array(nZmx)
824	altitudes_GCM = make_array(nZmx)
825	; Build scale heights
826	;FOR k=0,nZmx-1 DO BEGIN
827	; Hgcm(k)=H_low+(H_high-H_low)0.5(1.0+tanh(6.*(k-lev_trans)/trans_window))
828	;ENDFOR
829
830	FOR k=0,nZmx-1 DO BEGIN
831	Hgcm(k)=R*temp_gcm(k,lt_plotindex_gcm)/grav
832	ENDFOR
833	print, 'Hgcm'
834	print, Hgcm
835	; Compute altitudes_GCM
836	FOR k=0,nZmx-1 DO BEGIN
837	altitudes_GCM(k)=-Hgcm(k)*alog(pplay(k,lt_plotindex_gcm)/pGround)
838	ENDFOR
839	Hgcm=0.
840
841	teta_gcm = temp_gcm * (p0/pplay)^r_cp
842	if (overplot_convadj eq 'true') then begin
843	teta_gcm_convadj = temp_gcm_convadj * (p0/pplay_convadj)^r_cp
844	endif
845
846	OPENR, 1, gcm_path+'/profile'
847	data=FLTARR(nZmx+1)
848	READF, 1, data
849	temp_gcm_0_ground = data(0)
850	temp_gcm_0 = data(1:nZmx-1)
851	data = 0.
852	CLOSE, 1
853
854	teta_gcm_0 = temp_gcm_0 * (p0/pplay)^r_cp
855	approx_zdz_gcm = make_array(nZmx)
856	approx_zdz_gcm(0)=altitudes_GCM(1)
857	FOR k=1, nZmx-2 DO BEGIN
858	approx_zdz_gcm(k) = altitudes_GCM(k+1) - altitudes_GCM(k)
859	ENDFOR
860	approx_zdz_gcm(nZmx-1)=approx_zdz_gcm(nZmx-2)
861
862	print, 'approx zdz gcm'
863	print, approx_zdz_gcm
864
865	print, '****************************************************************************************************'
866	print, 'altitudes LES based on phtot : inter-levels'
867	print, altitudes_LES
868	print, '****************************************************************************************************'
869	print, 'altitudes GCM based on pplay : inter-levels'
870	print, altitudes_GCM
871	print, '****************************************************************************************************'
872
873	; Compute tracer deviation :
874
875	co2_col = co2_col/co2_col(0)
876	;ar_col = ar_col/ar_col(0)
877	tke_col = tke_col+1.
878
879	; Compute <teta> les
880
881	teta_les = temporary(t)
882
883
884	; ========================================================================
885	; ========================================================================
886
887	IF (visualization_mode eq 'true') THEN BEGIN
888
889	print,' *****************************************-----------------------------------'
890	print,' ********** PLUME ********************-----------------------------------'
891	print,' *****************************************-----------------------------------'
892
893	; We are evaluating the first time-step element of the file number 'loop-1' :
894	; file 1 starts at 8h (loop =0, loop2 =0)
895	; file 6 starts at 13h (roughly) (loop =5, loops2=0)
896	; file 12 starts at 18h (roughly) (loop = 11,loop2 = 0)
897
898	loop=uint(loop_special)-1
899	;loop2=34
900	loop2=10
901	domain='d01'
902	filesWRF = FindFile(les_path+'/wrfout_'+domain+'_????-??-??_??:??:??')
903	anomalw=1.
904
905	zqtrac1 = dblarr(nx,ny,nz) & zqtrac2 = dblarr(nx,ny,nz)
906	sigmazqtrac1 = fltarr(nz) & sigmazqtrac2 = fltarr(nz)
907	sigmazminqtrac1 = fltarr(nz) & sigmazminqtrac2 = fltarr(nz)
908	anomalqtrac1 = fltarr(nx,ny,nz) & anomalqtrac2 = fltarr(nx,ny,nz)
909	zqtrac1 = getget(filesWRF(loop), 'qtrac1', count=[0,0,0,1], offset=[0,0,0,loop2])
910	zqtrac2 = getget(filesWRF(loop), 'qtrac2', count=[0,0,0,1], offset=[0,0,0,loop2])
911	wprime = getget(filesWRF(loop), 'W', anomaly = anomalw, count=[0,0,0,1], offset=[0,0,0,loop2])
912	supermask1 = fltarr(nx,ny,nz)
913	supermask2 = fltarr(nx,ny,nz)
914	k_out_histo = 8
915	k_out_hist = [1,10,25,50,70,85]
916	nbtest=n_elements(k_out_hist)
917	b=0
918	FOR k=0,nz-1 DO BEGIN
919	mask1=fltarr(nx*ny)
920	mask2=fltarr(nx*ny)
921	anomalqtrac1(,,k) = zqtrac1(,,k) - TOTAL(TOTAL(REFORM(zqtrac1(,,k)),1),1)/ float(nx) / float(ny)
922	sigmazqtrac1(k) = STDDEV(REFORM(zqtrac1(,,k)))
923	IF (k ne 0) THEN BEGIN
924	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
925	sigmazminqtrac1(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac1(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
926	ENDIF ELSE BEGIN
927	sigmazminqtrac1(k) = sigmazqtrac1(k)
928	ENDELSE
929	print, sigmazqtrac1(k),sigmazminqtrac1(k)
930	plumeIndex1 = WHERE((anomalqtrac1(,,k) GT scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) AND ((anomalw(k)+wprime(,*,k)) GT 0.))
931	; anomalqtrac2(,,k) = zqtrac2(,,k) - TOTAL(TOTAL(REFORM(zqtrac2(,,k)),1),1)/ float(nx) / float(ny)
932	; sigmazqtrac2(k) = STDDEV(REFORM(zqtrac2(,,k)))
933	IF (k ne 0) THEN BEGIN
934	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
935	; sigmazminqtrac2(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac2(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
936	ENDIF ELSE BEGIN
937	; sigmazminqtrac2(k) = sigmazqtrac2(k)
938	ENDELSE
939	; plumeIndex2 = WHERE((anomalqtrac2(,,k) GT scaleMAX([sigmazqtrac2(k),sigmazminqtrac2(k)])) AND ((anomalw(k)+wprime(,*,k)) GT 0.))
940	IF(plumeIndex1(0) NE -1 ) THEN BEGIN
941	mask1(plumeIndex1) = 1.
942	ENDIF ELSE BEGIN
943	mask1(*)=0.
944	ENDELSE
945	IF(plumeIndex2(0) NE -1 ) THEN BEGIN
946	mask2(plumeIndex2) = 1.
947	ENDIF ELSE BEGIN
948	mask2(*)=0.
949	ENDELSE
950	mask1 = reform(mask1,[nx,ny])
951	supermask1(,,k) = mask1(,)
952	mask2 = reform(mask2,[nx,ny])
953	supermask2(,,k) = mask2(,)
954	; IF (k eq k_out_histo) THEN BEGIN
955	; plume1_out = plumeIndex1
956	; plume2_out = plumeIndex2
957	; ENDIF
958
959	IF (k eq k_out_hist(0)) THEN BEGIN
960	c1=plumeIndex1
961	c2=plumeIndex2
962	ENDIF
963	IF (k eq k_out_hist(1)) THEN BEGIN
964	d1=plumeIndex1
965	d2=plumeIndex2
966	ENDIF
967	IF (k eq k_out_hist(2)) THEN BEGIN
968	e1=plumeIndex1
969	e2=plumeIndex2
970	ENDIF
971	IF (k eq k_out_hist(3)) THEN BEGIN
972	f1=plumeIndex1
973	f2=plumeIndex2
974	ENDIF
975	IF (k eq k_out_hist(4)) THEN BEGIN
976	g1=plumeIndex1
977	g2=plumeIndex2
978	ENDIF
979	IF (k eq k_out_hist(5)) THEN BEGIN
980	h1=plumeIndex1
981	h2=plumeIndex2
982	ENDIF
983	ENDFOR
984
985	IF (Histo eq 'false') THEN BEGIN
986	IVOLUME, supermask1
987	IVOLUME, supermask2
988	ENDIF ELSE BEGIN
989
990	; -------------------------------------------------------------------------------
991	; THIS IS THE ULTRA-GORE UN-ESTHETIC UGLY-AS-HELL LOOP FOR CONCENTRATION PLOTTING
992	; but well, this is just because idl cannot handle arrays as well as I would like
993	; -------------------------------------------------------------------------------
994	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Distrib.ps'
995	PSOPEN, THICK=100, CHARSIZE=60, FILE = filename, FONT = 5, TFONT = 5,XPLOTS=2,YPLOTS=3,MARGIN=2000
996
997	FOR n=0, nbtest-1 DO BEGIN
998
999	CASE n OF
1000	0:BEGIN
1001	plume1_out=c1 & plume2_out=c2
1002	END
1003	1:BEGIN
1004	plume1_out=d1 & plume2_out=d2
1005	END
1006	2:BEGIN
1007	plume1_out=e1 & plume2_out=e2
1008	END
1009	3:BEGIN
1010	plume1_out=f1 & plume2_out=f2
1011	END
1012	4:BEGIN
1013	plume1_out=g1 & plume2_out=g2
1014	END
1015	5:BEGIN
1016	plume1_out=h1 & plume2_out=h2
1017	END
1018	ENDCASE
1019
1020	ToBin1 = reform(zqtrac1(,,k_out_hist(n)),[nx*ny,1])
1021	;ToBin2 = reform(zqtrac2(,,k_out_hist(n)),[nx*ny,1])
1022	svmin1=min(ToBin1) & svmin2=min(ToBin2)
1023	svmax1=max(ToBin1) & svmax2=max(ToBin2)
1024	NBINS=100
1025	ds1=(svmax1-svmin1+1.)/(NBINS-1) & ds2=(svmax2-svmin2+1.)/(NBINS-1)
1026	Xaxis1 = svmin1+((svmax1-svmin1)/(NBINS-1))*indgen(NBINS)
1027	Xaxis2 = svmin2+((svmax2-svmin2)/(NBINS-1))*indgen(NBINS)
1028	Bin1=HISTOGRAM(ToBin1,nbins=NBINS)
1029	Bin2=INTERPOL(HISTOGRAM(ToBin2,nbins=NBINS),Xaxis2,Xaxis1,/SPLINE)
1030
1031	what_I_plot = [[Bin1],[Bin2]]
1032	labels=['LES tracer 1 conc. distrib.','LES tracer 2 conc. distrib.']
1033
1034	title_user = TestCase+SubCase+' LES tracer 1&2 concentration distribution at '+string(altitudes_LES(k_out_hist(n)))+'m AGL'
1035	IF (n lt 3) THEN BEGIN
1036	POS, XPOS=1, YPOS=uint(n+1)
1037	ENDIF ELSE BEGIN
1038	POS, XPOS=2, YPOS=uint(n+1)-3
1039	ENDELSE
1040	MAP
1041	CS, SCALE=28
1042	GSET, XMIN=0, XMAX=20, YMIN=0, YMAX=300, TITLE=title_user
1043	cols=INDGEN(2)+2
1044	GPLOT, X=Xaxis1, Y=what_I_plot, /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
1045	AXES, XSTEP = 4, XTITLE='Trac concentration (kg/kg)', YSTEP=50, YTITLE='Counts',NDECS=1
1046
1047	ToBin1 = reform(zqtrac1(,,k_out_hist(n)),[nx*ny,1])
1048	ToBin2 = reform(zqtrac2(,,k_out_hist(n)),[nx*ny,1])
1049	ToBin1 = ToBin1(plume1_out)
1050	ToBin2 = ToBin2(plume2_out)
1051	svmin1=min(ToBin1) & svmin2=min(ToBin2)
1052	svmax1=max(ToBin1) & svmax2=max(ToBin2)
1053	NBINS=50
1054	ds1=(svmax1-svmin1+1.)/(NBINS-1) & ds2=(svmax2-svmin2+1.)/(NBINS-1)
1055	Xaxis1 = svmin1+((svmax1-svmin1)/(NBINS-1))*indgen(NBINS)
1056	Xaxis2 = svmin2+((svmax2-svmin2)/(NBINS-1))*indgen(NBINS)
1057	Bin1=HISTOGRAM(ToBin1,nbins=NBINS)
1058	Bin2=HISTOGRAM(ToBin2,nbins=NBINS)
1059	oplot, Xaxis1, Bin1, psym=4
1060	oplot, Xaxis2, Bin2, psym=5
1061
1062	ENDFOR
1063
1064	PSCLOSE, /NOVIEW
1065	spawn, 'ps2png '+filename
1066
1067	ENDELSE
1068
1069	ENDIF ELSE BEGIN
1070
1071	; =========================================================================
1072	; =========================================================================
1073
1074	print, '........ ALPHA'
1075
1076	alpha_interlay_gcm = make_array(nZmx)
1077	FOR k=0, nZmx-2 DO BEGIN
1078	alpha_interlay_gcm(k) = (alpha_gcm_lev(k,lt_plotindex_gcm)+alpha_gcm_lev(k+1,lt_plotindex_gcm))/2.
1079	ENDFOR
1080	alpha_interlay_gcm(nZmx-1)=0.
1081
1082	smoothed_alpha1_les = make_array(nz)
1083	smoothed_alpha2_les = make_array(nz)
1084	smoothed_beta1_les = make_array(nz)
1085	FOR t=-ns,ns DO BEGIN
1086	smoothed_alpha1_les = smoothed_alpha1_les + REFORM(alpha1out(*,lt_plotindex_les+t))
1087	smoothed_alpha2_les = smoothed_alpha2_les + REFORM(alpha2out(*,lt_plotindex_les+t))
1088	smoothed_beta1_les = smoothed_beta1_les + REFORM(beta1out(*,lt_plotindex_les+t))
1089	ENDFOR
1090
1091	smoothed_alpha1_les = smoothed_alpha1_les/nstot
1092	smoothed_alpha2_les = smoothed_alpha2_les/nstot
1093	smoothed_beta1_les = smoothed_beta1_les/nstot
1094
1095	; =========================================================================
1096
1097	; * Temperatures *
1098
1099	print, '........ TEMPERATURES'
1100
1101	xmin = 190
1102	xmax = 250
1103	if (TestCase eq 'Case_Z') then begin
1104	xmin = 170
1105	xmax = 250
1106	endif
1107
1108
1109	what_I_plot = [[reform(temp_gcm(,lt_plotindex_gcm_ini))],[reform(temp_gcm(,lt_plotindex_gcm0))],[reform(temp_gcm(,lt_plotindex_gcm))],[reform(temp_gcm(,lt_plotindex_gcm2))],[reform(temp_gcm(,lt_plotindex_gcm3))],[reform(temp_gcm(,lt_plotindex_gcm4))]]
1110	labels=['TH temp 1d, lt='+string(lt_plot_ini),'TH temp 1d, lt='+string(lt_plot0),'TH temp 1d, lt='+string(lt_plot),'TH temp 1d, lt='+string(lt_plot2),'TH temp 1d, lt='+string(lt_plot3),'TH temp 1d, lt='+string(lt_plot4)]
1111	title_user = TestCase+SubCase+LayerCase+' Temperatures Comparison'
1112	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Temperature.ps'
1113	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1114	CS, SCALE=28
1115	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=12, TITLE=title_user
1116	cols=INDGEN(6)+2
1117	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1118	AXES, XSTEP = 5, XTITLE='Temperature (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
1119
1120	oplot, temp_les(*,lt_plotindex_les_ini), altitudes_LES/1000., psym=4
1121	oplot, temp_les(*,lt_plotindex_les0), altitudes_LES/1000., psym=4
1122	oplot, temp_les(*,lt_plotindex_les), altitudes_LES/1000., psym=4
1123	oplot, temp_les(*,lt_plotindex_les2), altitudes_LES/1000., psym=4
1124	oplot, temp_les(*,lt_plotindex_les3), altitudes_LES/1000., psym=4
1125	oplot, temp_les(*,lt_plotindex_les4), altitudes_LES/1000., psym=4
1126
1127	if(overplot_convadj EQ 'true') then begin
1128	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj0), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1129	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1130	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj2), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1131	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj3), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1132	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj4), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1133	endif
1134
1135
1136	PSCLOSE, /NOVIEW
1137
1138	spawn, 'ps2png '+filename
1139
1140	; * Pressions *
1141
1142	print, '........ PRESSURES'
1143
1144	;what_I_plot = make_array(nZmx)
1145	;labels=['TH P 1d, lt='+string(lt_plot)]
1146	;title_user = TestCase+SubCase+' Pressure Comparisons'
1147	;filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Pressure.ps'
1148	;PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1149	;CS, SCALE=28
1150	;GSET, XMIN=600, XMAX=900, YMIN=0, YMAX=0.5, TITLE=title_user
1151	;cols=INDGEN(1)+2
1152	;GPLOT, X=pplay(,lt_plotindex_gcm), Y=-alog(pplay(,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1153	;AXES, XSTEP = 25, XTITLE='Log P', YSTEP=0.1, YTITLE='Altitude (km)',NDECS=1
1154	;
1155	;oplot, pt(,lt_plotindex_les),-alog(pt(,lt_plotindex_les)/pGround), psym=4
1156	;
1157	;PSCLOSE, /NOVIEW
1158	;
1159	;spawn, 'ps2png '+filename
1160
1161	what_I_plot = make_array(nZmx)
1162	labels=['TH P 1d, lt='+string(lt_plot)]
1163	title_user = TestCase+SubCase+LayerCase+' Pressure Comparisons'
1164	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Pressure.ps'
1165	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1166	CS, SCALE=28
1167	GSET, XMIN=400, XMAX=900, YMIN=0, YMAX=10, TITLE=title_user
1168	cols=INDGEN(1)+2
1169	GPLOT, X=pplay(*,lt_plotindex_gcm), Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1170	AXES, XSTEP = 25, XTITLE='Pression (Pa)', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
1171
1172	oplot, pt(*,lt_plotindex_les),altitudes_LES/1000., psym=4
1173
1174	PSCLOSE, /NOVIEW
1175
1176	spawn, 'ps2png '+filename
1177
1178
1179	; * Temperatures potentielles *
1180	if(full eq 'true') then begin
1181
1182	xmin = 185
1183	xmax = 240
1184	if (TestCase eq 'Case_C') then begin
1185	xmin = 200
1186	xmax = 255
1187	endif
1188	if (TestCase eq 'Case_I') then begin
1189	xmin = 195
1190	xmax = 250
1191	endif
1192	if (TestCase eq 'Case_Z') then begin
1193	xmin = 230
1194	xmax = 290
1195	endif
1196	if (TestCase eq 'ExtremeCase') then begin
1197	xmin = 195
1198	xmax = 260
1199	endif
1200
1201	print, '........ POTENTIAL TEMPERATURES'
1202	what_I_plot = [[reform(teta_gcm(,lt_plotindex_gcm_ini))],[reform(teta_gcm(,lt_plotindex_gcm0))],[reform(teta_gcm(,lt_plotindex_gcm))],[reform(teta_gcm(,lt_plotindex_gcm2))],[reform(teta_gcm(,lt_plotindex_gcm3))],[reform(teta_gcm(,lt_plotindex_gcm4))]]
1203	labels=['TH teta 1d, lt='+string(lt_plot_ini),'TH teta 1d, lt='+string(lt_plot0),'TH teta 1d, lt='+string(lt_plot),'TH teta 1d, lt='+string(lt_plot2),'TH teta 1d, lt='+string(lt_plot3),'TH teta 1d, lt='+string(lt_plot4)]
1204	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1205	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta.ps'
1206	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1207	CS, SCALE=28
1208	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=1.4, TITLE=title_user
1209	cols=INDGEN(6)+2
1210	GPLOT, X=what_I_plot, Y=-alog(pplay(*,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1211	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.2, YTITLE='-Log(P/P0) ',NDECS=1
1212
1213	oplot, teta_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1214	oplot, teta_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1215	oplot, teta_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1216	oplot, teta_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1217	oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1218	oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1219	if(overplot_convadj EQ 'true') then begin
1220	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj0), -alog(pplay_convadj(,lt_plotindex_gcm_convadj0)/pGround), thick=0.1,color=8,linestyle=3
1221	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj), -alog(pplay_convadj(,lt_plotindex_gcm_convadj)/pGround), thick=0.1,color=8,linestyle=3
1222	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj2), -alog(pplay_convadj(,lt_plotindex_gcm_convadj2)/pGround), thick=0.1,color=8,linestyle=3
1223	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj3), -alog(pplay_convadj(,lt_plotindex_gcm_convadj3)/pGround), thick=0.1,color=8,linestyle=3
1224	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj4), -alog(pplay_convadj(,lt_plotindex_gcm_convadj4)/pGround), thick=0.1,color=8,linestyle=3
1225	endif
1226	oplot, teta_gcm_0(), -alog(pplay(,lt_plotindex_gcm)/pGround)
1227
1228	PSCLOSE, /NOVIEW
1229
1230	spawn, 'ps2png '+filename
1231
1232	xmin = 230
1233	xmax = 245
1234	what_I_plot = [[reform(teta_gcm(,lt_plotindex_gcm)),tsurf_gcm(lt_plotindex_gcm)],[reform(teta_gcm(,lt_plotindex_gcm2)),tsurf_gcm(lt_plotindex_gcm2)]]
1235	labels=['TH teta 1d, lt='+string(lt_plot),'TH teta 1d, lt='+string(lt_plot2)]
1236	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1237	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta_zoom.ps'
1238	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1239	CS, SCALE=28
1240	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=0.05, TITLE=title_user
1241	cols=INDGEN(2)+2
1242	GPLOT, X=what_I_plot, Y=[0.,-alog(pplay(*,lt_plotindex_gcm)/pGround)], /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1243	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.005, YTITLE='-Log(P/P0) ',NDECS=1
1244
1245	;oplot, teta_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1246	;oplot, teta_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1247	oplot, teta_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1248	oplot, teta_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1249	;oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1250	;oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1251	if(overplot_convadj EQ 'true') then begin
1252	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj0), -alog(pplay_convadj(,lt_plotindex_gcm_convadj0)/pGround), thick=0.1,color=8,linestyle=3
1253	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj), -alog(pplay_convadj(,lt_plotindex_gcm_convadj)/pGround), thick=0.1,color=8,linestyle=3
1254	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj2), -alog(pplay_convadj(,lt_plotindex_gcm_convadj2)/pGround), thick=0.1,color=8,linestyle=3
1255	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj3), -alog(pplay_convadj(,lt_plotindex_gcm_convadj3)/pGround), thick=0.1,color=8,linestyle=3
1256	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj4), -alog(pplay_convadj(,lt_plotindex_gcm_convadj4)/pGround), thick=0.1,color=8,linestyle=3
1257	endif
1258	;oplot, teta_gcm_0(), -alog(pplay(,lt_plotindex_gcm)/pGround)
1259
1260	PSCLOSE, /NOVIEW
1261
1262	spawn, 'ps2png '+filename
1263
1264
1265	if (plot_3d eq 'true') then begin
1266	what_I_plot = [[reform(teta_gcm(,lt_plotindex_gcm_ini))],[reform(teta_gcm(,lt_plotindex_gcm0))],[reform(teta_gcm(,lt_plotindex_gcm))],[reform(teta_gcm(,lt_plotindex_gcm2))],[reform(teta_gcm(,lt_plotindex_gcm3))],[reform(teta_gcm(,lt_plotindex_gcm4))]]
1267	labels=['TH teta 1d, lt='+string(lt_plot_ini),'TH teta 1d, lt='+string(lt_plot0),'TH teta 1d, lt='+string(lt_plot),'TH teta 1d, lt='+string(lt_plot2),'TH teta 1d, lt='+string(lt_plot3),'TH teta 1d, lt='+string(lt_plot4)]
1268	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1269	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta.ps'
1270	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1271	CS, SCALE=28
1272	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=1.4, TITLE=title_user
1273	cols=INDGEN(6)+2
1274	GPLOT, X=what_I_plot, Y=-alog(pplay_3d(*,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1275	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.2, YTITLE='-Log(P/P0) ',NDECS=1
1276
1277	oplot, teta_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1278	oplot, teta_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1279	oplot, teta_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1280	oplot, teta_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1281	oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1282	oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1283
1284	PSCLOSE, /NOVIEW
1285
1286	spawn, 'ps2png '+filename
1287	endif
1288
1289	endif else begin
1290
1291	print, '........ POTENTIAL TEMPERATURES'
1292
1293	what_I_plot = reform(teta_gcm(*,lt_plotindex_gcm))
1294	labels=['TH teta 1d, lt='+string(lt_plot)]
1295	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1296	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta.ps'
1297	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1298	CS, SCALE=28
1299	GSET, XMIN=210, XMAX=240, YMIN=0, YMAX=2, TITLE=title_user
1300	cols=INDGEN(1)+2
1301	GPLOT, X=what_I_plot, Y=-alog(pplay(*,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1302	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.2, YTITLE='-Log(P/P0) ',NDECS=1
1303
1304	oplot, teta_les(,lt_plotindex_les),-alog(pt(,lt_plotindex_les)/pGround), psym=4
1305	PSCLOSE, /NOVIEW
1306
1307	spawn, 'ps2png '+filename
1308	endelse
1309
1310	print, '........ SURFACE TEMPERATURES'
1311	what_I_plot = tsurf_gcm
1312	labels=['TH 1d tsurf']
1313	title_user = TestCase+SubCase+LayerCase+' Surface temperatures (recomputed from T and P)'
1314	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_tsurf.ps'
1315	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1316	CS, SCALE=28
1317	GSET, XMIN=6, XMAX=18, YMIN=200, YMAX=300, TITLE=title_user
1318	cols=INDGEN(1)+2
1319	GPLOT, X=localtime_gcm, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1320	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=5., YTITLE='Surface temperature',NDECS=1
1321
1322	oplot, localtime, tsurf, psym=4
1323	PSCLOSE, /NOVIEW
1324
1325	spawn, 'ps2png '+filename
1326
1327
1328	; ---------------------- * Vitesses verticales * -------------------------------
1329	; ------------ Verification de l'approx terrestre wmax = vmoy dans la couche instable
1330
1331	print, '........ CHECKING wmax = vmoy in unstable layer'
1332
1333	what_I_plot = uv_moy(*,lt_plotindex_les)
1334	labels=['LES uv_moy']
1335	title_user = TestCase+SubCase+LayerCase+' LES mean UV comp to max W in plume trac1'
1336	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_UV.ps'
1337	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1338	CS, SCALE=28
1339	GSET, XMIN=0, XMAX=10, YMIN=0, YMAX=10, TITLE=title_user
1340	cols=INDGEN(1)+2
1341	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1342	AXES, YSTEP = 1, YTITLE='Altitude (km)', XSTEP=1, XTITLE='Mean horizontal velocity inside domain (m/s)',NDECS=1
1343
1344	oplot, make_array(nz,value=wmax(lt_plotindex_les)), altitudes_LES/1000., psym=4
1345
1346	PSCLOSE, /NOVIEW
1347
1348	spawn, 'ps2png '+filename
1349
1350	; ------------ Profil de vitesse
1351
1352	print, '........ VERTICAL VELOCITY'
1353
1354	what_I_plot = make_array(nZmx)
1355	FOR k=0, nZmx-2 DO BEGIN
1356	what_I_plot(k) = 0.5*(zw2_lev(k,lt_plotindex_gcm) + zw2_lev(k+1,lt_plotindex_gcm))
1357	ENDFOR
1358	what_I_plot(nZmx-1) = 0.
1359
1360	smoothed_w_mean1_les = make_array(nz)
1361	smoothed_w_mean2_les = make_array(nz)
1362	smoothed_w_mean1_down_les = make_array(nz)
1363	FOR t=-ns,ns DO BEGIN
1364	smoothed_w_mean1_les = smoothed_w_mean1_les + REFORM(w_mean1(*,lt_plotindex_les+t))
1365	smoothed_w_mean2_les = smoothed_w_mean2_les + REFORM(w_mean2(*,lt_plotindex_les+t))
1366	smoothed_w_mean1_down_les = smoothed_w_mean1_down_les + REFORM(w_mean1_down(*,lt_plotindex_les+t))
1367	ENDFOR
1368
1369	smoothed_w_mean1_les = smoothed_w_mean1_les/nstot
1370	smoothed_w_mean2_les = smoothed_w_mean2_les/nstot
1371	smoothed_w_mean1_down_les = smoothed_w_mean1_down_les/nstot
1372
1373	ratio = make_array(nz)
1374	FOR k=0, nz-1 DO BEGIN
1375	IF(smoothed_w_mean1_les(k) ne 0.) then ratio(k) = smoothed_w_mean1_down_les(k)/smoothed_w_mean1_les(k) else ratio(k)=0.
1376	ENDFOR
1377
1378	labels=['TH 1d w, lt='+string(lt_plot)]
1379	title_user = TestCase+SubCase+LayerCase+' Vertical velocity comparisons (inside thermal)'
1380	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Wprofile.ps'
1381	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1382	CS, SCALE=28
1383	GSET, XMIN=-6, XMAX=8, YMIN=0, YMAX=10, TITLE=title_user
1384	cols=INDGEN(1)+2
1385	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1386	AXES, YSTEP = 2, YTITLE='Altitude (km)', XSTEP=1, XTITLE='Vertical velocity inside thermal (m/s)',NDECS=1
1387
1388	oplot, smoothed_w_mean1_les, altitudes_LES/1000., psym=4
1389	oplot, smoothed_w_mean2_les, altitudes_LES/1000., psym=5
1390	oplot, smoothed_w_mean1_down_les, altitudes_LES/1000., psym=6
1391
1392	PSCLOSE, /NOVIEW
1393
1394	spawn, 'ps2png '+filename
1395
1396
1397	; * Static stability *
1398
1399	print, '........ STATIC STABILITY'
1400
1401	dteta_dz_gcm = deriv(altitudes_GCM,reform(teta_gcm(*,lt_plotindex_gcm)))
1402	dteta_dz_les = deriv(altitudes_LES,reform(teta_les(*,lt_plotindex_les)))
1403
1404	what_I_plot = dteta_dz_gcm
1405	labels=['TH static stability 1d']
1406	title_user = TestCase+SubCase+LayerCase+' Static stability comparison'
1407	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_dTetadz.ps'
1408	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1409	CS, SCALE=28
1410	GSET, XMIN=-0.002, XMAX=0.006, YMIN=0, YMAX=10, TITLE=title_user
1411	cols=INDGEN(1)+2
1412	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1413	AXES, XSTEP = 0.0005 , XTITLE='Static stability (K.m-1)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
1414
1415	oplot, dteta_dz_les, altitudes_LES/1000., psym=4
1416
1417	PSCLOSE, /NOVIEW
1418
1419	spawn, 'ps2png '+filename
1420
1421	print,' -----------------------------------------------------------------------------------------------------------------------'
1422	print,' * LES diagnostics of the PLUME * MUAHAHAHAHAHA'
1423	print,' V2 with E and D computed for a UDE plume'
1424	print,' -----------------------------------------------------------------------------------------------------------------------'
1425
1426	print, '........ EXTRACTING DATA'
1427
1428	; --- Reinterpolation of F
1429
1430	fm_therm_gcm_interlay = make_array(nZmx,nTmx)
1431
1432	FOR k=0, nZmx-2 DO BEGIN
1433	fm_therm_gcm_interlay(k,) = (fm_therm_gcm_lev(k,) + fm_therm_gcm_lev(k+1,*))/2.
1434	ENDFOR
1435	fm_therm_gcm_interlay(nZmx-1,*)=0.
1436
1437	; --- Calculation of gcm df/dz using entrainment and detrainments and NOT F
1438
1439	df_dz_gcm = deriv(altitudes_GCM,reform(fm_therm_gcm_interlay(*,lt_plotindex_gcm)))
1440
1441	; --- Smoothing of the mass flux on a user-defined window
1442
1443	smoothed_fm_trac1_les = make_array(nz)
1444	smoothed_fm_trac2_les = make_array(nz)
1445	smoothed_downward_fm_trac1_les = make_array(nz)
1446	FOR t=-ns,ns DO BEGIN
1447	smoothed_fm_trac1_les = smoothed_fm_trac1_les + REFORM(fm_trac1_les(*,lt_plotindex_les+t))
1448	smoothed_fm_trac2_les = smoothed_fm_trac2_les + REFORM(fm_trac2_les(*,lt_plotindex_les+t))
1449	smoothed_downward_fm_trac1_les = smoothed_downward_fm_trac1_les + REFORM(downward_flux1(*,lt_plotindex_les+t))
1450	ENDFOR
1451
1452	smoothed_fm_trac1_les = smoothed_fm_trac1_les/nstot
1453	smoothed_fm_trac2_les = smoothed_fm_trac2_les/nstot
1454	smoothed_downward_fm_trac1_les = smoothed_downward_fm_trac1_les/nstot
1455
1456	; --- Calculation of the entrainement rate according to Rio(2010)
1457	; done in the heavy part at the begining (reeaaaally heavy)
1458
1459	; --- Smoothing of the entrainment on a ~30min window
1460
1461	; term 1
1462
1463
1464	smoothed_e_term1_trac1_les = make_array(nz)
1465	smoothed_e_term1_trac2_les = make_array(nz)
1466	FOR t=-ns,ns DO BEGIN
1467	smoothed_e_term1_trac1_les = smoothed_e_term1_trac1_les + REFORM(e_trac1_les(*,lt_plotindex_les+t))
1468	smoothed_e_term1_trac2_les = smoothed_e_term1_trac2_les + REFORM(e_trac2_les(*,lt_plotindex_les+t))
1469	ENDFOR
1470
1471	smoothed_e_term1_trac1_les = smoothed_e_term1_trac1_les/nstot
1472	smoothed_e_term1_trac2_les = smoothed_e_term1_trac2_les/nstot
1473
1474	smoothed_e_rate_term1_trac1_les = make_array(nz)
1475	smoothed_e_rate_trac2_les = smoothed_e_term1_trac2_les
1476
1477	; it already is an entrainment rate ! KIND OF : it is E/Mc, and Mc is not F !! NOW it is Mc/deltaTeta * dchi/dz
1478	FOR k=0, nz-1 DO BEGIN
1479	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_e_rate_term1_trac1_les(k) = smoothed_e_term1_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_e_rate_term1_trac1_les(k) =0.
1480	; IF (smoothed_fm_trac2_les(k) ne 0.) THEN smoothed_e_rate_trac2_les(k) = smoothed_e_trac2_les(k)/smoothed_fm_trac2_les(k) ELSE smoothed_e_rate_trac2_les(k) =0.
1481	ENDFOR
1482
1483	; term 2 & 3
1484
1485
1486	smoothed_e_term2_trac1_les = make_array(nz)
1487	smoothed_e_term3_trac1_les = make_array(nz)
1488	smoothed_e_term1_ude_trac1_les = make_array(nz)
1489	smoothed_e_term2_ude_trac1_les = make_array(nz)
1490	smoothed_e_term3_ude_trac1_les = make_array(nz)
1491
1492	FOR t=-ns,ns DO BEGIN
1493	smoothed_e_term2_trac1_les = smoothed_e_term2_trac1_les + REFORM(e1_term2(*,lt_plotindex_les+t))
1494	smoothed_e_term2_trac1_les = smoothed_e_term2_trac1_les + REFORM(e1_term2(*,lt_plotindex_les+t))
1495	smoothed_e_term1_ude_trac1_les = smoothed_e_term1_ude_trac1_les + REFORM(e1_term1_ude(*,lt_plotindex_les+t))
1496	smoothed_e_term2_ude_trac1_les = smoothed_e_term2_ude_trac1_les + REFORM(e1_term2_ude(*,lt_plotindex_les+t))
1497	smoothed_e_term3_ude_trac1_les = smoothed_e_term3_ude_trac1_les + REFORM(e1_term3_ude(*,lt_plotindex_les+t))
1498	ENDFOR
1499
1500	smoothed_e_term2_trac1_les = smoothed_e_term2_trac1_les/nstot
1501	smoothed_e_term3_trac1_les = smoothed_e_term3_trac1_les/nstot
1502	smoothed_e_term1_ude_trac1_les = smoothed_e_term1_ude_trac1_les/nstot
1503	smoothed_e_term2_ude_trac1_les = smoothed_e_term2_ude_trac1_les/nstot
1504	smoothed_e_term3_ude_trac1_les = smoothed_e_term3_ude_trac1_les/nstot
1505
1506	smoothed_e_rate_term2_trac1_les = make_array(nz)
1507	smoothed_e_rate_term3_trac1_les = make_array(nz)
1508	smoothed_e_rate_term1_ude_trac1_les = make_array(nz)
1509	smoothed_e_rate_term2_ude_trac1_les = make_array(nz)
1510	smoothed_e_rate_term3_ude_trac1_les = make_array(nz)
1511
1512	FOR k=0, nz-1 DO BEGIN
1513	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_e_rate_term2_trac1_les(k) = smoothed_e_term2_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_e_rate_term2_trac1_les(k) =0.
1514	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_e_rate_term3_trac1_les(k) = smoothed_e_term3_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_e_rate_term3_trac1_les(k) =0.
1515	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_e_rate_term1_ude_trac1_les(k) = smoothed_e_term1_ude_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_e_rate_term1_ude_trac1_les(k) =0.
1516	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_e_rate_term2_ude_trac1_les(k) = smoothed_e_term2_ude_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_e_rate_term2_ude_trac1_les(k) =0.
1517	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_e_rate_term3_ude_trac1_les(k) = smoothed_e_term3_ude_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_e_rate_term3_ude_trac1_les(k) =0.
1518
1519	ENDFOR
1520
1521	; --- Summing...
1522
1523	smoothed_e_rate_trac1_les = smoothed_e_rate_term1_trac1_les + smoothed_e_rate_term2_trac1_les + smoothed_e_rate_term3_trac1_les
1524	smoothed_e_rate_ude_trac1_les = smoothed_e_rate_term1_ude_trac1_les + smoothed_e_rate_term2_ude_trac1_les + smoothed_e_rate_term3_ude_trac1_les
1525
1526	;print, 'ommiting term3'
1527	;smoothed_e_rate_trac1_les = smoothed_e_rate_term1_trac1_les + smoothed_e_rate_term2_trac1_les
1528
1529	; --- Smoothing of the detrainment rate
1530
1531	smoothed_d_term1_trac1_les = make_array(nz)
1532	smoothed_d_term2_trac1_les = make_array(nz)
1533	smoothed_d_term3_trac1_les = make_array(nz)
1534	smoothed_d_term1_ude_trac1_les = make_array(nz)
1535	smoothed_d_term2_ude_trac1_les = make_array(nz)
1536	smoothed_d_term3_ude_trac1_les = make_array(nz)
1537
1538	FOR t=-ns,ns DO BEGIN
1539	smoothed_d_term1_trac1_les = smoothed_d_term1_trac1_les + REFORM(d1_term1(*,lt_plotindex_les+t))
1540	smoothed_d_term2_trac1_les = smoothed_d_term2_trac1_les + REFORM(d1_term2(*,lt_plotindex_les+t))
1541	smoothed_d_term3_trac1_les = smoothed_d_term3_trac1_les + REFORM(d1_term3(*,lt_plotindex_les+t))
1542	smoothed_d_term1_ude_trac1_les = smoothed_d_term1_ude_trac1_les + REFORM(d1_term1_ude(*,lt_plotindex_les+t))
1543	smoothed_d_term2_ude_trac1_les = smoothed_d_term2_ude_trac1_les + REFORM(d1_term2_ude(*,lt_plotindex_les+t))
1544	smoothed_d_term3_ude_trac1_les = smoothed_d_term3_ude_trac1_les + REFORM(d1_term3_ude(*,lt_plotindex_les+t))
1545	ENDFOR
1546
1547	smoothed_d_term1_trac1_les = smoothed_d_term1_trac1_les/nstot
1548	smoothed_d_term2_trac1_les = smoothed_d_term2_trac1_les/nstot
1549	smoothed_d_term3_trac1_les = smoothed_d_term3_trac1_les/nstot
1550	smoothed_d_term1_ude_trac1_les = smoothed_d_term1_ude_trac1_les/nstot
1551	smoothed_d_term2_ude_trac1_les = smoothed_d_term2_ude_trac1_les/nstot
1552	smoothed_d_term3_ude_trac1_les = smoothed_d_term3_ude_trac1_les/nstot
1553
1554	smoothed_d_rate_term1_trac1_les = make_array(nz)
1555	smoothed_d_rate_term2_trac1_les = make_array(nz)
1556	smoothed_d_rate_term3_trac1_les = make_array(nz)
1557	smoothed_d_rate_term1_ude_trac1_les = make_array(nz)
1558	smoothed_d_rate_term2_ude_trac1_les = make_array(nz)
1559	smoothed_d_rate_term3_ude_trac1_les = make_array(nz)
1560
1561	FOR k=0, nz-1 DO BEGIN
1562	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_d_rate_term1_trac1_les(k) = smoothed_d_term1_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_d_rate_term1_trac1_les(k) =0.
1563	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_d_rate_term2_trac1_les(k) = smoothed_d_term2_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_d_rate_term2_trac1_les(k) =0.
1564	IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_d_rate_term3_trac1_les(k) = smoothed_d_term3_trac1_les(k)/smoothed_fm_trac1_les(k) ELSE smoothed_d_rate_term3_trac1_les(k) =0.
1565	IF (smoothed_fm_trac1_les(k) ne 0.) THEN BEGIN
1566	smoothed_d_rate_term1_ude_trac1_les(k) = smoothed_d_term1_ude_trac1_les(k)/smoothed_fm_trac1_les(k)
1567	smoothed_d_rate_term2_ude_trac1_les(k) = smoothed_d_term2_ude_trac1_les(k)/smoothed_fm_trac1_les(k)
1568	smoothed_d_rate_term3_ude_trac1_les(k) = smoothed_d_term3_ude_trac1_les(k)/smoothed_fm_trac1_les(k)
1569	ENDIF ELSE BEGIN
1570	smoothed_d_rate_term1_ude_trac1_les(k)=0.
1571	smoothed_d_rate_term2_ude_trac1_les(k)=0.
1572	smoothed_d_rate_term3_ude_trac1_les(k)=0.
1573	ENDELSE
1574	ENDFOR
1575
1576	; --- Summing...
1577
1578	smoothed_d_rate_trac1_les = smoothed_d_rate_term1_trac1_les+smoothed_d_rate_term2_trac1_les+smoothed_d_rate_term3_trac1_les
1579	smoothed_d_rate_ude_trac1_les = smoothed_d_rate_term1_ude_trac1_les+smoothed_d_rate_term2_ude_trac1_les+smoothed_d_rate_term3_ude_trac1_les
1580	;print, 'ommiting term3'
1581	;smoothed_d_rate_trac1_les = smoothed_d_rate_term1_trac1_les+smoothed_d_rate_term2_trac1_les
1582
1583	; --- PLOTTING : BUOYANCY TERM
1584
1585	smoothed_buoyancy_trac1_les = make_array(nz)
1586	smoothed_buoyancy_ude_trac1_les = make_array(nz)
1587	smoothed_buoyancy_trac2_les = make_array(nz)
1588	smoothed_buoyancy_downdraft1_les_ude = make_array(nz)
1589
1590	FOR t=-ns,ns DO BEGIN
1591	smoothed_buoyancy_trac1_les = smoothed_buoyancy_trac1_les + REFORM(buoyancy1_les(*,lt_plotindex_les+t))
1592	smoothed_buoyancy_ude_trac1_les = smoothed_buoyancy_ude_trac1_les + REFORM(grav(tplume1moy(,lt_plotindex_les+t)/tenv1moy_ude(*,lt_plotindex_les+t)-1.))
1593	smoothed_buoyancy_trac2_les = smoothed_buoyancy_trac2_les + REFORM(buoyancy2_les(*,lt_plotindex_les+t))
1594	smoothed_buoyancy_downdraft1_les_ude = smoothed_buoyancy_downdraft1_les_ude + REFORM(grav(tdown1moy(,lt_plotindex_les+t)/tenv1moy_ude(*,lt_plotindex_les+t)-1.))
1595	ENDFOR
1596
1597	smoothed_buoyancy_trac1_les = smoothed_buoyancy_trac1_les/nstot
1598	smoothed_buoyancy_ude_trac1_les = smoothed_buoyancy_ude_trac1_les/nstot
1599	smoothed_buoyancy_trac2_les = smoothed_buoyancy_trac2_les/nstot
1600	smoothed_buoyancy_downdraft1_les_ude = smoothed_buoyancy_downdraft1_les_ude/nstot
1601
1602	print, '........ BUOYANCY'
1603
1604	what_I_plot = [[buoyancy_gcm(,lt_plotindex_gcm)],[buoyancy_est_gcm(,lt_plotindex_gcm)]]
1605	labels=['TH buoyancy term','TH estimated buoyancy in plume']
1606	title_user = TestCase+SubCase+LayerCase+' UDE plume buoyancy'
1607	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_B.ps'
1608	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1609	CS, SCALE=28
1610	GSET, XMIN=-0.06, XMAX=0.06, YMIN=0, YMAX=10, TITLE=title_user
1611	cols=INDGEN(2)+2
1612	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1613	AXES, XSTEP = 0.01 , XTITLE='N.m-1', YSTEP=1, YTITLE=' Altitude (km)',NDECS=3
1614
1615	;oplot, smoothed_buoyancy_trac1_les, altitudes_LES/1000., psym=4
1616	;oplot, smoothed_buoyancy_trac2_les, altitudes_LES/1000., psym=5
1617	print, smoothed_buoyancy_ude_trac1_les
1618	oplot, smoothed_buoyancy_ude_trac1_les, altitudes_LES/1000., psym=4
1619	oplot, smoothed_buoyancy_downdraft1_les_ude, altitudes_LES/1000., psym=7
1620
1621	PSCLOSE, /NOVIEW
1622
1623	spawn, 'ps2png '+filename
1624
1625
1626	; --- PLOTTING : MASS FLUX
1627
1628	print, '........ MASS FLUX'
1629
1630	f_gcm = fm_therm_gcm_interlay(*,lt_plotindex_gcm)
1631	what_I_plot = f_gcm
1632	labels=['TH mass flux']
1633	title_user = TestCase+SubCase+LayerCase+' mass flux comparison, average over '+taverage+' mn'
1634	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_f.ps'
1635	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1636	CS, SCALE=28
1637	GSET, XMIN=-0.008, XMAX=0.008, YMIN=0, YMAX=10, TITLE=title_user
1638	cols=INDGEN(1)+2
1639	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1640	AXES, XSTEP = 0.002 , XTITLE='Kg.m-2.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
1641
1642	oplot, smoothed_fm_trac1_les, altitudes_LES/1000., psym=4
1643	oplot, smoothed_fm_trac2_les, altitudes_LES/1000., psym=5
1644	oplot, smoothed_downward_fm_trac1_les, altitudes_LES/1000., psym=6
1645
1646	PSCLOSE, /NOVIEW
1647
1648	spawn, 'ps2png '+filename
1649
1650
1651	; --- PLOTTING : MASS FLUX DERIVATIVE
1652
1653	print, '........ MASS FLUX DERIVATIVE'
1654
1655	df_dz_les1 = deriv(altitudes_LES,reform(smoothed_fm_trac1_les))
1656	df_dz_les2 = deriv(altitudes_LES,reform(smoothed_fm_trac2_les))
1657
1658	what_I_plot = df_dz_gcm
1659	labels=['TH mass flux vertical derivative']
1660	title_user = TestCase+SubCase+LayerCase+' mass flux vertical derivative comparison, average over '+taverage+' mn'
1661	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_dfdz.ps'
1662	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1663	CS, SCALE=28
1664	GSET, XMIN=-0.00002, XMAX=0.00002, YMIN=0, YMAX=10, TITLE=title_user
1665	cols=INDGEN(1)+2
1666	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1667	AXES, XSTEP = 0.000005 , XTITLE='Kg.m-3.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=6
1668	oplot, df_dz_les1, altitudes_LES/1000., psym=4
1669	oplot, df_dz_les2, altitudes_LES/1000., psym=5
1670
1671	PSCLOSE, /NOVIEW
1672
1673	spawn, 'ps2png '+filename
1674
1675	; --- PLOTTING : ENTRAINMENT RATE e = E/f
1676
1677	print, '........ ENTRAINMENT RATE'
1678
1679	e_gcm = make_array(nZmx)
1680
1681	FOR k=0, nZmx-1 DO BEGIN
1682	IF (fm_therm_gcm_interlay(k,lt_plotindex_gcm) ne 0.) THEN BEGIN
1683	e_gcm(k) = zdz_entr_therm_gcm(k,lt_plotindex_gcm)/(approx_zdz_gcm(k)*fm_therm_gcm_interlay(k,lt_plotindex_gcm))
1684	ENDIF ELSE BEGIN
1685	e_gcm(k) = 0.
1686	ENDELSE
1687	ENDFOR
1688
1689
1690	what_I_plot = e_gcm
1691	labels=['TH entrainment rate']
1692	title_user = TestCase+SubCase+LayerCase+' UDE entrainment rate comparison, average over '+taverage+' mn'
1693	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e.ps'
1694	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1695	CS, SCALE=28
1696	GSET, XMIN=-0.003, XMAX=0.003, YMIN=0, YMAX=10, TITLE=title_user
1697	cols=INDGEN(1)+2
1698	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1699	AXES, XSTEP = 0.0006 , XTITLE='entrainment rate m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
1700
1701	oplot, smoothed_e_rate_ude_trac1_les, altitudes_LES/1000., psym=4
1702	;oplot, smoothed_e_rate_trac1_les, altitudes_LES/1000., psym=4
1703	;oplot, smoothed_e_rate_trac2_les, altitudes_LES/1000., psym=5
1704
1705	PSCLOSE, /NOVIEW
1706
1707	spawn, 'ps2png '+filename
1708
1709	; --- PLOTTING : EXTENDED ENTRAINMENT RATE
1710
1711	print, '........ EXTENDED ENTRAINMENT RATE'
1712
1713	;what_I_plot = [[smoothed_e_rate_term1_trac1_les],[smoothed_e_rate_term2_trac1_les],[smoothed_e_rate_term3_trac1_les],[smoothed_e_rate_trac1_les]]
1714	what_I_plot = [[smoothed_e_rate_term1_ude_trac1_les],[smoothed_e_rate_term2_ude_trac1_les],[smoothed_e_rate_term3_ude_trac1_les],[smoothed_e_rate_ude_trac1_les]]
1715	labels=['LES base entrainment rate','LES term2 e rate','LES term3 e rate','LES total e rate']
1716	title_user = TestCase+SubCase+LayerCase+' UDE entrainment rate comparison, average over '+taverage+' mn'
1717	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e_terms.ps'
1718	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1719	CS, SCALE=28
1720	GSET, XMIN=-0.01, XMAX=0.01, YMIN=0, YMAX=10, TITLE=title_user
1721	cols=INDGEN(4)+2
1722	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1723	AXES, XSTEP = 0.005 , XTITLE='entrainment rate m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=3
1724
1725	PSCLOSE, /NOVIEW
1726
1727	spawn, 'ps2png '+filename
1728
1729
1730	; --- PLOTTING : EXTENDED DETRAINMENT RATE
1731
1732	print, '........ EXTENDED DETRAINMENT RATE'
1733
1734	what_I_plot = [[smoothed_d_rate_term1_ude_trac1_les],[smoothed_d_rate_term2_ude_trac1_les],[smoothed_d_rate_term3_ude_trac1_les],[smoothed_d_rate_ude_trac1_les]]
1735	labels=['LES term 1 detrainment rate','LES term2 d rate','LES term3 d rate','LES Total d rate']
1736	title_user = TestCase+SubCase+LayerCase+' UDE detrainment rate comparison, average over '+taverage+' mn'
1737	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d_terms.ps'
1738	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1739	CS, SCALE=28
1740	GSET, XMIN=-0.01, XMAX=0.01, YMIN=0, YMAX=10, TITLE=title_user
1741	cols=INDGEN(4)+2
1742	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1743	AXES, XSTEP = 0.002 , XTITLE='detrainment rate m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
1744
1745	PSCLOSE, /NOVIEW
1746
1747	spawn, 'ps2png '+filename
1748
1749	; --- PLOTTING : DETRAINMENT RATE d = D/f
1750
1751	print, '........ DETRAINMENT RATE'
1752
1753	;smoothed_d_rate_trac2_les = make_array(nz)
1754	;
1755	;FOR k=0, nz-1 DO BEGIN
1756	; IF (smoothed_fm_trac1_les(k) ne 0.) THEN smoothed_d_rate_trac1_les(k) = smoothed_e_rate_trac1_les(k) - df_dz_les1(k)/smoothed_fm_trac1_les(k) ELSE smoothed_d_rate_trac1_les(k) =0.
1757	; IF (smoothed_fm_trac2_les(k) ne 0.) THEN smoothed_d_rate_trac2_les(k) = smoothed_e_rate_trac2_les(k) - df_dz_les2(k)/smoothed_fm_trac2_les(k) ELSE smoothed_d_rate_trac2_les(k) =0.
1758	;ENDFOR
1759	;
1760	d_gcm = make_array(nZmx)
1761	FOR k=0, nZmx-1 DO BEGIN
1762	IF (fm_therm_gcm_interlay(k,lt_plotindex_gcm) ne 0.) THEN BEGIN
1763	d_gcm(k) = zdz_detr_therm_gcm(k,lt_plotindex_gcm)/(approx_zdz_gcm(k)*fm_therm_gcm_interlay(k,lt_plotindex_gcm))
1764	ENDIF ELSE BEGIN
1765	d_gcm(k) = 0.
1766	ENDELSE
1767	ENDFOR
1768
1769	what_I_plot = d_gcm
1770	labels=['TH detrainment rate']
1771	title_user = TestCase+SubCase+LayerCase+' UDE detrainment rate comparison, average over '+taverage+' mn'
1772	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d.ps'
1773	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1774	CS, SCALE=28
1775	GSET, XMIN=0., XMAX=0.03, YMIN=0, YMAX=10, TITLE=title_user
1776	cols=INDGEN(1)+2
1777	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1778	AXES, XSTEP = 0.005 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=2
1779
1780	oplot, smoothed_d_rate_ude_trac1_les, altitudes_LES/1000., psym=4
1781	;oplot, smoothed_d_rate_trac1_les, altitudes_LES/1000., psym=4
1782	;oplot, smoothed_d_rate_trac2_les, altitudes_LES/1000., psym=5
1783
1784	PSCLOSE, /NOVIEW
1785	spawn, 'ps2png '+filename
1786
1787	; --- PLOTTING : FRACTION COVERAGE
1788
1789	print, '........ EXTENDED ALPHA'
1790
1791	what_I_plot = alpha_interlay_gcm
1792	labels=['TH alpha']
1793	title_user = TestCase+SubCase+LayerCase+' fraction coverage comparison, average over '+taverage+' mn'
1794	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_alpha.ps'
1795	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1796	CS, SCALE=28
1797	GSET, XMIN=0., XMAX=1., YMIN=0, YMAX=10, TITLE=title_user
1798	cols=INDGEN(1)+2
1799	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1800	AXES, XSTEP = 0.1 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
1801
1802	oplot, smoothed_alpha1_les, altitudes_LES/1000., psym=4
1803	oplot, smoothed_alpha2_les, altitudes_LES/1000., psym=5
1804	oplot, smoothed_beta1_les, altitudes_LES/1000., psym=6
1805
1806	PSCLOSE, /NOVIEW
1807	spawn, 'ps2png '+filename
1808
1809	; --- PLOTTING : THEORETICAL ENTRAINMENT RATE FROM LES DATA
1810
1811	print, '........ PARAMETRIZED RATES'
1812
1813	approx_zdz_les = make_array(nz)
1814
1815	approx_zdz_les(0)=altitudes_LES(1)
1816	FOR k=1, nz-2 DO BEGIN
1817	approx_zdz_les(k) = altitudes_LES(k+1) - altitudes_LES(k)
1818	ENDFOR
1819	approx_zdz_les(nz-1)=approx_zdz_les(nz-2)
1820
1821
1822	theoretical_e_trac1_les = make_array(nz)
1823	theoretical_e_trac2_les = make_array(nz)
1824
1825
1826	FOR k=0, nz-1 DO BEGIN
1827	theoretical_e_trac1_les(k) = MAX([0.,(betalpha/(1.+betalpha))((afactsmoothed_buoyancy_trac1_les(k)/((smoothed_w_mean1_les(k))^2.)) - fact_epsilon)])
1828	theoretical_e_trac2_les(k) = MAX([0.,(betalpha/(1.+betalpha))((afactsmoothed_buoyancy_trac2_les(k)/((smoothed_w_mean2_les(k))^2.)) - fact_epsilon)])
1829	ENDFOR
1830
1831
1832	what_I_plot = [[theoretical_e_trac1_les],[theoretical_e_trac2_les]]
1833	labels=['LES TH theo e rate trac1','LES TH theo e rate trac2']
1834	title_user = TestCase+SubCase+LayerCase+' comp. theor. entr. rate comparison, average over '+taverage+' mn'
1835	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e_theoretical.ps'
1836	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1837	CS, SCALE=28
1838	GSET, XMIN=-0.015, XMAX=0.03, YMIN=0, YMAX=10, TITLE=title_user
1839	cols=INDGEN(2)+2
1840	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1841	AXES, XSTEP = 0.003 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=3
1842
1843	oplot, smoothed_e_rate_trac1_les, altitudes_LES/1000., psym=4
1844	oplot, smoothed_e_rate_trac2_les, altitudes_LES/1000., psym=5
1845
1846	PSCLOSE, /NOVIEW
1847
1848	spawn, 'ps2png '+filename
1849
1850	; --- PLOTTING : THEORETICAL DETRAINMENT RATE FROM LES DATA
1851	; ZDZ STUFF REMOVED
1852
1853	print, '........ PARAMETRIZED DETRAINMENT'
1854
1855	theoretical_d_trac1_les = make_array(nz)
1856	theoretical_d_trac2_les = make_array(nz)
1857
1858	FOR k=0, nz-1 DO BEGIN
1859	theoretical_d_trac1_les(k) = MAX([detr_min,-afact(betalpha/(1.+betalpha))(smoothed_buoyancy_trac1_les(k)/((smoothed_w_mean1_les(k))^2.))])
1860	theoretical_d_trac2_les(k) = MAX([detr_min,-afact(betalpha/(1.+betalpha))(smoothed_buoyancy_trac2_les(k)/((smoothed_w_mean2_les(k))^2.))])
1861	ENDFOR
1862
1863	what_I_plot = [[theoretical_d_trac1_les],[theoretical_d_trac2_les]]
1864	labels=['LES TH theo d rate trac1','LES TH theo d rate trac2']
1865	title_user = TestCase+SubCase+LayerCase+' comp. theor. detr. rate comparison, average over '+taverage+' mn'
1866	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d_theoretical.ps'
1867	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1868	CS, SCALE=28
1869	GSET, XMIN=-0.1, XMAX=0.1, YMIN=0, YMAX=10, TITLE=title_user
1870	cols=INDGEN(2)+2
1871	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1872	AXES, XSTEP = 0.01 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=2
1873
1874	oplot, smoothed_d_rate_trac1_les, altitudes_LES/1000., psym=4
1875	;oplot, smoothed_d_rate_trac2_les, altitudes_LES/1000., psym=5
1876
1877	PSCLOSE, /NOVIEW
1878
1879	spawn, 'ps2png '+filename
1880
1881	; --- PLOTTING : THEORETICAL E-D RATE FROM LES DATA
1882
1883	print, '........ PARAMETRIZED MASS FLUX DERIVATIVE'
1884
1885	theoretical_dfdz_f_trac1_les = make_array(nz)
1886	theoretical_dfdz_f_trac2_les = make_array(nz)
1887
1888	theoretical_dfdz_f_trac1_les = theoretical_e_trac1_les - theoretical_d_trac1_les
1889	theoretical_dfdz_f_trac2_les = theoretical_e_trac2_les - theoretical_d_trac2_les
1890
1891	df_dz_f_les1 = make_array(nz)
1892	df_dz_f_les2 = make_array(nz)
1893
1894	FOR k=0, nz-1 DO BEGIN
1895	IF (smoothed_fm_trac1_les(k) ne 0.) THEN df_dz_f_les1(k) = df_dz_les1(k)/smoothed_fm_trac1_les(k) ELSE df_dz_f_les1(k)=0.
1896	IF (smoothed_fm_trac2_les(k) ne 0.) THEN df_dz_f_les2(k) = df_dz_les2(k)/smoothed_fm_trac2_les(k) ELSE df_dz_f_les2(k)=0.
1897	ENDFOR
1898
1899	what_I_plot = [[theoretical_dfdz_f_trac1_les],[theoretical_dfdz_f_trac2_les]]
1900	labels=['LES TH theo 1/f df/dz trac1','LES TH theo 1/f df/dz trac2']
1901	title_user = TestCase+SubCase+LayerCase+' comp. theor. entr. - detr. rate comparison, average over '+taverage+' mn'
1902	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_dfdzf_theoretical.ps'
1903	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1904	CS, SCALE=28
1905	GSET, XMIN=-0.02, XMAX=0.02, YMIN=0, YMAX=10, TITLE=title_user
1906	cols=INDGEN(2)+2
1907	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1908	AXES, XSTEP = 0.01 , XTITLE='entr - detr (rates) m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
1909
1910	oplot, df_dz_f_les1, altitudes_LES/1000., psym=4
1911	oplot, df_dz_f_les2, altitudes_LES/1000., psym=5
1912
1913	PSCLOSE, /NOVIEW
1914
1915	spawn, 'ps2png '+filename
1916
1917	; --- PLOTTING : e versus B/w2
1918
1919	print, '........ BUOYANCY AND VERTICAL VELOCITY ENTRAINMENT RATE DEPENDENCY'
1920
1921	B_w2_trac1 = make_array(nz)
1922	B_w2_trac2 = make_array(nz)
1923	;dwdz_trac1 = deriv(altitudes_LES,smoothed_w_mean1_les)
1924	;dwdz_trac2 = deriv(altitudes_LES,smoothed_w_mean2_les)
1925	;full_dwdz_trac1 = make_array(nz,nttot)
1926	;full_dadz_trac1 = make_array(nz,nttot)
1927	;FOR l=0, nttot -1 DO BEGIN
1928	; full_dwdz_trac1(,l) = deriv(altitudes_LES,w_mean1(,l))
1929	; full_dadz_trac1(,l) = deriv(altitudes_LES,alpha1out(,l))
1930	;ENDFOR
1931	;alpha = 0.
1932
1933	FOR k=0, nz-1 DO BEGIN
1934	IF (smoothed_e_rate_trac1_les(k) ne 0.) THEN B_w2_trac1(k) = smoothed_buoyancy_ude_trac1_les(k)/(smoothed_w_mean1_les(k))^2 ELSE B_w2_trac1(k)=-0.
1935	; IF (smoothed_e_rate_trac2_les(k) ne 0.) THEN B_w2_trac2(k) = smoothed_buoyancy_trac2_les(k)/(smoothed_w_mean2_les(k))^2 ELSE B_w2_trac2(k)=-0.
1936	ENDFOR
1937
1938	;FOR zzz=0.,30 DO BEGIN
1939
1940	;alpha = zzz/10.
1941	;
1942	;FOR k=0, nz-1 DO BEGIN
1943	; IF (smoothed_e_rate_trac1_les(k) ne 0. and smoothed_w_mean1_les(k) ne 0.) THEN B_w2_trac1(k) = 0.5(smoothed_buoyancy_trac1_les(k)/(smoothed_w_mean1_les(k))^2 - alpha(1./smoothed_w_mean1_les(k))*dwdz_trac1(k)) ELSE B_w2_trac1(k)=0.
1944	; IF (smoothed_e_rate_trac2_les(k) ne 0. and smoothed_w_mean2_les(k) ne 0.) THEN B_w2_trac2(k) = 0.5(smoothed_buoyancy_trac2_les(k)/(smoothed_w_mean2_les(k))^2 - alpha(1./smoothed_w_mean2_les(k))*dwdz_trac2(k)) ELSE B_w2_trac2(k)=0.
1945	;ENDFOR
1946
1947	;print, smoothed_buoyancy_trac1_les()/(smoothed_w_mean1_les())^2
1948	;print, (1./smoothed_w_mean1_les())dwdz_trac1(*)
1949
1950	full_e1 = make_array(nz,nttot)
1951	full_bw2 = make_array(nz,nttot)
1952	FOR k=0, nz-1 DO BEGIN
1953	FOR l=0, nttot-1 DO BEGIN
1954	if(fm_trac1_les(k,l) ne 0.) then full_e1(k,l)=(e1_term1_ude(k,l)+e1_term2_ude(k,l)+e1_term3_ude(k,l))/fm_trac1_les(k,l) else full_e1(k,l)=0.
1955	if(w_mean1(k,l) ne 0.) then full_bw2(k,l)=grav*(tplume1moy(k,l)/tenv1moy_ude(k,l) -1.)/(w_mean1(k,l)^2) else full_bw2(k,l)=-0.
1956	; if(w_mean1(k,l) ne 0.) then full_bw2(k,l)=0.5(alphabuoyancy1_les(k,l)/(w_mean1(k,l)^2) - full_dwdz_trac1(k,l)/w_mean1(k,l)) else full_bw2(k,l)=0.
1957	ENDFOR
1958	ENDFOR
1959
1960	what_I_plot = smoothed_e_rate_ude_trac1_les
1961	labels=['e_rate trac1']
1962	title_user = TestCase+SubCase+LayerCase+' LES UDE entrainment rate dep with B/w2, average over '+taverage+' mn,'
1963	;filename = TestCase+SubCase+LayerCase+string(alpha,format='(F3.1)')+'Gcm_Les_Comp_e_Bw2.ps'
1964	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e_Bw2.ps'
1965	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1966	CS, SCALE=28
1967	GSET, XMIN=0., XMAX=0.4, YMIN=0., YMAX=0.4, TITLE=title_user
1968	cols=INDGEN(1)+2
1969	GPLOT, X=what_I_plot, Y=B_w2_trac1, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30, SYM=5, /NOLINES
1970	AXES, XSTEP = 0.05 , XTITLE='LES Entrainment rate m-1', YSTEP=0.05, YTITLE='Parametrized entrainement rate m-1',NDECS=4
1971
1972	;oplot, smoothed_e_rate_trac2_les, B_w2_trac2, psym=5
1973	FOR l=0, nttot-1 DO BEGIN
1974	oplot, full_e1(,l),full_bw2(,l),thick=0.05,psym=1
1975	ENDFOR
1976	;mean_full_e1 = make_array(nz) & mean_full_bw2 = make_array(nz)
1977	;FOR k=0, nz-1 DO BEGIN
1978	; mean_full_e1(k) = MEAN(reform(full_e1(k,*)))
1979	; mean_full_bw2(k) = MEAN(reform(full_bw2(k,*)))
1980	;ENDFOR
1981	;oplot, mean_full_e1, mean_full_bw2, thick=0.3, psym = 2,color=5
1982	;oplot, theoretical_e_trac1_les, B_w2_trac1,psym=2,thick=0.8,color=7
1983	oplot,(B_w2_trac1)/2.2222 + 0.0005,B_w2_trac1,thick=0.3,color=7
1984	;oplot, 0.0118*(B_w2_trac1/0.043)^(1./1.65),B_w2_trac1,thick=0.3,color=7
1985	oplot, 0.0113*(B_w2_trac1/0.043)^(1./1.65),B_w2_trac1,thick=0.3,color=7
1986	oplot, 0.0105*(B_w2_trac1/0.048)^(1./1.7),B_w2_trac1,thick=0.3,color=6
1987	;oplot, alog((B_w2_trac1 - 0.000942361)/0.0444855) - 3.85453, B_w2_trac1, thick=0.3,color=7
1988
1989	;print, alog((B_w2_trac1)/0.0444855) - 3.85453
1990
1991	PSCLOSE, /NOVIEW
1992
1993	spawn, 'ps2png '+filename
1994
1995	;ENDFOR
1996
1997	what_I_plot = full_bw2(*,lt_plotindex_les)
1998	labels=['B/w2']
1999	title_user = TestCase+SubCase+LayerCase+' LES UDE B/w2, average over '+taverage+' mn,'
2000	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Bw2.ps'
2001	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2002	CS, SCALE=28
2003	GSET, XMIN=-0.01, XMAX=0.01, YMIN=0., YMAX=6., TITLE=title_user
2004	cols=INDGEN(1)+2
2005	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2006	AXES, XSTEP = 0.002 , XTITLE='B/w2 term in LES UDE', YSTEP=0.5, YTITLE='Altitude (km)',NDECS=4
2007
2008	;FOR l=0, nttot-1 DO BEGIN
2009	; oplot, full_bw2(*,l),altitudes_LES/1000.,thick=0.05,psym=1
2010	;ENDFOR
2011
2012	PSCLOSE, /NOVIEW
2013
2014	spawn, 'ps2png '+filename
2015
2016
2017	print, '........ BUOYANCY AND VERTICAL VELOCITY DETRAINMENT RATE DEPENDENCY'
2018
2019	full_d1 = make_array(nz,nttot)
2020	;full_dSiebesma = make_array(nz,nttot)
2021
2022
2023	FOR k=0, nz-1 DO BEGIN
2024	FOR l=0, nttot-1 DO BEGIN
2025	if(fm_trac1_les(k,l) ne 0.) then full_d1(k,l)=(d1_term1_ude(k,l)+d1_term2_ude(k,l)+d1_term3_ude(k,l))/fm_trac1_les(k,l) else full_d1(k,l)=-0.
2026	; if(w_mean1(k,l) ne 0.) then full_dSiebesma(k,l)=0.5buoyancy1_les(k,l)/(w_mean1(k,l)^2) -1.5full_dwdz_trac1(k,l)/w_mean1(k,l) - full_dadz_trac1(k,l)/alpha1out(k,l) else full_dSiebesma(k,l)=-0.
2027	ENDFOR
2028	ENDFOR
2029
2030	what_I_plot = smoothed_d_rate_ude_trac1_les
2031	labels=['d_rate trac1']
2032	title_user = TestCase+SubCase+LayerCase+' LES UDE detrainment rate dep with B/w2, average over '+taverage+' mn'
2033	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d_Bw2.ps'
2034	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2035	CS, SCALE=28
2036	GSET, XMIN=-0, XMAX=0.4, YMIN=0., YMAX=0.4, TITLE=title_user
2037	cols=INDGEN(1)+2
2038	GPLOT, X=what_I_plot, Y=full_bw2(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30, SYM=5, /NOLINES
2039	AXES, XSTEP = 0.05 , XTITLE='Detrainment rate m-1', YSTEP=0.05, YTITLE='B/w²',NDECS=4
2040
2041	FOR l=0, nttot-1 DO BEGIN
2042	oplot, full_d1(,l),full_bw2(,l),thick=0.05,psym=1
2043	ENDFOR
2044	;oplot, theoretical_d_trac1_les, full_bw2(*,lt_plotindex_les),psym=2,thick=0.8,color=7
2045	;oplot,B_w2_trac1/2.7 + 0.0002,B_w2_trac1,thick=0.3,color=7
2046	oplot,B_w2_trac1/2.222 + 0.0002,B_w2_trac1,thick=0.3,color=7
2047	oplot, 0.0105*(B_w2_trac1/0.048)^(1./1.7),B_w2_trac1,thick=0.3,color=7
2048	oplot, 0.0118*(B_w2_trac1/0.043)^(1./1.65),B_w2_trac1,thick=0.3,color=6
2049
2050	PSCLOSE, /NOVIEW
2051
2052	spawn, 'ps2png '+filename
2053
2054
2055	; --- PLOTTING : 0.5*dwu2/dz
2056	print, '........ ///////////// starting local thermal model ///////////'
2057
2058	print, ' -> alimentation'
2059	a_star = make_array(nz, value=0.)
2060	f_star = make_array(nz, value=0.)
2061	f_star(0) = 1.
2062	teta_est = make_array(nz, value=0.)
2063	teta_p = make_array(nz, value=0.)
2064	zw2 = make_array(nz,value=0.)
2065	w_est = make_array(nz,value=0.)
2066	entr_star = make_array(nz,value=0.)
2067	detr_star = make_array(nz,value=0.)
2068	zbuoy_est = make_array(nz,value=0.)
2069	zbuoy = make_array(nz,value=0.)
2070	a_star_tot = 0.
2071
2072	zw2(1)= 0.40.38115522grav(teta_les(0,lt_plotindex_les)/teta_les(1,lt_plotindex_les) -1.)*approx_zdz_les(0)
2073	w_est(1) = zw2(1)
2074	FOR k=0, nz-2 DO BEGIN
2075	if ((teta_les(k,lt_plotindex_les) GT (teta_les(k+1,lt_plotindex_les) +0.1)) AND (teta_les(0,lt_plotindex_les) GE teta_les(k,lt_plotindex_les))) then begin
2076	a_star(k) = MAX([(teta_les(k,lt_plotindex_les)-teta_les(k+1,lt_plotindex_les)),0.])*sqrt(altitudes_LES(k))
2077	a_star_tot = a_star_tot + a_star(k)
2078	lalim = k+1.
2079	endif
2080	ENDFOR
2081	FOR k=0, nz-1 DO BEGIN
2082	a_star(k) = a_star(k)/a_star_tot
2083	ENDFOR
2084	print, 'alimentation :'
2085	;print, a_star
2086	f_star(0)=0.
2087	f_star(1) = a_star(0)
2088	teta_p(0) = teta_les(0,lt_plotindex_les)
2089	teta_est(0) = teta_les(0,lt_plotindex_les)
2090	print, ' -> plume'
2091	FOR k=1, nz-2 DO BEGIN
2092	if (k LT lalim) then begin
2093	teta_est(k) = (f_star(k)teta_p(k-1)+a_star(k)0.5*(teta_les(k,lt_plotindex_les) + teta_p(k-1)))/(f_star(k) + a_star(k))
2094	endif else begin
2095	teta_est(k) = teta_p(k-1)
2096	endelse
2097	zbuoy_est(k) = grav*(teta_est(k)/teta_les(k,lt_plotindex_les) -1.)
2098	zw2fact=fact_epsilon2.approx_zdz_les(k)/(1.+betalpha)
2099	zdw2=afact*zbuoy_est(k)/fact_epsilon
2100	w_est(k+1) = MAX([0.0001,exp(-zw2fact)*(w_est(k)-zdw2)+zdw2])
2101	if (w_est(k+1) lt 0.) then begin
2102	w_est(k+1)=zw2(k)
2103	endif
2104	if (w_est(k+1) gt 0.001) then begin
2105	entr_star(k)=f_star(k)approx_zdz_les(k)(betalpha/(1.+betalpha))MAX([0.,afactzbuoy_est(k)/w_est(k+1)])
2106	detr_star(k)=f_star(k)approx_zdz_les(k)MAX([detr_min,-afact(betalpha/(1.+betalpha))zbuoy_est(k)/w_est(k+1)])
2107	endif
2108	if (k lt lalim) then begin
2109	a_star(k)=max([a_star(k),entr_star(k)])
2110	entr_star(k)=0.
2111	endif
2112	if (w_est(k+1) gt 0.001) then begin
2113	f_star(k+1)=f_star(k)+a_star(k)+entr_star(k)-detr_star(k)
2114	if (k lt lalim) then begin
2115	teta_p(k)=(f_star(k)teta_p(k-1)+(a_star(k)+entr_star(k))0.5*(teta_p(k-1) + teta_les(k,lt_plotindex_les)))/(f_star(k+1)+detr_star(k))
2116	endif else begin
2117	teta_p(k)=(f_star(k)teta_p(k-1)+(a_star(k)+entr_star(k))teta_les(k,lt_plotindex_les))/(f_star(k+1)+detr_star(k))
2118	endelse
2119	zbuoy(k) = grav*(teta_p(k)/teta_les(k,lt_plotindex_les) -1.)
2120	zw2fact=fact_epsilon2.approx_zdz_les(k)/(1.+betalpha)
2121	zdw2=afact*zbuoy(k)/fact_epsilon
2122	zw2(k+1) = MAX([0.0001,exp(-zw2fact)*(zw2(k)-zdw2)+zdw2])
2123	endif
2124	ENDFOR
2125	print, ' -> done'
2126
2127	print, '........ CHECKING VERTICAL VELOCITY FORMULATION'
2128	what_I_plot = make_array(nz,value=0.)
2129	what_I_overplot = make_array(nz,value=0.)
2130	FOR k=0, nz-2 DO BEGIN
2131	what_I_plot(k) = 0.5*(sqrt(zw2(k)) + sqrt(zw2(k+1)))
2132	ENDFOR
2133	FOR k=0, nZmx-2 DO BEGIN
2134	what_I_overplot(k) = 0.5*(zw2_lev(k,lt_plotindex_gcm) + zw2_lev(k+1,lt_plotindex_gcm))
2135	ENDFOR
2136	labels=['zw2 in les calc as in TH']
2137	title_user = TestCase+SubCase+LayerCase+' LES vertical velocity formulation check'
2138	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_w2_check.ps'
2139	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2140	CS, SCALE=28
2141	GSET, XMIN=0., XMAX=8., YMIN=0., YMAX=7, TITLE=title_user
2142	cols=INDGEN(1)+2
2143	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2144	AXES, XSTEP = 1 , XTITLE='w2 in LES from TH calc m/s', YSTEP=1., YTITLE='Altitude (km)',NDECS=3
2145
2146	oplot, w_mean1(*,lt_plotindex_les), altitudes_LES/1000.
2147	oplot, what_I_overplot, altitudes_GCM/1000.,psym =4
2148
2149	PSCLOSE, /NOVIEW
2150
2151	spawn, 'ps2png '+filename
2152
2153	print, '........ CHECKING TETA ESTIMATIONS FORMULATION'
2154
2155	what_I_plot = [[teta_est],[teta_p],[tplume1moy(*,lt_plotindex_les)]]
2156	labels=['LES estimated teta','LES teta plume calc as in TH','LES teta plume']
2157	title_user = TestCase+SubCase+LayerCase+' LES estimated teta formulation check'
2158	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_teta_check.ps'
2159	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2160	CS, SCALE=28
2161	GSET, XMIN=214., XMAX=220., YMIN=0., YMAX=7, TITLE=title_user
2162	cols=INDGEN(3)+2
2163	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2164	AXES, XSTEP = 1 , XTITLE='Teta plume and Est in LES from TH calc K', YSTEP=1., YTITLE='Altitude (km)',NDECS=3
2165
2166	oplot, teta_gcm(,lt_plotindex_gcm)(buoyancy_gcm(*,lt_plotindex_gcm)/grav +1.), altitudes_GCM/1000.
2167	oplot, teta_gcm(,lt_plotindex_gcm)(buoyancy_est_gcm(*,lt_plotindex_gcm)/grav +1.), altitudes_GCM/1000.
2168
2169	PSCLOSE, /NOVIEW
2170
2171	spawn, 'ps2png '+filename
2172
2173	print, '........ CHECKING MASS FLUX FORMULATION'
2174
2175	what_I_plot = [[f_star/MAX(f_star)],[smoothed_fm_trac1_les(,lt_plotindex_les)/MAX(smoothed_fm_trac1_les(,lt_plotindex_les))]]
2176	labels=['LES normalized f_star ','LES normalized updraft mass flux']
2177	title_user = TestCase+SubCase+LayerCase+' LES normalized f_star formulation check'
2178	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fm_check.ps'
2179	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2180	CS, SCALE=28
2181	GSET, XMIN=0., XMAX=1., YMIN=0., YMAX=7, TITLE=title_user
2182	cols=INDGEN(2)+2
2183	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2184	AXES, XSTEP = 0.1 , XTITLE='f/max(f) in LES from TH calc', YSTEP=1., YTITLE='Altitude (km)',NDECS=3
2185
2186	oplot,fm_therm_gcm_interlay(,lt_plotindex_gcm)/MAX(fm_therm_gcm_interlay(,lt_plotindex_gcm)), altitudes_GCM/1000.,psym=4
2187
2188	PSCLOSE, /NOVIEW
2189
2190	spawn, 'ps2png '+filename
2191
2192	; COMPUTING THE CONTINUITY EQUATION IN THE QUASI-BOUSSINESQ APPROX
2193
2194	da_dt = make_array(nz,n_elements(localtime))
2195	smoothed_da_dt = make_array(nz)
2196	FOR k=0, nz-1 DO BEGIN
2197	da_dt(k,) = deriv(localtime,reform(alpha1out(k,)))/3700.
2198	ENDFOR
2199	FOR t=-ns,ns DO BEGIN
2200	smoothed_da_dt = smoothed_da_dt + REFORM(da_dt(*,lt_plotindex_les+t))
2201	ENDFOR
2202	smoothed_da_dt = smoothed_da_dt/nstot
2203
2204	rho = pt/(R*temp_les)
2205	smoothed_rho = make_array(nz)
2206	FOR t=-ns,ns DO BEGIN
2207	smoothed_rho = smoothed_rho + REFORM(rho(*,lt_plotindex_les+t))
2208	ENDFOR
2209	smoothed_rho = smoothed_rho/nstot
2210
2211	continuity1 = smoothed_rhosmoothed_da_dt + df_dz_les1 - smoothed_e_rate_ude_trac1_lessmoothed_fm_trac1_les + smoothed_d_rate_ude_trac1_les*smoothed_fm_trac1_les
2212	print, '........ CONTINUITY CHECK'
2213
2214	what_I_plot = [[continuity1],[smoothed_rhosmoothed_da_dt],[df_dz_les1],[-smoothed_e_rate_ude_trac1_lessmoothed_fm_trac1_les],[smoothed_d_rate_ude_trac1_les*smoothed_fm_trac1_les]]
2215	labels=['total continuity','rho*da/dt','df/dz','-E','D']
2216	title_user = TestCase+SubCase+LayerCase+' LES UDE continuity check, average over '+taverage+' mn'
2217	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_continuity.ps'
2218	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2219	CS, SCALE=28
2220	GSET, XMIN=-0.0005, XMAX=0.0005, YMIN=0, YMAX=10, TITLE=title_user
2221	cols=INDGEN(5)+2
2222	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2223	AXES, XSTEP = 0.0001 , XTITLE='kg.m-2.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2224	PSCLOSE, /NOVIEW
2225
2226	spawn, 'ps2png '+filename
2227
2228	; COMPUTING THE E-D TERM FROM THE CONTINUITY EQUATION
2229
2230	eminusd1=make_array(nz)
2231	FOR k=0, nz-1 DO BEGIN
2232	IF(smoothed_fm_trac1_les(k) ne 0.) THEN eminusd1(k) = (smoothed_rho(k)*smoothed_da_dt(k) - df_dz_les1(k))/smoothed_fm_trac1_les(k) ELSE eminusd1(k)=0.
2233	ENDFOR
2234	what_I_plot = eminusd1
2235	labels=['e-d']
2236	title_user = TestCase+SubCase+LayerCase+' LES e-d, average over '+taverage+' mn'
2237	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_EminusD.ps'
2238	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2239	CS, SCALE=28
2240	GSET, XMIN=-0.002, XMAX=0.002, YMIN=0, YMAX=10, TITLE=title_user
2241	cols=INDGEN(1)+2
2242	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2243	AXES, XSTEP = 0.0005 , XTITLE='kg.m-2.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2244	PSCLOSE, /NOVIEW
2245
2246	spawn, 'ps2png '+filename
2247
2248	; COMPUTING THE TURBULENT FLUX DECOMPOSITION IN PASSIVE ENV AND ACTIVE PLUME
2249	; TO CHECK CONSISTENCY
2250
2251	print, '........ CHECKING CONSISTENCY OF UPDRAFT/ENV DECOMPOSITION'
2252
2253	smoothed_hf1_term1 = make_array(nz)
2254	smoothed_hf1_term2 = make_array(nz)
2255	smoothed_hf1_term3 = make_array(nz)
2256	smoothed_wt = make_array(nz)
2257	FOR t=-ns,ns DO BEGIN
2258	smoothed_hf1_term1 = smoothed_hf1_term1 + REFORM(hf1_term1(*,lt_plotindex_les+t))
2259	smoothed_hf1_term2 = smoothed_hf1_term2 + REFORM(hf1_term2(*,lt_plotindex_les+t))
2260	smoothed_hf1_term3 = smoothed_hf1_term3 + REFORM(hf1_term3(*,lt_plotindex_les+t))
2261	smoothed_wt = smoothed_wt + REFORM(wt(*,lt_plotindex_les+t))
2262	ENDFOR
2263	smoothed_hf1_term1 = smoothed_hf1_term1/nstot
2264	smoothed_hf1_term2 = smoothed_hf1_term2/nstot
2265	smoothed_hf1_term3 = smoothed_hf1_term3/nstot
2266	smoothed_wt = smoothed_wt/nstot
2267
2268	what_I_plot = [[smoothed_hf1_term1],[smoothed_hf1_term2],[smoothed_hf1_term3],[smoothed_hf1_term1+smoothed_hf1_term2+smoothed_hf1_term3]]
2269	labels=['within plume turbulence','within env. turbulence','organized turbulence','TOTAL']
2270	title_user = TestCase+SubCase+LayerCase+' LES turbulence decomposition check, average over '+taverage+' mn'
2271	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_turbu.ps'
2272	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2273	CS, SCALE=28
2274	GSET, XMIN=-1, XMAX=1.5, YMIN=0, YMAX=6, TITLE=title_user
2275	cols=INDGEN(4)+2
2276	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
2277	AXES, XSTEP = 0.5 , XTITLE='m.K/s', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2278	oplot, smoothed_wt, altitudes_LES/1000.,psym=3
2279	PSCLOSE, /NOVIEW
2280	spawn, 'ps2png '+filename
2281
2282
2283	; COMPUTING THE TURBULENT FLUX DECOMPOSITION IN PASSIVE ENV ,ACTIVE PLUME and ACTIVE DOWNDRAFT
2284	; TO CHECK CONSISTENCY
2285
2286	print, '........ CHECKING CONSISTENCY OF UPDRAFT/DOWNDRAFT/ENV DECOMPOSITION'
2287
2288	smoothed_hf1_ude_term1 = make_array(nz)
2289	smoothed_hf1_ude_term2 = make_array(nz)
2290	smoothed_hf1_ude_term3 = make_array(nz)
2291	smoothed_hf1_ude_term4 = make_array(nz)
2292	FOR t=-ns,ns DO BEGIN
2293	smoothed_hf1_ude_term1 = smoothed_hf1_ude_term1 + REFORM(hf1_ude_term1(*,lt_plotindex_les+t))
2294	smoothed_hf1_ude_term2 = smoothed_hf1_ude_term2 + REFORM(hf1_ude_term2(*,lt_plotindex_les+t))
2295	smoothed_hf1_ude_term3 = smoothed_hf1_ude_term3 + REFORM(hf1_ude_term3(*,lt_plotindex_les+t))
2296	smoothed_hf1_ude_term4 = smoothed_hf1_ude_term4 + REFORM(hf1_ude_term4(*,lt_plotindex_les+t))
2297	ENDFOR
2298	smoothed_hf1_ude_term1 = smoothed_hf1_ude_term1/nstot
2299	smoothed_hf1_ude_term2 = smoothed_hf1_ude_term2/nstot
2300	smoothed_hf1_ude_term3 = smoothed_hf1_ude_term3/nstot
2301	smoothed_hf1_ude_term4 = smoothed_hf1_ude_term4/nstot
2302
2303	what_I_plot = [[smoothed_hf1_ude_term1],[smoothed_hf1_ude_term2],[smoothed_hf1_ude_term3],[smoothed_hf1_ude_term4],[smoothed_hf1_ude_term1+smoothed_hf1_ude_term2+smoothed_hf1_ude_term3+smoothed_hf1_ude_term4]]
2304	labels=['within plume turbulence','within downdraft turbulence','within env. turbulence','organized turbulence','TOTAL']
2305	title_user = TestCase+SubCase+LayerCase+' LES UDE turbulence decomposition check, average over '+taverage+' mn'
2306	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_turbu_ude.ps'
2307	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2308	CS, SCALE=28
2309	GSET, XMIN=-1, XMAX=1.5, YMIN=0, YMAX=6, TITLE=title_user
2310	cols=INDGEN(5)+2
2311	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
2312	AXES, XSTEP = 0.25 , XTITLE='m.K/s', YSTEP=0.5, YTITLE='Altitude (km)',NDECS=3
2313	oplot, smoothed_wt, altitudes_LES/1000.,psym=3
2314	PSCLOSE, /NOVIEW
2315
2316	spawn, 'ps2png '+filename
2317
2318	; CHECK CONSISTENCY OF We = W and THETA e = THETA approximation
2319
2320	print, '........ CHECKING CONSISTENCY OF env variable (w_e,theta_e) = mean variable (w_overbar,theta_overbar)'
2321	print, 'as well as mean(w) = 0, in the UDE decomposition'
2322
2323	smoothed_delta_theta_ude = make_array(nz)
2324	smoothed_delta_w_ude = make_array(nz)
2325	smoothed_w_mean1_full = make_array(nz)
2326
2327	FOR t=-ns,ns DO BEGIN
2328	smoothed_delta_theta_ude = smoothed_delta_theta_ude + REFORM(w_mean1_env_ude(,lt_plotindex_les+t)-w_mean1_full(,lt_plotindex_les+t))
2329	smoothed_delta_w_ude = smoothed_delta_w_ude + REFORM(tenv1moy_ude(,lt_plotindex_les+t)-tmoy_full(,lt_plotindex_les+t))
2330	smoothed_w_mean1_full = smoothed_w_mean1_full + REFORM(w_mean1_full(*,lt_plotindex_les+t))
2331	ENDFOR
2332
2333	smoothed_delta_theta_ude = smoothed_delta_theta_ude/nstot
2334	smoothed_delta_w_ude = smoothed_delta_w_ude/nstot
2335	smoothed_w_mean1_full = smoothed_w_mean1_full/nstot
2336
2337	what_I_plot = [[smoothed_delta_theta_ude],[smoothed_delta_w_ude],[smoothed_w_mean1_full]]
2338	labels=['theta env_ude - theta moy','w env_ude - w moy','mean w over domain (,,k)']
2339	title_user = TestCase+SubCase+LayerCase+' LES UDE env/mean approximation check, average over '+taverage+' mn'
2340	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_approx_ude.ps'
2341	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2342	CS, SCALE=28
2343	GSET, XMIN=-2, XMAX=2, YMIN=0, YMAX=10, TITLE=title_user
2344	cols=INDGEN(3)+2
2345	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2346	AXES, XSTEP = 0.5 , XTITLE='(m/s) and (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2347	PSCLOSE, /NOVIEW
2348
2349	spawn, 'ps2png '+filename
2350
2351	; GETTING SOME INSIGHT ON PLUME'S INSIDE TEMPERATURES
2352
2353	print, '........ STRUCTURE POTENTIAL TEMPERATURES'
2354	xmin = 210
2355	xmax = 220
2356	if (TestCase eq 'Case_Z') then begin
2357	xmin = 260
2358	xmax = 270
2359	endif
2360	ztva = teta_gcm*(buoyancy_gcm/grav +1.)
2361	ztva_est = teta_gcm*(buoyancy_est_gcm/grav +1.)
2362	what_I_plot = [[tplume1moy(,lt_plotindex_les)],[tenv1moy(,lt_plotindex_les)],[teta_les(*,lt_plotindex_les)]]
2363	labels=['Teta updraft','Teta env ','Teta moy']
2364	title_user = TestCase+SubCase+LayerCase+' LES Teta in the structures, no average'
2365	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fullTeta.ps'
2366	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2367	CS, SCALE=28
2368	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=6, TITLE=title_user
2369	cols=INDGEN(3)+2
2370	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2371	AXES, XSTEP = 1 , XTITLE='Potential Temperature (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2372	oplot, ztva(*,lt_plotindex_gcm), altitudes_GCM/1000., thick=0.3
2373	oplot, ztva_est(*,lt_plotindex_gcm), altitudes_GCM/1000., thick=0.3
2374	oplot, teta_gcm(*,lt_plotindex_gcm), altitudes_GCM/1000., thick=0.3
2375	PSCLOSE, /NOVIEW
2376	spawn, 'ps2png '+filename
2377
2378
2379	print, '........ UDE STRUCTURE POTENTIAL TEMPERATURES'
2380
2381	what_I_plot = [[tplume1moy(,lt_plotindex_les)],[tdown1moy(,lt_plotindex_les)],[tenv1moy_ude(,lt_plotindex_les)],[teta_les(,lt_plotindex_les)]]
2382	labels=['Teta updraft','Teta downdraft','Teta env (UDE)','Teta moy']
2383	title_user = TestCase+SubCase+LayerCase+' LES UDE Teta in the structures, no average'
2384	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fullTeta_ude.ps'
2385	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2386	CS, SCALE=28
2387	GSET, XMIN=210, XMAX=220, YMIN=0, YMAX=6, TITLE=title_user
2388	cols=INDGEN(4)+2
2389	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2390	AXES, XSTEP = 1 , XTITLE='Potential Temperature (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2391	PSCLOSE, /NOVIEW
2392
2393	spawn, 'ps2png '+filename
2394
2395	; -----------------------------------------------------------------------------------------------------------------------
2396	; End of PLUME diagnostics
2397	; -----------------------------------------------------------------------------------------------------------------------
2398
2399	; * TKE *
2400
2401	print, '........ TKE'
2402
2403	what_I_plot = reform(tke_gcm(*,lt_plotindex_gcm))
2404	labels=['TH tke 1d']
2405	title_user = TestCase+SubCase+LayerCase+' TKE comparison'
2406	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_tke.ps'
2407	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2408	CS, SCALE=28
2409	GSET, XMIN=-1, XMAX=8, YMIN=0, YMAX=10, TITLE=title_user
2410	cols=INDGEN(1)+2
2411	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2412	AXES, XSTEP = 1, XTITLE='Turbulent kinetic energy (kg.m-3)', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
2413
2414	oplot, tke_les(*,lt_plotindex_les), altitudes_LES/1000., psym=4
2415
2416	PSCLOSE, /NOVIEW
2417
2418	spawn, 'ps2png '+filename
2419
2420	; * HEAT FLUX *
2421
2422	print, '........ HEAT FLUX'
2423
2424	lay_heatFlux_up = make_array(nZmx,nTmx)
2425	lay_heatFlux_down = make_array(nZmx,nTmx)
2426
2427	FOR k=1, nZmx-1 DO BEGIN
2428	lay_heatFlux_up(k,) = 0.5(heatFlux_up(k,) + heatFlux_up(k-1,))
2429	lay_heatFlux_down(k,) = 0.5(heatFlux_down(k,) + heatFlux_down(k-1,))
2430	ENDFOR
2431	lay_heatFlux_up(0,)=0.5(heatFlux_up(0,*))
2432	lay_heatFlux_down(0,)=0.5(heatFlux_down(0,*))
2433	zkh_gcm_int = make_array(nZmx)
2434
2435	FOR k=0, nZmx-2 DO BEGIN
2436	zkh_gcm_int(k) = 0.5*(zkh(k,lt_plotindex_gcm) + zkh(k+1,lt_plotindex_gcm))
2437	ENDFOR
2438	MY_gcm = -zkh_gcm_int(deriv(altitudes_GCM,reform(zh(,lt_plotindex_gcm)))); - 0.025max(lay_heatFlux_up(,lt_plotindex_gcm)))
2439
2440
2441	;what_I_plot = [[lay_heatFlux_up(,lt_plotindex_gcm)],[MY_gcm],[lay_heatFlux_up(,lt_plotindex_gcm)+MY_gcm]]
2442	what_I_plot = [[lay_heatFlux_up(,lt_plotindex_gcm)],[MY_gcm],[lay_heatFlux_down(,lt_plotindex_gcm)],[lay_heatFlux_down(,lt_plotindex_gcm)+lay_heatFlux_up(,lt_plotindex_gcm)+MY_gcm]]
2443
2444	;labels=['TH updraft heat flux','Mellor and Yamada gcm heat flux','Total']
2445	labels=['TH updraft heat flux','Mellor and Yamada gcm heat flux','TH downdraft heat flux','Total']
2446
2447	title_user = TestCase+SubCase+LayerCase+' TH vertical heat flux'
2448	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_WT.ps'
2449	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2450	CS, SCALE=28
2451	GSET, XMIN=-2, XMAX=3, YMIN=0, YMAX=10, TITLE=title_user
2452	;cols=INDGEN(3)+2
2453	cols=INDGEN(4)+2
2454	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2455	AXES, XSTEP = 1, XTITLE='vertical turbulent heat flux', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
2456
2457	;hf1_ude_term4 = alpha1out(w_mean1 - w_mean1_full)(tplume1moy - tmoy_full) + beta1out(w_mean1_down - w_mean1_full)(tdown1moy - tmoy_full) + (1.- (alpha1out+beta1out))(w_mean1_env_ude - w_mean1_full)(tenv1moy_ude - tmoy_full)
2458
2459	oplot, alpha1out(,lt_plotindex_les)(w_mean1(,lt_plotindex_les) - w_mean1_full(,lt_plotindex_les))(tplume1moy(,lt_plotindex_les) - tmoy_full(*,lt_plotindex_les)), altitudes_LES/1000., color=2
2460
2461	oplot, beta1out(,lt_plotindex_les)(w_mean1_down(,lt_plotindex_les) - w_mean1_full(,lt_plotindex_les))(tdown1moy(,lt_plotindex_les) - tmoy_full(*,lt_plotindex_les)), altitudes_LES/1000.,color=6
2462
2463	oplot, smoothed_hf1_ude_term4, altitudes_LES/1000., color = 5
2464
2465	;oplot, smoothed_hf1_term1, altitudes_LES/1000.,thick=0.1,LINESTYLE = 5
2466	;oplot, smoothed_hf1_term2, altitudes_LES/1000.,color=2,thick=0.1,LINESTYLE = 5
2467	;oplot, smoothed_hf1_term3, altitudes_LES/1000.,color=8,thick=0.1,LINESTYLE = 5
2468
2469	oplot, wt(*,lt_plotindex_les), altitudes_LES/1000.
2470
2471	;oplot, smoothed_hf1_term3(*,lt_plotindex_les), altitudes_LES/1000., psym=7
2472
2473	PSCLOSE, /NOVIEW
2474
2475	spawn, 'ps2png '+filename
2476
2477	; * TRACERS *
2478
2479	print, '........ TRACER PLOTS DEACTIVATED'
2480
2481	; trying stuff
2482
2483	buoyancy_downdraft = grav*(tdown1moy/tenv1moy_ude-1.)
2484	lmix = make_array(nttot,value=-1.)
2485	altitudes_rel_LES = make_array(nz,nttot)
2486	FOR l=0, nttot-1 DO BEGIN
2487	; kmax = where(w_mean1(,l) eq max(w_mean1(,l)))
2488	; if (kmax(0) ne -1) then lmix(l) = altitudes_LES(kmax(0)) else lmix(l) = -1.
2489	;FOR k=nz-2, 1,-1 DO BEGIN
2490	; if ((buoyancy_downdraft(k,l) gt 0.) and (buoyancy_downdraft(k-1,l) lt 0.)) then lmix(l) = 0.5*(altitudes_LES(k)+altitudes_LES(k+1))
2491	;ENDFOR
2492	FOR k=nz-2, 1,-1 DO BEGIN
2493	if (tdown1moy(k,l) eq 0.) then lmix(l) = altitudes_LES(k)
2494	ENDFOR
2495	ENDFOR
2496
2497	FOR l=0, nttot-1 DO BEGIN
2498	FOR k=0, nz-1 DO BEGIN
2499	altitudes_rel_LES(k,l) = altitudes_LES(k)/lmix(l)
2500	ENDFOR
2501	ENDFOR
2502
2503	print, '........ Teta down / Teta up in UDE'
2504
2505	stuff2=tdown1moy/tenv1moy_ude
2506
2507	what_I_plot = stuff2(*,lt_plotindex_les)
2508	labels=['Teta d/Teta env 12h']
2509	title_user = TestCase+SubCase+LayerCase+' TH trying stuff'
2510	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuff2.ps'
2511	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2512	CS, SCALE=28
2513	GSET, XMIN=0.992, XMAX=1.004, YMIN=0, YMAX=1, TITLE=title_user
2514	cols=INDGEN(1)+2
2515	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2516	AXES, XSTEP = 0.002, XTITLE='Teta d/ Teta env', YSTEP=0.2, YTITLE='Altitude/zi ',NDECS=3
2517
2518	FOR i=0,nttot-1 DO BEGIN
2519	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
2520	ENDFOR
2521	oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.075)/187.931 + 0.9977, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2522	oplot, (altitudes_rel_LES(,lt_plotindex_les))/19.231 + 0.9938, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2523	oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.60)/(-1333) + 1.00025, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2524	PSCLOSE, /NOVIEW
2525	spawn, 'ps2png '+filename
2526
2527	print, '........ Teta down / Teta up in UDE'
2528
2529	stuff2=tdown1moy/tplume1moy
2530
2531	what_I_plot = stuff2(*,lt_plotindex_les)
2532	labels=['Teta d/Teta u 12h']
2533	title_user = TestCase+SubCase+LayerCase+' TH trying stuff'
2534	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuff2.5.ps'
2535	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2536	CS, SCALE=28
2537	GSET, XMIN=0.95, XMAX=1.1, YMIN=0, YMAX=1, TITLE=title_user
2538	cols=INDGEN(1)+2
2539	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2540	AXES, XSTEP = 0.01, XTITLE='Teta d/ Teta u', YSTEP=0.2, YTITLE='Altitude/zi ',NDECS=3
2541
2542	FOR i=0,nttot-1 DO BEGIN
2543	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
2544	ENDFOR
2545	;oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.075)/187.931 + 0.9977, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2546	;oplot, (altitudes_rel_LES(,lt_plotindex_les))/19.231 + 0.9938, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2547	;oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.60)/(-1333) + 1.00025, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2548	PSCLOSE, /NOVIEW
2549	spawn, 'ps2png '+filename
2550
2551	print, '........ B down / B up in UDE'
2552
2553	stuff2 = (tdown1moy/tenv1moy_ude -1.)/(tplume1moy/tenv1moy_ude -1.)
2554	what_I_plot = stuff2(*,lt_plotindex_les)
2555	labels=['B down/B up 12h']
2556	title_user = TestCase+SubCase+LayerCase+' TH trying stuff Bd/Bu'
2557	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuffBuBd.ps'
2558	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2559	CS, SCALE=28
2560	GSET, XMIN=-1, XMAX=1., YMIN=0, YMAX=1, TITLE=title_user
2561	cols=INDGEN(1)+2
2562	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2563	AXES, XSTEP = 0.1, XTITLE='B down/ B up', YSTEP=0.1, YTITLE='Altitude/zi ',NDECS=1
2564
2565	FOR i=0,nttot-1 DO BEGIN
2566	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
2567	ENDFOR
2568	;oplot, ((altitudes_rel_LES(,lt_plotindex_les)-0.06)/0.839841)^2 - 0.3, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2569	oplot, ((altitudes_rel_LES(,lt_plotindex_les)-0.06)/1.16847)^2 - 0.3, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2570	oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.7)/1., altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2571	;oplot, (altitudes_rel_LES(,lt_plotindex_les))/0.08333333-1., altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2572	oplot, sqrt(altitudes_rel_LES(,lt_plotindex_les)/0.122449)-1., altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2573	PSCLOSE, /NOVIEW
2574	spawn, 'ps2png '+filename
2575
2576	print, '........ F down / F up in UDE'
2577
2578	stuff2 = downward_flux1/fm_trac1_les
2579	what_I_plot = stuff2(*,lt_plotindex_les)
2580	labels=['f down/f up 12h']
2581	title_user = TestCase+SubCase+LayerCase+' TH trying stuff f down/f up'
2582	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stufffufd.ps'
2583	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2584	CS, SCALE=28
2585	GSET, XMIN=-8, XMAX=0.5, YMIN=0, YMAX=1, TITLE=title_user
2586	cols=INDGEN(1)+2
2587	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2588	AXES, XSTEP = 0.5, XTITLE='f down/ f up', YSTEP=0.1, YTITLE='Altitude/zi ',NDECS=1
2589
2590	FOR i=0,nttot-1 DO BEGIN
2591	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
2592	ENDFOR
2593	oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.0149259)/0.00333)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2594	;oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.02)/0.006)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2595	PSCLOSE, /NOVIEW
2596	spawn, 'ps2png '+filename
2597
2598	print, '........ dFdz down / dzFdz up in UDE'
2599	stuff3=make_array(nz,nttot)
2600	FOR l=0, nttot-1 DO BEGIN
2601	stuff3(,l) = deriv(altitudes_LES,downward_flux1(,l))/deriv(altitudes_LES,fm_trac1_les(*,l))
2602	ENDFOR
2603	what_I_plot = stuff3(*,lt_plotindex_les)
2604	labels=['dfdz down/dfdz up 12h']
2605	title_user = TestCase+SubCase+LayerCase+' TH trying stuff f down/f up'
2606	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuffdfudfd.ps'
2607	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2608	CS, SCALE=28
2609	GSET, XMIN=-30, XMAX=30, YMIN=0, YMAX=1, TITLE=title_user
2610	cols=INDGEN(1)+2
2611	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2612	AXES, XSTEP = 5, XTITLE='dfdz down/ dfdz up', YSTEP=0.1, YTITLE='Altitude/zi ',NDECS=1
2613
2614	FOR i=0,nttot-1 DO BEGIN
2615	if(lmix(i) ne -1) then oplot, stuff3(,i), altitudes_rel_LES(,i), thick=0.1
2616	ENDFOR
2617	;oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.0149259)/0.00333)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2618	;oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.02)/0.006)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
2619	PSCLOSE, /NOVIEW
2620	spawn, 'ps2png '+filename
2621	;
2622	;what_I_plot = [[ar_col],[co2_col],[tke_col]]
2623	;labels=['Ar deviation','Co2 deviation','TKE deviation']
2624	;title_user = TestCase+SubCase+' TH 1d tracer conservation'
2625	;filename = TestCase+SubCase+'Gcm_Les_Comp_tracer.ps'
2626	;PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2627	;CS, SCALE=28
2628	;GSET, XMIN=localtime_gcm(0), XMAX=localtime_gcm(nTmx-1), YMIN=0, YMAX=2, TITLE=title_user
2629	;cols=INDGEN(3)+2
2630	;GPLOT, X=localtime_gcm, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2631	;AXES, XSTEP = 2, XTITLE='Local time (h)', YSTEP=0.1, YTITLE='Tracer integrated column mass deviation from origin (%)',NDECS=1
2632	;
2633	;PSCLOSE, /NOVIEW
2634	;
2635	;spawn, 'ps2png '+filename
2636	;
2637	;title_user = TestCase+SubCase+' TH argon propagation from first layer'
2638	;PS_START, file = TestCase+SubCase+'Gcm_Les_Comp_Argon.ps'
2639	;
2640	;what_I_plot = transpose(ar(0:6,*))
2641	;
2642	;maxfield_init = 0.05
2643	;minfield_init = 0
2644	;pal=33
2645	;lim_max = maxfield_init & w=where((what_I_plot ge lim_max) and (what_I_plot le 1e9)) & if (w[0] ne -1) then what_I_plot[w]=lim_max
2646	;lim_min = minfield_init & w=where(what_I_plot le lim_min) & if (w[0] ne -1) then what_I_plot[w]=lim_min
2647	;
2648	;section, $
2649	; what_I_plot, $ ; 2D field
2650	; localtime_gcm, $ ; horizontal coordinate
2651	; altitudes_gcm(0:6), $ ; altitude coordinate
2652	; minfield=minfield_init, $ ; minimum value of plotted field (=0: calculate)
2653	; maxfield=maxfield_init, $ ; maximum value of plotted field (=0: calculate)
2654	;; minspace=minspace, $ ; minimum value of space window (=0: calculate)
2655	;; maxspace=maxspace, $ ; maximum value of space window (=0: calculate)
2656	;; overcontour=overcontour, $ ; another 2D field to overplot with contour lines (=0: no)
2657	;; overvector_x=overvector_x, $ ; wind vector - x component (=0: no)
2658	;; overvector_y=overvector_y, $ ; wind vector - y component (=0: no)
2659	;; colors=colors, $ ; number of colors/levels (32 is default)
2660	; title_plot=title_user, $ ; title of the plot ('Profile' is default)
2661	; title_axis=['Martian hour (h)','Height above ground (m)'], $ ; title of the [x,y] axis (['Field','Altitude'] is default)
2662	; ct=pal, $ ; color table (33-rainbow is default)
2663	;; topo=topography, $
2664	; format=format ; format of colorbar annotations ('(F6.2)' is default)
2665	;
2666	;PS_END, /PNG
2667	;
2668	;INTERVAL_VOLUME, supermask1, 0.5, 1.,verts, conn
2669	;conn = TETRA_SURFACE(verts, conn)
2670	;oRain = OBJ_NEW('IDLgrPolygon', verts, POLYGONS=conn, $
2671	; COLOR=[255,255,255], SHADING=1)
2672	;XOBJVIEW, oRain, BACKGROUND=[150,200,255]
2673
2674	;INTERVAL_VOLUME, supermask2, 0.5, 1.5,verts, conn
2675	;conn = TETRA_SURFACE(verts, conn)
2676	;oRain = OBJ_NEW('IDLgrPolygon', verts, POLYGONS=conn, $
2677	; COLOR=[255,255,255], SHADING=1)
2678	;XOBJVIEW, oRain, BACKGROUND=[150,200,255]
2679
2680	ENDELSE
2681
2682	print, ''
2683	print, '........ ALL DONE'
2684	print, ''
2685
2686	END

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

Download in other formats:

Original Format