Context Navigation

gettherm.pro @ 638

Last change on this file since 638 was 211, checked in by aslmd, 14 years ago
MESOSCALE_DEV: thermal graphical stuff was lost in the clean-up. this is fixed. plus checked if no file is missing anymore, OK.
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