Context Navigation

thermiques.pro @ 824

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

Line
1	FUNCTION myfunct, Bw2, A
2	G = A[0]*Bw2 + A[1]
3	RETURN,[ [G], [Bw2] , [1.0] ]
4	END
5
6	FUNCTION myfunct2, Gammaw2, A
7	G = A[0]*(Gammaw2)^A[1]
8	RETURN,[ [G], [(Gammaw2)^A[1]] , [G*alog(Gammaw2)] ]
9	END
10
11	FUNCTION myfunct3, Gammaw2, A
12	G = A[0]*(Gammaw2)^A[1]+A[2]
13	RETURN,[ [G], [(Gammaw2)^A[1]] , [G*alog(Gammaw2)] , [1.0] ]
14	END
15
16	FUNCTION myfunct4, Bw2, A
17	G = A*Bw2
18	RETURN,[[G], [Bw2]]
19	END
20
21	PRO thermiques
22
23	spawn, 'clear'
24	print, ''
25	print, ' Thermals Analysis '
26	print, ' (usine à gaz) '
27	print, ''
28
29	full='true'
30	f_offset='false'
31	overplot_convadj='false'
32	plot_3d = 'false'
33	les_special=''
34	lctu_gcm = 6.
35	got_swdownz='false'
36	got_tracer_flux='false'
37	got_hfx='false'
38
39	;got_swdownz='true'
40	;got_hfx='true'
41	;got_tracer_flux='true'
42
43	; *********** Best values for LES thermals with tau =0.5
44	betalpha = 1.3
45	afact = 1.8
46	fact_epsilon = 0.0008
47	detr_min = 0.0007
48
49
50	;betalpha = 1.
51	;afact = 2.4
52	;fact_epsilon = 0.0007
53	;detr_min = 0.0007
54
55	; --------
56
57	; *********** Best values for LES thermals with tau =1.
58	;betalpha = 1.
59	;afact = 1.1
60	;fact_epsilon = 0.00025
61	;detr_min = 0.0007
62	; --------
63
64	;datname='thermiques.dat.scale1.2'
65	;datname='thermiques.dat.scale1.4'
66	;datname='thermiques.dat.scale0.6'
67	datname='thermiques.dat'
68	;datname='thermiques.dat.scale0.8'
69	;datname='thermiques.dat' ;scale =1.0, sigmao =1.0
70
71	ns = 0 ; number of points for time-smoothing of LES data : 2*ns+1 points, ns = 9 eq to 30mn (-15mn//+15mn)
72	nstot = float(2.*ns+1.)
73
74	GcmSubCase = ''
75	LayerCase=''
76	s_trac1 = 'qtrac1'
77	s_trac2 = 'qtrac2'
78	got_pdt = 'true'
79	TestCase = 'Case_A'
80	SubCase = '_11_shorter'
81	Histo = 'true'
82	newtest = ''
83	visualization_mode = 'false'
84	got_updrafts='false'
85
86	label_init:
87	spawn, 'clear'
88	print, ' Available simulations :'
89	print, '----'
90	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'
91	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'
92	print, ' 3/ Case_A_4_shorter : Case A4 with Dtrac1 = 5 mn and Dtrac2 = 10 mn (compared to 20 and 100)'
93	print, ' 4/ Case_A_4_shorter_winds : Case A4_shorter with bckgrnd wind u=10 m/s'
94	print, ' 5/ Case_A_4_shorter_winds_30 : Case A4_shorter with bckgrnd wind u=30 m/s'
95	print, ' 6/ Case_A_4_shorter_winds_tau1 : Case A4_shorter with bckgrnd wind u=10 m/s and tau=1'
96	print, ' 7/ Case_A_4_shorter_winds_tau2 : Case A4_shorter with bckgrnd wind u=10 m/s and tau=2'
97	print, ' 8/ Case_ExtremeCase : 101x101x201 ztop=15km,dx=100m,dz=75m,Ls=0°,(0.N;0.E),50tiu,A=0.1,Tau=0.05'
98	print, '----'
99	print, ' 9/ Case_C_4_shorter_winds : '
100	print, ' 10/ Case_I_4_shorter_winds : '
101	print, ' 11/ Case_Z_4_shorter_winds : '
102	print, ' 12/ 1D : 124 layers'
103	print, ' 13/ 1D : 32 layers'
104	print, ' 14/ 1D : low dt'
105	print, ' 15/ LES HR run (257x257x301)'
106	print, ' 16/ 1D : 13 layers : gcm to 15km'
107	print, ''
108	print, ' 0/ PLUME VISUALISATION : '+visualization_mode
109	print, ' 999/ CLEAR thermiques.dat for considered case'
110	print, ' 100/ OVERPLOT CONVADJ ONLY RESULTS : '+overplot_convadj
111	print, ''
112	print, ' SIMULATION NUMBER : '
113	print, ' '+TestCase+SubCase+LayerCase+' '
114	print, ''
115	print, 'Any change ? (number of new case to change, or any other key to continue)'
116	read, newtest
117	if (newtest eq '1') then begin
118	TestCase = 'Case_A'
119	SubCase = '_11'
120	pGround = 867.5594
121	goto,label_init
122	endif
123	if (newtest eq '2') then begin
124	TestCase = 'Case_A'
125	SubCase = '_4'
126	got_pdt = 'false'
127	pGround = 867.5594
128	goto,label_init
129	endif
130	if (newtest eq '3') then begin
131	TestCase = 'Case_A'
132	SubCase = '_4_shorter'
133	s_trac1 = 'qtrac2'
134	s_trac2 = 'qtrac1'
135	full = 'true'
136	f_offset='false'
137	pGround = 867.5594
138	goto,label_init
139	endif
140	if (newtest eq '4') then begin
141	TestCase = 'Case_A'
142	SubCase = '_4_shorter_wind'
143	s_trac1 = 'qtrac2'
144	s_trac2 = 'qtrac1'
145	full = 'true'
146	f_offset='false'
147	lctu_gcm = 6.
148	pGround = 856.997 ;6h
149	;pGround = 867.5594 ;8h
150	;got_swdownz='true'
151	;got_hfx='true'
152	;got_tracer_flux='true'
153	goto,label_init
154	endif
155	if (newtest eq '5') then begin
156	TestCase = 'Case_A'
157	SubCase = '_4_shorter_wind_30'
158	GcmSubCase = '_wind_30'
159	s_trac1 = 'qtrac2'
160	s_trac2 = 'qtrac1'
161	full = 'true'
162	pGround = 856.997 ;6h
163	f_offset='false'
164	goto,label_init
165	endif
166	if (newtest eq '6') then begin
167	TestCase = 'Case_A'
168	SubCase = '_4_shorter_wind_tau1'
169	GcmSubCase = '_tau1'
170	s_trac1 = 'qtrac2'
171	s_trac2 = 'qtrac1'
172	full = 'true'
173	pGround = 856.997 ;6h
174	f_offset='false'
175	goto,label_init
176	endif
177	if (newtest eq '7') then begin
178	TestCase = 'Case_A'
179	SubCase = '_4_shorter_wind_tau2'
180	GcmSubCase = '_tau2'
181	s_trac1 = 'qtrac2'
182	s_trac2 = 'qtrac1'
183	full = 'true'
184	f_offset='false'
185	pGround = 856.997 ;6h
186	goto,label_init
187	endif
188	if (newtest eq '8') then begin
189	TestCase = 'ExtremeCase'
190	SubCase = ''
191	GcmSubCase = ''
192	s_trac1 = 'qtrac2'
193	s_trac2 = 'qtrac1'
194	full = 'true'
195	f_offset='false'
196	pGround = 670.018 ;6h
197	;pGround = 677.722 ;8h
198	lctu_gcm = 6.
199	;got_swdownz='true'
200	;got_hfx='true'
201	;got_tracer_flux='true'
202	goto,label_init
203	endif
204	if (newtest eq '9') then begin
205	TestCase = 'Case_C'
206	SubCase = '_4_shorter_wind'
207	s_trac1 = 'qtrac2'
208	s_trac2 = 'qtrac1'
209	full = 'true'
210	f_offset='false'
211	pGround = 480.6
212	goto,label_init
213	endif
214	if (newtest eq '10') then begin
215	TestCase = 'Case_I'
216	SubCase = '_4_shorter_wind'
217	s_trac1 = 'qtrac2'
218	s_trac2 = 'qtrac1'
219	full = 'true'
220	f_offset='false'
221	pGround = 630.
222	goto,label_init
223	endif
224	if (newtest eq '11') then begin
225	TestCase = 'Case_Z'
226	SubCase = '_4_shorter_wind'
227	s_trac1 = 'qtrac2'
228	s_trac2 = 'qtrac1'
229	full = 'true'
230	f_offset='false'
231	pGround = 266.
232	goto,label_init
233	endif
234	if (newtest eq '12') then begin
235	LayerCase=''
236	goto,label_init
237	endif
238	if (newtest eq '13') then begin
239	LayerCase='_32lev'
240	goto,label_init
241	endif
242	if (newtest eq '14') then begin
243	LayerCase='_low_dt'
244	goto,label_init
245	endif
246	if (newtest eq '15') then begin
247	les_special='_HR'
248	got_swdownz='true'
249	got_hfx='true'
250	got_tracer_flux='true'
251	goto,label_init
252	endif
253	if (newtest eq '16') then begin
254	LayerCase='_13lev'
255	got_swdownz='true'
256	goto,label_init
257	endif
258	if (newtest eq '17') then begin
259	TestCase = 'Exomars'
260	SubCase = ''
261	GcmSubCase = ''
262	s_trac1 = 'qtrac2'
263	s_trac2 = 'qtrac1'
264	full = 'true'
265	f_offset='false'
266	;pGround = 717.899 ;6h
267	pGround = 718.
268	lctu_gcm = 6.
269	;got_swdownz='true'
270	;got_hfx='true'
271	;got_tracer_flux='true'
272	goto,label_init
273	endif
274
275	if (newtest eq '100') then begin
276	if (overplot_convadj eq 'true') then overplot_convadj = 'false' else overplot_convadj = 'true'
277	goto,label_init
278	endif
279	if (newtest eq '0') then begin
280	visualization_mode = 'true'
281	spawn, 'clear'
282	print, 'The first timestep of the considered file will be shown.'
283	print, 'Defaut file is : file 6 (lt ~ 13h10)'
284	print, 'Which file do you want ? (lt ~= file_number + 7)'
285	loop_special = '6'
286	read, loop_special
287	print, 'Do you wish to plot histograms or manipulate volumic data ?'
288	print, '1/ Histogram'
289	print, '2/ Volumic data'
290	x=''
291	read, x
292	if (x eq '1') then Histo='true' else Histo='false'
293	goto,label_init
294	endif
295
296	les_path='/san0/acolmd/SIMUS/LES_'+TestCase+SubCase+les_special
297	;gcm_path='/san0/acolmd/SIMUS/GCM_'+TestCase+LayerCase+'_2'
298	gcm_path='/san0/acolmd/SIMUS/GCM_'+TestCase+GcmSubCase+LayerCase
299	gcm_convadj_path=gcm_path+'_convadj'
300
301	if (newtest eq '999') then spawn, 'rm -f '+les_path+'/'+datname
302	print, ''
303	print, ' -- Loading LES data -- '
304
305	print, 'LES DATA IN : '
306	print, les_path
307	print, 'GCM DATA IN : '
308	print, gcm_path
309
310	p0=610. & t0=220. & r_cp=1./3.89419 & grav=3.72 & R=191.182
311	history_interval_s = 100.
312	;lctu_gcm = 8. ; Initial local time of gcm 1d simu
313	scale = 1. ; Scaling factor for conditional sampling
314	decimate = 10. ; Coeff for subsampling the data for sigma integral
315	sigmao= 1. ; multiplicative coeff for the computation of Sigma0 in the CS
316	sigmao_ude = 0.3 ; number of standard deviation away from mean for the selection of downdraft in UDE
317
318	openr,unit,les_path+'/'+datname,/get_lun,error=err
319	IF (err ne 0) THEN BEGIN
320
321	OPENR, 22, les_path+'/input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22
322	OPENR, 23, les_path+'/input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
323
324	domain='d01' & filesWRF = FindFile(les_path+'/wrfout_'+domain+'_????-??-??_??:??:??') & nf=n_elements(filesWRF)
325	;print, filesWRF
326	;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)
327	;ce fichier utilise aussi offset_localtime = 3.05
328
329	; 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
330	; we switch the names of trac1 and trac2 in the initialization of this routine, in the "case".
331
332	id=ncdf_open(filesWRF(0))
333	NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east' ), dummy, nx & NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north' ), dummy, ny
334	NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top' ), dummy, nz & NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), dummy, nt
335	NCDF_CLOSE, id
336	id=ncdf_open(filesWRF(nf-1)) ;; for interrupted runs
337	NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), dummy, ntlast
338	NCDF_CLOSE, id
339	nttot = (nf-1)*nt + ntlast
340	wt = fltarr(nz,nttot)
341	tke_les = fltarr(nz,nttot) & ztke = fltarr(nz,nttot) & t = fltarr(nz,nttot)
342	p = fltarr(nz) & ph = fltarr(nz) & pht = fltarr(nz,nttot) & pt = fltarr(nz,nttot)
343	xtke = fltarr(nz,nttot) & ytke = fltarr(nz,nttot) & temp_les = fltarr(nz,nttot)
344	wmax = fltarr(nttot)
345	alpha1 = fltarr(nz) & alpha2 = fltarr(nz)
346	alpha1out = fltarr(nz,nttot) & alpha2out = fltarr(nz,nttot)
347	zqtrac1 = dblarr(nx,ny,nz) & zqtrac2 = dblarr(nx,ny,nz)
348	sigmazqtrac1 = fltarr(nz) & sigmazqtrac2 = fltarr(nz)
349	sigmazminqtrac1 = fltarr(nz) & sigmazminqtrac2 = fltarr(nz)
350	fm_trac1_les = fltarr(nz,nttot) & fm_trac2_les = fltarr(nz,nttot)
351	anomalqtrac1 = fltarr(nx,ny,nz) & anomalqtrac2 = fltarr(nx,ny,nz)
352	e_trac1_les = fltarr(nz,nttot) & e_trac2_les = fltarr(nz,nttot)
353	dtempdztmp = fltarr(nx,ny,nz)
354	localtime = lctu + history_interval_s*findgen(nttot)/3700.
355	w_mean1 = fltarr(nz,nttot) & w_mean2 = fltarr(nz,nttot)
356	w_mean1_env = fltarr(nz,nttot) & w_mean1_down = fltarr(nz,nttot)
357	w_mean1_env_ude = fltarr(nz,nttot) & w_mean1_full = fltarr(nz,nttot)
358	buoyancy1_les = fltarr(nz,nttot) & buoyancy2_les = fltarr(nz,nttot)
359	e1_term2 = fltarr(nz,nttot) & e1_term3 = fltarr(nz,nttot)
360	dtetadttmp = fltarr(nx,ny,nz)
361	rhomoy1 = fltarr(nz,nttot)
362	plumeIndex1out = make_array(nxny,nz,VALUE=-1.) & envIndex1out = make_array(nxny,nz,VALUE=-1.)
363	hf1tmp = fltarr(nz,nttot) & hf1tmpenv = fltarr(nz,nttot)
364	tplume1moy = fltarr(nz,nttot) & tenv1moy = fltarr(nz,nttot) & tenv1moy_ude = fltarr(nz,nttot)
365	tmoy_full = fltarr(nz,nttot) & tdown1moy = fltarr(nz,nttot)
366	dteta1moydt_entr = fltarr(nz,nttot) & dteta1moydt_detr = fltarr(nz,nttot)
367	d1_term1 = fltarr(nz,nttot) & d1_term2 = fltarr(nz,nttot) & d1_term3=fltarr(nz,nttot)
368	hf1_term1 = fltarr(nz,nttot) & hf1_term2 = fltarr(nz,nttot) & hf1_term3 = fltarr(nz,nttot)
369	d1_term1_ude = fltarr(nz,nttot) & d1_term2_ude = fltarr(nz,nttot) & d1_term3_ude=fltarr(nz,nttot)
370	e1_term1_ude = fltarr(nz,nttot) & e1_term2_ude = fltarr(nz,nttot) & e1_term3_ude=fltarr(nz,nttot)
371	downward_flux1 = fltarr(nz,nttot) & beta1out = fltarr(nz,nttot)
372	hf1_ude_term1 = fltarr(nz,nttot) & hf1_ude_term2 = fltarr(nz,nttot) & hf1_ude_term3 = fltarr(nz,nttot) & hf1_ude_term4 = fltarr(nz,nttot)
373	hf1tmpenv_ude = fltarr(nz,nttot) & hf1tmp_down = fltarr(nz,nttot)
374	dTeta_phys = make_array(nz,nttot)
375	exner = fltarr(nz,nttot)
376	uv_moy = fltarr(nz,nttot)
377	Gamma_1 = fltarr(nz,nttot) & Gamma_2 = fltarr(nz,nttot) & Gamma_3 = fltarr(nz,nttot)
378	Gamma_1_tmp = fltarr(nz,nttot)
379	ptotprime = fltarr(nx,ny,nz) & anomalptot = fltarr(nx,ny,nz) & dptotprimedztmp = fltarr(nx,ny,nz)
380	l=0
381	tsurf_les = fltarr(nttot)
382
383	FOR loop = 0, nf-1 DO BEGIN
384	timetime = SYSTIME(1)
385	if (loop ne nf-1) then nloop2=nt else nloop2=ntlast
386	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
387	FOR loop2 = loop2_init, nloop2-1 DO BEGIN
388	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
389	ph = TEMPORARY(ph) + pht(*,l) / (nttot-1)
390	ENDFOR
391	print, 'computing altitudes, file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s'
392	ENDFOR
393	altitudes_LES = 1000.*(TEMPORARY(ph) - hgtu/1000.) ;; altitude above ground
394	pht = fltarr(nz,nttot)
395	ph = fltarr(nz)
396
397	FOR loop = 0, nf-1 DO BEGIN
398	timetime = SYSTIME(1)
399	if (loop ne nf-1) then nloop2=nt else nloop2=ntlast
400	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
401	FOR loop2 = loop2_init, nloop2-1 DO BEGIN
402
403	anomalt = 1. & anomalu = 1. & anomalv = 1. & anomalw = 1.
404	; --------------------------------------------------------
405	; u' = u and v' = v (car PAS de background wind !)
406	; tke = 0.5 ( <u'^2> + <v'^2> + <w'^2> ) ; u' = u ; v' = v
407	; --------------------------------------------------------
408
409	tprime = getget(filesWRF(loop), 'T', anomaly=anomalt, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; t' = t - <t>
410	t(*,l) = t0 + temporary(anomalt)
411	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)
412	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)
413	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)
414	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)
415	tke_les(,l) = xtke(,l) + ytke(,l) + ztke(,l)
416	wprime = getget(filesWRF(loop), 'W', anomaly=anomalw, count=[0,0,0,1], offset=[0,0,0,loop2])
417	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
418	pt(*,l) = TOTAL(TOTAL(getget(filesWRF(loop), 'PTOT' , count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny)
419	temp_les(,l) = t(,l)(pt(,l)/p0)^r_cp
420	IF (got_pdt eq 'true') then begin
421	exner(,l) = (pt(,l)/p0)^r_cp
422	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)
423	ENDIF
424	ph = TEMPORARY(ph) + pht(*,l) / (nttot-1)
425	p = TEMPORARY(p ) + pt(*,l) / nttot
426	tsurf_les(l)=TOTAL(TOTAL(getget(filesWRF(loop), 'TSURF', count=[0,0,1], offset=[0,0,loop2]),1),1)/ float(nx) / float(ny)
427
428	IF (got_updrafts EQ 'true') THEN BEGIN
429
430	ptotprime(,,*) = getget(filesWRF(loop), 'PTOT', count=[0,0,0,1], offset=[0,0,0,loop2])
431	FOR k=0, nz-1 DO BEGIN
432	rhomoy1(,l) = TOTAL(TOTAL(reform(((ptotprime(,,k)/(R(t(k,l)+tprime(,,k))))(p0/ptotprime(,,k))^r_cp)),1),1)/(float(nx)float(ny))
433	anomalptot(,,k) = ptotprime(,,k) - pt(k,l)
434	ENDFOR
435	zqtrac1 = getget(filesWRF(loop), s_trac1, count=[0,0,0,1], offset=[0,0,0,loop2])
436	FOR i=0,nx-1 DO BEGIN
437	FOR j=0, ny-1 DO BEGIN
438	dtempdztmp(i,j,) = deriv(altitudes_LES, tprime(i,j,) + t(*,l))
439	dptotprimedztmp(i,j,) = deriv(altitudes_LES, anomalptot(i,j,))
440	ENDFOR
441	ENDFOR
442	FOR k=0, nz-1 DO BEGIN
443	Gamma_2(k,l) = -(TOTAL(TOTAL(dptotprimedztmp(,,k),1),1)/float(nx) /float(ny))/rhomoy1(k,l)
444	Gamma_3(k,l) = -grav(TOTAL(TOTAL(anomalptot(,*,k),1),1)/float(nx) /float(ny))/pt(k,l)
445	ENDFOR
446	FOR k=0,nz-1 DO BEGIN
447	anomalqtrac1(,,k) = zqtrac1(,,k) - TOTAL(TOTAL(REFORM(zqtrac1(,,k)),1),1)/ float(nx) / float(ny)
448	sigmazqtrac1(k) = STDDEV(REFORM(zqtrac1(,,k)))
449	IF (k ne 0) THEN BEGIN
450	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
451	sigmazminqtrac1(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac1(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
452	ENDIF ELSE BEGIN
453	sigmazminqtrac1(k) = sigmazqtrac1(k)
454	ENDELSE
455	; plumeIndex1 = WHERE((anomalqtrac1(,,k) GT scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) AND ((anomalw(k)+wprime(,*,k)) GT 0.))
456	; envIndex1 = WHERE((anomalqtrac1(,,k) LE scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) OR ((anomalw(k)+wprime(,*,k)) LE 0.))
457	plumeIndex1 = WHERE(anomalqtrac1(,,k) GT scale*MAX([sigmazqtrac1(k),sigmazminqtrac1(k)]))
458	envIndex1 = WHERE(anomalqtrac1(,,k) LE scale*MAX([sigmazqtrac1(k),sigmazminqtrac1(k)]))
459	IF(plumeIndex1(0) EQ -1) THEN BEGIN
460	fm_trac1_les(k,l)=0.
461	e_trac1_les(k,l)=0.
462	alpha1out(k,l)=0.
463	buoyancy1_les(k,l)=0.
464	w_mean1(k,l)=0.
465	w_mean1_env(k,l)=0.
466	w_mean1_down(k,l)=0.
467	w_mean1_full(k,l)=0.
468	w_mean1_env_ude(k,l)=0.
469	e1_term2(k,l)=0.
470	e1_term3(k,l)=0.
471	e1_term1_ude(k,l)=0.
472	e1_term2_ude(k,l)=0.
473	e1_term3_ude(k,l)=0.
474	hf1tmp(k,l)=0.
475	hf1tmpenv(k,l)=0.
476	plumeIndex1out(*,k)=-1.
477	envIndex1out(*,k)=-1.
478	d1_term1(k,l)=0.
479	d1_term2(k,l)=0.
480	d1_term3(k,l)=0.
481	d1_term1_ude(k,l)=0.
482	d1_term2_ude(k,l)=0.
483	d1_term3_ude(k,l)=0.
484	downward_flux1(k,l)=0.
485	beta1out(k,l)=0.
486	tmoy_full(k,l)=0.
487	tdown1moy(k,l)=0.
488	ENDIF ELSE BEGIN
489	FOR n=0,n_elements(plumeIndex1)-1 DO BEGIN
490	plumeIndex1out(n,k)=plumeIndex1(n)
491	ENDFOR
492	FOR n=0,n_elements(envIndex1)-1 DO BEGIN
493	envIndex1out(n,k)=envIndex1(n)
494	ENDFOR
495	alpha1(k) = n_elements(plumeIndex1) / float(nx) / float(ny)
496	wprimetmp = reform(reform((anomalw(k)+wprime(,,k))),[nx*ny,1])
497	w_mean1_full(k,l) = mean(wprimetmp)
498	w_mean1(k,l) = mean(wprimetmp(plumeIndex1))
499	w_mean1_env(k,l) = mean(wprimetmp(envIndex1))
500	downdraft_index1 = WHERE((abs(anomalw(k)+wprime(,,k)) gt sigmao_udeSTDDEV(wprimetmp(envIndex1))) and (anomalw(k)+wprime(,*,k) lt 0.))
501
502	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.)))
503	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.
504	if (downdraft_index1(0) ne -1) then begin
505	w_mean1_down(k,l)=mean(wprimetmp(downdraft_index1))
506	wprimetmp=0.
507	beta1 = n_elements(downdraft_index1) / float(nx) / float(ny)
508	beta1out(k,l)=beta1
509	downward_flux1(k,l) = beta1rhomoy1(k,l)w_mean1_down(k,l)
510	endif else begin
511	downward_flux1(k,l)=0.
512	beta1out(k,l)=0.
513	w_mean1_down(k,l)=0.
514	tdown1moy(k,l)=0.
515	endelse
516	fm_trac1_les(k,l) = alpha1(k)rhomoy1(k,l)w_mean1(k,l)
517	dtempdztmplin = reform(reform(dtempdztmp(,,k)),[nx*ny,1])
518	alpha1out(k,l)=alpha1(k)
519	tfull=reform(tprime(,,k)+t(k,l),[nx*ny,1])
520	if (downdraft_index1(0) ne -1) then tdown1moy(k,l)=mean(tfull(downdraft_index1))
521	tplume1moy(k,l)=mean(tfull(plumeIndex1))
522	tenv1moy(k,l)=mean(tfull(envIndex1))
523	if (envIndex1_ude(0) ne -1) then tenv1moy_ude(k,l) = mean(tfull(envIndex1_ude)) else tenv1moy_ude(k,l)=0.
524	tmoy_full(k,l) = mean(tfull)
525	buoyancy1_les(k,l)=grav*(tplume1moy(k,l)/tenv1moy(k,l)-1.)
526	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))
527	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))
528	if (envIndex1_ude(0) ne -1) then begin
529	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))
530	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))
531	endif else begin
532	e1_term1_ude(k,l) = 0.
533	d1_term1_ude(k,l) = 0.
534	endelse
535	wtmp=reform(wprime(,,k)+anomalw(k),[nx*ny,1])
536	ttmp=reform(tprime(,,k)+t(k,l),[nx*ny,1])
537	Gamma_1_tmp (k,l) = alpha1out(k,l)rhomoy1(k,l)(wtmp(plumeIndex1)-w_mean1(k,l))^2
538	hf1tmp(k,l) = TOTAL((wtmp(plumeIndex1)-w_mean1(k,l))*(ttmp(plumeIndex1)-tplume1moy(k,l)),1) / float(n_elements(plumeIndex1))
539	hf1tmpenv(k,l) = TOTAL((wtmp(envIndex1)-w_mean1_env(k,l))*(ttmp(envIndex1)-tenv1moy(k,l)),1) / float(n_elements(envIndex1))
540	if (envIndex1_ude(0) ne -1) then 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)) else hf1tmpenv_ude(k,l) =0.
541	if (downdraft_index1(0) ne -1) then 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)) else hf1tmp_down(k,l)=0.
542	IF((n_elements(plumeIndex1) + n_elements(envIndex1)) ne float(nx*ny)) then print, 'WARNING : INDEX PROBLEM : plume / env : ', n_elements(plumeIndex1), n_elements(envIndex1)
543	; 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)
544	ENDELSE
545	ENDFOR
546
547	Gamma_1(,l) = -(1./(alpha1out(,l)rhomoy1(,l)))deriv(altitudes_LES,Gamma_1_tmp(,l))
548	drhoahfdztmp = deriv(altitudes_LES,rhomoy1(,l)alpha1out(,l)hf1tmp(*,l))
549	drhoahfdztmpDetr = deriv(altitudes_LES,rhomoy1(,l)(1.-alpha1out(,l))hf1tmpenv(*,l))
550	drhoahfdztmpDetr_ude = deriv(altitudes_LES,rhomoy1(,l)(1.-alpha1out(,l)-beta1out(,l))hf1tmpenv_ude(,l))
551
552	wtmp=0.
553	ttmp=0.
554
555	FOR k=0,nz-1 DO BEGIN
556	IF(plumeIndex1out(0,k) eq -1) THEN BEGIN
557	e1_term2(k,l)=0.
558	e1_term2_ude(k,l)=0.
559	ENDIF ELSE BEGIN
560	e1_term2(k,l)=drhoahfdztmp(k)/(tenv1moy(k,l)-tplume1moy(k,l))
561	e1_term2_ude(k,l)=drhoahfdztmp(k)/(tenv1moy_ude(k,l)-tplume1moy(k,l))
562	ENDELSE
563
564	IF(envIndex1out(0,k) eq -1) THEN BEGIN
565	d1_term2(k,l)=0.
566	d1_term2_ude(k,l)=0.
567	ENDIF ELSE BEGIN
568	d1_term2(k,l)=-drhoahfdztmpDetr(k)/(tenv1moy(k,l)-tplume1moy(k,l))
569	d1_term2_ude(k,l)=-drhoahfdztmpDetr_ude(k)/(tenv1moy_ude(k,l)-tplume1moy(k,l))
570	ENDELSE
571	ENDFOR
572
573
574	tfull1=0.
575
576	zqtrac1=0.
577	zqtrac2 = getget(filesWRF(loop), s_trac2, count=[0,0,0,1], offset=[0,0,0,loop2])
578	FOR k=0,nz-1 DO BEGIN
579	anomalqtrac2(,,k) = zqtrac2(,,k) - TOTAL(TOTAL(REFORM(zqtrac2(,,k)),1),1)/ float(nx) / float(ny)
580	sigmazqtrac2(k) = STDDEV(zqtrac2(,,k))
581	IF (k ne 0) THEN BEGIN
582	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
583	sigmazminqtrac2(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac2(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
584	ENDIF ELSE BEGIN
585	sigmazminqtrac2(k) = sigmazqtrac2(k)
586	ENDELSE
587	plumeIndex2 = WHERE((anomalqtrac2(,,k) GT scaleMAX([sigmazqtrac2(k),sigmazminqtrac2(k)])) AND ((wprime(,*,k)+anomalw(k)) GT 0.))
588	envIndex2 = WHERE((anomalqtrac2(,,k) LE scaleMAX([sigmazqtrac2(k),sigmazminqtrac2(k)])) OR ((wprime(,*,k)+anomalw(k)) LE 0.))
589	IF(plumeIndex2(0) EQ -1) THEN BEGIN
590	fm_trac2_les(k,l)=0.
591	e_trac2_les(k,l)=0.
592	alpha2out(k,l)=0.
593	buoyancy2_les(k,l)=0.
594	w_mean2(k,l)=0.
595	ENDIF ELSE BEGIN
596	alpha2(k) = n_elements(plumeIndex2) / float(nx) / float(ny)
597	wprimetmp = reform(reform((anomalw(k)+wprime(,,k))),[nx*ny,1])
598	w_mean2(k,l) = mean(wprimetmp(plumeIndex2))
599	wprimetmp=0.
600	fm_trac2_les(k,l) = alpha2(k)rhomoy1(k,l)w_mean2(k,l)
601	tprimetmp = reform(reform(-tprime(,,k)),[nx*ny,1])
602	dtempdztmplin = reform(reform(dtempdztmp(,,k)),[nx*ny,1])
603	e_trac2_les(k,l) = TOTAL((1./(temporary(tprimetmp(plumeIndex2))))*(temporary(dtempdztmplin(plumeIndex2))),1)/n_elements(plumeIndex2)
604	alpha2out(k,l)=alpha2(k)
605	tfull=reform(tprime(,,k)+t(k,l),[nx*ny,1])
606	tplume2moy=mean(tfull(plumeIndex2))
607	tenv2moy=mean(tfull(envIndex2))
608	buoyancy2_les(k,l)=grav*(tplume2moy/tenv2moy-1.)
609	ENDELSE
610	ENDFOR
611	zqtrac2=0.
612
613	ENDIF
614
615	wt(,l) = TOTAL(TOTAL(TEMPORARY(tprime)TEMPORARY(wprime),1),1) / float(nx) / float(ny)
616	; wmax(l) = max(w_mean1(*,l))
617	l=l+1
618	ENDFOR
619	print, 'file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s'
620
621	ENDFOR
622
623	IF (got_updrafts EQ 'true') THEN BEGIN
624
625
626	hf1_term1 = hf1tmp*alpha1out
627	hf1_term2 = temporary(hf1tmpenv)*(1.-alpha1out)
628	hf1_term3 = alpha1out(1.-alpha1out)(w_mean1 - temporary(w_mean1_env))*(tplume1moy - tenv1moy)
629
630	hf1_ude_term1 = temporary(hf1tmp)*alpha1out
631	hf1_ude_term2 = temporary(hf1tmp_down)*beta1out
632	hf1_ude_term3 = temporary(hf1tmpenv_ude)*(1.-(alpha1out+beta1out))
633	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)
634
635	FOR k=0, nz-1 DO BEGIN
636	; dteta1moydt_entr(k,) = deriv(localtime,tplume1moy(k,))/3700. - dTeta_phys(k,*)
637	; dteta1moydt_detr(k,) = deriv(localtime,tplume1moy(k,))/3700. + dTeta_phys(k,*)
638	dteta1moydt_entr(k,) = deriv(localtime,tplume1moy(k,))/3700. - dTeta_phys(k,*)
639	dteta1moydt_detr(k,) = dTeta_phys(k,) - deriv(localtime,tplume1moy(k,*))/3700.
640	ENDFOR
641
642	FOR k=0, nz-1 DO BEGIN
643	FOR l=0, nttot-1 DO BEGIN
644	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.
645	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
646	; 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.
647	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.
648	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.
649	ENDFOR
650	ENDFOR
651
652	ENDIF
653
654	ht = TEMPORARY(pht) - hgtu/1000.
655	save, tsurf_les, w_mean1_env, 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
656
657	nz = n_elements(altitudes_LES)
658
659	ENDIF ELSE BEGIN
660
661	print, 'OK, file is here'
662	restore, filename=les_path+'/'+datname
663	nz = n_elements(altitudes_LES)
664	nttot = n_elements(tmoy_full(0,*))
665
666	OPENR, 23, les_path+'/input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
667
668	ENDELSE
669
670	tenv1moy = tplume1moy/((buoyancy1_les/grav)+1.)
671
672	taverage = string((localtime(nstot)-localtime(1))*3700./60.)
673	print, ''
674	print, ' -- Loading testphys1d data -- '
675
676	file1=gcm_path+'/diagfi.nc'
677	file2=gcm_convadj_path+'/diagfi.nc'
678	file3='/san0/acolmd/SIMUS/GCM3D_TestBed/diagfi.nc'
679
680	getcdf, file=file1, charvar='q2', invar=tke_gcm
681	getcdf, file=file1, charvar='aps', invar=aps
682	getcdf, file=file1, charvar='bps', invar=bps
683	getcdf, file=file1, charvar='co2col', invar=co2_col
684	;getcdf, file=file1, charvar='arcol', invar=ar_col
685	;getcdf, file=file1, charvar='ar', invar=ar
686	getcdf, file=file1, charvar='heatFlux_up', invar=heatFlux_up
687	getcdf, file=file1, charvar='heatFlux_down', invar=heatFlux_down
688	getcdf, file=file1, charvar='pplay', invar=pplay
689	getcdf, file=file1, charvar='pplev', invar=pplev
690	getcdf, file=file1, charvar='temp', invar=temp_gcm
691	getcdf, file=file1, charvar='zw2', invar=zw2_lev
692	getcdf, file=file1, charvar='fm_therm', invar=fm_therm_gcm_lev
693	getcdf, file=file1, charvar='entr_therm', invar=zdz_entr_therm_gcm
694	getcdf, file=file1, charvar='detr_therm', invar=zdz_detr_therm_gcm
695	getcdf, file=file1, charvar='fraca', invar=alpha_gcm_lev
696	getcdf, file=file1, charvar='buoyancyOut', invar=buoyancy_gcm
697	getcdf, file=file1, charvar='buoyancyEst', invar=buoyancy_est_gcm
698	getcdf, file=file1, charvar='zkh', invar=zkh
699	getcdf, file=file1, charvar='zh', invar=zh
700	getcdf, file=file1, charvar='tsurf', invar=tsurf_gcm
701	;getcdf, file=file1, charvar='zmax', invar=zi_gcm
702	getcdf, file=file1, charvar='lmax_th', invar=lmax_gcm
703	getcdf, file=file1, charvar='hfmax_th', invar=hfmax_th1d
704	getcdf, file=file1, charvar='wmax_th', invar=wmax_th1d
705
706	if (overplot_convadj eq 'true') then begin
707	getcdf, file=file2, charvar='temp', invar=temp_gcm_convadj
708	getcdf, file=file2, charvar='pplay', invar=pplay_convadj
709	endif
710	if (plot_3d eq 'true') then begin
711	getcdf, file=file3, charvar='temp', invar=temp_gcm_3d
712	getcdf, file=file3, charvar='pplay', invar=pplay_3d
713	getcdf, file=file3, charvar='latitude', invar=latitude_3d
714	getcdf, file=file3, charvar='longitude', invar=longitude_3d
715	nWEmx_3d = n_elements(reform(temp_gcm_3d(*,0,0,0)))
716	nNSmx_3d = n_elements(reform(temp_gcm_3d(0,*,0,0)))
717	nZmx_3d = n_elements(reform(temp_gcm_3d(0,0,*,0)))
718	nTmx_3d = n_elements(reform(temp_gcm_3d(0,0,0,*)))
719	ndays_3d = 1.
720	lctu_gcm_3d = 0.
721	history_interval_s_gcm_3d = ndays_3d*88800./float(nTmx_3d) ; Timestep interval of gcm 1d simu in sec
722	localtime_lon0 = lctu_gcm_3d + history_interval_s_gcm_3d*findgen(nTmx_3d)/3700.
723	Xc = 205.
724	Yc = 21.8
725	plot_index_x = (Xc-longitude_3d(0))/(longitude_3d(1)-longitude_3d(0))
726	plot_index_y = (Yc-latitude_3d(0))/(latitude_3d(1)-latitude_3d(0))
727	localtime_true = localtime_lon0 -(12./180.)*Xc
728	endif
729
730
731	nTmx = n_elements(reform(temp_gcm(0,*)))
732	if (overplot_convadj eq 'true') then begin
733	nTmx_convadj = n_elements(reform(temp_gcm_convadj(0,*)))
734	endif else begin
735	nTmx_convadj = 10000.
736	endelse
737	ndays = 1.
738	print, ''
739	print, 'WARNING ----------------------- '
740	print, 'CONFIGURATION : '+string(ndays,format='(I0)')+' days simulation'
741	print, ''
742	history_interval_s_gcm = ndays*88800./float(nTmx) ; Timestep interval of gcm 1d simu in sec
743	history_interval_s_gcm_convadj = ndays*88800./float(nTmx_convadj)
744	localtime_gcm = lctu_gcm + history_interval_s_gcm*findgen(nTmx)/3700.
745	localtime_gcm_convadj = lctu_gcm + history_interval_s_gcm_convadj*findgen(nTmx_convadj)/3700.
746	; **********************************
747	; ****** PLOTS ****************
748
749	if (f_offset eq 'true') then begin
750	offset_localtime = 3.108100
751	lt_plot=11.
752	endif else begin
753	offset_localtime = 0.
754	lt_plot=12.
755	endelse
756	localtime=localtime+offset_localtime
757	;localtime_gcm=localtime_gcm+history_interval_s_gcm/3700. ; we add the offset from the fact that we output (non sense :) )
758	localtime_gcm=localtime_gcm
759	; par contre il faut prendre en compte le fait que la premiere frame du les est decalee de 1 !
760	localtime=localtime
761
762	print, '****************************************************************************************************'
763	print, 'local time LES'
764	print, localtime
765	print, '****************************************************************************************************'
766	print, 'local time GCM'
767	print, localtime_gcm
768	print, '****************************************************************************************************'
769
770	time_offset = (ndays-1.)*24.
771
772	print, 'TIME STEP LES : ',(localtime(1)-localtime(0))*3700.
773
774
775	lt_plot_ini = 6.
776	lt_plotindex_les_ini = where(localtime eq lt_plot_ini)
777	lt_plotindex_gcm_ini = where(localtime_gcm eq (lt_plot_ini+time_offset))
778	;lt_plotindex_gcm_ini = where(localtime_gcm eq (lt_plot_ini+time_offset+history_interval_s_gcm/3700.))
779
780	lt_plot0 = 10.
781	lt_plotindex_les0 = where(localtime eq lt_plot0)
782	lt_plotindex_gcm0 = where(localtime_gcm eq (lt_plot0+time_offset))
783	lt_plotindex_gcm_convadj0 = where(localtime_gcm_convadj eq (lt_plot0+time_offset))
784
785	lt_plot0a = 11.
786	lt_plotindex_les0a = where(localtime eq lt_plot0a)
787	lt_plotindex_gcm0a = where(localtime_gcm eq (lt_plot0a+time_offset))
788
789	lt_plotindex_les = where((localtime lt lt_plot+0.01) and (localtime gt lt_plot-0.01))
790	lt_plotindex_gcm = where(localtime_gcm eq (lt_plot+time_offset))
791	lt_plotindex_gcm_convadj = where(localtime_gcm_convadj eq (lt_plot+time_offset))
792	print, 'lt plotindex les 12h'
793	print, lt_plotindex_les
794
795	lt_plota = 13.
796	lt_plotindex_lesa = where(localtime eq lt_plota)
797	lt_plotindex_gcma = where(localtime_gcm eq (lt_plota+time_offset))
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_plot2a = 15.
805	lt_plotindex_les2a = where(localtime eq lt_plot2a)
806	lt_plotindex_gcm2a = where(localtime_gcm eq (lt_plot2a+time_offset))
807
808	lt_plot3 = 16.
809	lt_plotindex_les3 = where(localtime eq lt_plot3)
810	lt_plotindex_gcm3 = where(localtime_gcm eq (lt_plot3+time_offset))
811	lt_plotindex_gcm_convadj3 = where(localtime_gcm_convadj eq (lt_plot3+time_offset))
812
813	lt_plot3a = 17.
814	lt_plotindex_les3a = where(localtime eq lt_plot3a)
815	lt_plotindex_gcm3a = where(localtime_gcm eq (lt_plot3a+time_offset))
816
817	lt_plot4 = 18.
818	lt_plotindex_les4 = where(localtime eq lt_plot4)
819	lt_plotindex_gcm4 = where(localtime_gcm eq (lt_plot4+time_offset))
820	lt_plotindex_gcm_convadj4 = where(localtime_gcm_convadj eq (lt_plot4+time_offset))
821	;--------------------------------------------------------------------------------
822	;---------------------------------------------------------------------------------
823
824	nTmx_les=n_elements(reform(wt(0,*)))
825	nZmx=n_elements(aps) ; number of vertical levels
826	H_low=9650. ; scale height at low altitudes
827	H_high=15000. ; scale height at high altitudes
828	trans_window=10. ; # of levels over which H(:) goes from H_low to H_high
829	lev_trans=32.+trans_window/2. ; level at which H(lev_trans)=(H_low+H_high)/2
830	P_ref=p0 ; reference surface pressure used to build zsurface -610 Pa-
831	Hgcm = make_array(nZmx)
832	altitudes_GCM = make_array(nZmx)
833	; Build scale heights
834	;FOR k=0,nZmx-1 DO BEGIN
835	; Hgcm(k)=H_low+(H_high-H_low)0.5(1.0+tanh(6.*(k-lev_trans)/trans_window))
836	;ENDFOR
837
838	FOR k=0,nZmx-1 DO BEGIN
839	Hgcm(k)=R*temp_gcm(k,lt_plotindex_gcm)/grav
840	ENDFOR
841	print, 'Hgcm'
842	print, Hgcm
843	; Compute altitudes_GCM
844	FOR k=0,nZmx-1 DO BEGIN
845	altitudes_GCM(k)=-Hgcm(k)*alog(pplay(k,lt_plotindex_gcm)/pGround)
846	ENDFOR
847	Hgcm=0.
848
849	teta_gcm = temp_gcm * (p0/pplay)^r_cp
850	if (overplot_convadj eq 'true') then begin
851	teta_gcm_convadj = temp_gcm_convadj * (p0/pplay_convadj)^r_cp
852	endif
853
854	OPENR, 1, gcm_path+'/profile'
855	data=FLTARR(nZmx+1)
856	READF, 1, data
857	temp_gcm_0_ground = data(0)
858	temp_gcm_0 = data(1:nZmx-1)
859	data = 0.
860	CLOSE, 1
861
862	teta_gcm_0 = temp_gcm_0 * (p0/pplay)^r_cp
863	approx_zdz_gcm = make_array(nZmx)
864	approx_zdz_gcm(0)=altitudes_GCM(1)
865	FOR k=1, nZmx-2 DO BEGIN
866	approx_zdz_gcm(k) = altitudes_GCM(k+1) - altitudes_GCM(k)
867	ENDFOR
868	approx_zdz_gcm(nZmx-1)=approx_zdz_gcm(nZmx-2)
869
870	print, 'approx zdz gcm'
871	print, approx_zdz_gcm
872
873	print, '****************************************************************************************************'
874	print, 'altitudes LES based on phtot : inter-levels'
875	print, altitudes_LES
876	print, '****************************************************************************************************'
877	print, 'altitudes GCM based on pplay : inter-levels'
878	print, altitudes_GCM
879	print, '****************************************************************************************************'
880
881	; Compute tracer deviation :
882
883	co2_col = co2_col/co2_col(0)
884	;ar_col = ar_col/ar_col(0)
885	;tke_col = tke_col+1.
886
887	; Compute <teta> les
888
889	teta_les = temporary(t)
890
891	rho = pt/(R*temp_les)
892
893	;print, 'bidouille'
894	;
895	;FOR l=0,nTmx -1 DO BEGIn
896	;print, (1300.*hfmax_th1d(l)/(TOTAL(temp_gcm(0:lmax_gcm(l),l),1)/(lmax_gcm(l)+1.)))/wmax_th1d(l)
897	;ENDFOR
898
899
900	; ========================================================================
901	; ========================================================================
902
903	IF (visualization_mode eq 'true') THEN BEGIN
904
905	print,' *****************************************-----------------------------------'
906	print,' ********** PLUME ********************-----------------------------------'
907	print,' *****************************************-----------------------------------'
908
909	; We are evaluating the first time-step element of the file number 'loop-1' :
910	; file 1 starts at 8h (loop =0, loop2 =0)
911	; file 6 starts at 13h (roughly) (loop =5, loops2=0)
912	; file 12 starts at 18h (roughly) (loop = 11,loop2 = 0)
913
914	loop=uint(loop_special)-1
915	;loop2=34
916	loop2=10
917	domain='d01'
918	filesWRF = FindFile(les_path+'/wrfout_'+domain+'_????-??-??_??:??:??')
919	anomalw=1.
920
921	zqtrac1 = dblarr(nx,ny,nz) & zqtrac2 = dblarr(nx,ny,nz)
922	sigmazqtrac1 = fltarr(nz) & sigmazqtrac2 = fltarr(nz)
923	sigmazminqtrac1 = fltarr(nz) & sigmazminqtrac2 = fltarr(nz)
924	anomalqtrac1 = fltarr(nx,ny,nz) & anomalqtrac2 = fltarr(nx,ny,nz)
925	zqtrac1 = getget(filesWRF(loop), 'qtrac1', count=[0,0,0,1], offset=[0,0,0,loop2])
926	zqtrac2 = getget(filesWRF(loop), 'qtrac2', count=[0,0,0,1], offset=[0,0,0,loop2])
927	wprime = getget(filesWRF(loop), 'W', anomaly = anomalw, count=[0,0,0,1], offset=[0,0,0,loop2])
928	supermask1 = fltarr(nx,ny,nz)
929	supermask2 = fltarr(nx,ny,nz)
930	k_out_histo = 8
931	k_out_hist = [1,10,25,50,70,85]
932	nbtest=n_elements(k_out_hist)
933	b=0
934	FOR k=0,nz-1 DO BEGIN
935	mask1=fltarr(nx*ny)
936	mask2=fltarr(nx*ny)
937	anomalqtrac1(,,k) = zqtrac1(,,k) - TOTAL(TOTAL(REFORM(zqtrac1(,,k)),1),1)/ float(nx) / float(ny)
938	sigmazqtrac1(k) = STDDEV(REFORM(zqtrac1(,,k)))
939	IF (k ne 0) THEN BEGIN
940	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
941	sigmazminqtrac1(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac1(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
942	ENDIF ELSE BEGIN
943	sigmazminqtrac1(k) = sigmazqtrac1(k)
944	ENDELSE
945	print, sigmazqtrac1(k),sigmazminqtrac1(k)
946	plumeIndex1 = WHERE((anomalqtrac1(,,k) GT scaleMAX([sigmazqtrac1(k),sigmazminqtrac1(k)])) AND ((anomalw(k)+wprime(,*,k)) GT 0.))
947	anomalqtrac2(,,k) = zqtrac2(,,k) - TOTAL(TOTAL(REFORM(zqtrac2(,,k)),1),1)/ float(nx) / float(ny)
948	sigmazqtrac2(k) = STDDEV(REFORM(zqtrac2(,,k)))
949	IF (k ne 0) THEN BEGIN
950	subsampledAltitudes = INTERPOL(altitudes_LES(0:k),findgen(k+1),findgen(decimate*k+1)/decimate)
951	sigmazminqtrac2(k) = (sigmao/(altitudes_LES(k)-altitudes_LES(0)))*INT_TABULATED(subsampledAltitudes,INTERPOL(sigmazqtrac2(0:k),altitudes_LES(0:k),subsampledAltitudes),/DOUBLE)
952	ENDIF ELSE BEGIN
953	sigmazminqtrac2(k) = sigmazqtrac2(k)
954	ENDELSE
955	plumeIndex2 = WHERE((anomalqtrac2(,,k) GT scaleMAX([sigmazqtrac2(k),sigmazminqtrac2(k)])) AND ((anomalw(k)+wprime(,*,k)) GT 0.))
956	IF(plumeIndex1(0) NE -1 ) THEN BEGIN
957	mask1(plumeIndex1) = 1.
958	ENDIF ELSE BEGIN
959	mask1(*)=0.
960	ENDELSE
961	IF(plumeIndex2(0) NE -1 ) THEN BEGIN
962	mask2(plumeIndex2) = 1.
963	ENDIF ELSE BEGIN
964	mask2(*)=0.
965	ENDELSE
966	mask1 = reform(mask1,[nx,ny])
967	supermask1(,,k) = mask1(,)
968	mask2 = reform(mask2,[nx,ny])
969	supermask2(,,k) = mask2(,)
970	; IF (k eq k_out_histo) THEN BEGIN
971	; plume1_out = plumeIndex1
972	; plume2_out = plumeIndex2
973	; ENDIF
974
975	IF (k eq k_out_hist(0)) THEN BEGIN
976	c1=plumeIndex1
977	c2=plumeIndex2
978	ENDIF
979	IF (k eq k_out_hist(1)) THEN BEGIN
980	d1=plumeIndex1
981	d2=plumeIndex2
982	ENDIF
983	IF (k eq k_out_hist(2)) THEN BEGIN
984	e1=plumeIndex1
985	e2=plumeIndex2
986	ENDIF
987	IF (k eq k_out_hist(3)) THEN BEGIN
988	f1=plumeIndex1
989	f2=plumeIndex2
990	ENDIF
991	IF (k eq k_out_hist(4)) THEN BEGIN
992	g1=plumeIndex1
993	g2=plumeIndex2
994	ENDIF
995	IF (k eq k_out_hist(5)) THEN BEGIN
996	h1=plumeIndex1
997	h2=plumeIndex2
998	ENDIF
999	ENDFOR
1000
1001	IF (Histo eq 'false') THEN BEGIN
1002	IVOLUME, supermask1
1003	IVOLUME, supermask2
1004	ENDIF ELSE BEGIN
1005
1006	; -------------------------------------------------------------------------------
1007	; THIS IS THE ULTRA-GORE UN-ESTHETIC UGLY-AS-HELL LOOP FOR CONCENTRATION PLOTTING
1008	; but well, this is just because idl cannot handle arrays as well as I would like
1009	; -------------------------------------------------------------------------------
1010	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Distrib.ps'
1011	PSOPEN, THICK=100, CHARSIZE=60, FILE = filename, FONT = 5, TFONT = 5,XPLOTS=2,YPLOTS=3,MARGIN=2000
1012
1013	FOR n=0, nbtest-1 DO BEGIN
1014
1015	CASE n OF
1016	0:BEGIN
1017	plume1_out=c1 & plume2_out=c2
1018	END
1019	1:BEGIN
1020	plume1_out=d1 & plume2_out=d2
1021	END
1022	2:BEGIN
1023	plume1_out=e1 & plume2_out=e2
1024	END
1025	3:BEGIN
1026	plume1_out=f1 & plume2_out=f2
1027	END
1028	4:BEGIN
1029	plume1_out=g1 & plume2_out=g2
1030	END
1031	5:BEGIN
1032	plume1_out=h1 & plume2_out=h2
1033	END
1034	ENDCASE
1035
1036	ToBin1 = reform(zqtrac1(,,k_out_hist(n)),[nx*ny,1])
1037	ToBin2 = reform(zqtrac2(,,k_out_hist(n)),[nx*ny,1])
1038	svmin1=min(ToBin1) & svmin2=min(ToBin2)
1039	svmax1=max(ToBin1) & svmax2=max(ToBin2)
1040	NBINS=100
1041	ds1=(svmax1-svmin1+1.)/(NBINS-1) & ds2=(svmax2-svmin2+1.)/(NBINS-1)
1042	Xaxis1 = svmin1+((svmax1-svmin1)/(NBINS-1))*indgen(NBINS)
1043	Xaxis2 = svmin2+((svmax2-svmin2)/(NBINS-1))*indgen(NBINS)
1044	Bin1=HISTOGRAM(ToBin1,nbins=NBINS)
1045	Bin2=INTERPOL(HISTOGRAM(ToBin2,nbins=NBINS),Xaxis2,Xaxis1,/SPLINE)
1046
1047	what_I_plot = [[Bin1],[Bin2]]
1048	labels=['LES tracer 1 conc. distrib.','LES tracer 2 conc. distrib.']
1049
1050	title_user = TestCase+SubCase+' LES tracer 1&2 concentration distribution at '+string(altitudes_LES(k_out_hist(n)))+'m AGL'
1051	IF (n lt 3) THEN BEGIN
1052	POS, XPOS=1, YPOS=uint(n+1)
1053	ENDIF ELSE BEGIN
1054	POS, XPOS=2, YPOS=uint(n+1)-3
1055	ENDELSE
1056	MAP
1057	CS, SCALE=28
1058	GSET, XMIN=0, XMAX=20, YMIN=0, YMAX=300, TITLE=title_user
1059	cols=INDGEN(2)+2
1060	GPLOT, X=Xaxis1, Y=what_I_plot, /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
1061	AXES, XSTEP = 4, XTITLE='Trac concentration (kg/kg)', YSTEP=50, YTITLE='Counts',NDECS=1
1062
1063	ToBin1 = reform(zqtrac1(,,k_out_hist(n)),[nx*ny,1])
1064	ToBin2 = reform(zqtrac2(,,k_out_hist(n)),[nx*ny,1])
1065	ToBin1 = ToBin1(plume1_out)
1066	ToBin2 = ToBin2(plume2_out)
1067	svmin1=min(ToBin1) & svmin2=min(ToBin2)
1068	svmax1=max(ToBin1) & svmax2=max(ToBin2)
1069	NBINS=50
1070	ds1=(svmax1-svmin1+1.)/(NBINS-1) & ds2=(svmax2-svmin2+1.)/(NBINS-1)
1071	Xaxis1 = svmin1+((svmax1-svmin1)/(NBINS-1))*indgen(NBINS)
1072	Xaxis2 = svmin2+((svmax2-svmin2)/(NBINS-1))*indgen(NBINS)
1073	Bin1=HISTOGRAM(ToBin1,nbins=NBINS)
1074	Bin2=HISTOGRAM(ToBin2,nbins=NBINS)
1075	oplot, Xaxis1, Bin1, psym=4
1076	oplot, Xaxis2, Bin2, psym=5
1077
1078	ENDFOR
1079
1080	PSCLOSE, /NOVIEW
1081	spawn, 'ps2png '+filename
1082
1083	ENDELSE
1084
1085	ENDIF ELSE BEGIN
1086
1087	; =========================================================================
1088	; =========================================================================
1089	IF (got_updrafts EQ 'true') THEN BEGIN
1090
1091	print, '........ ALPHA'
1092
1093	alpha_interlay_gcm = make_array(nZmx)
1094	FOR k=0, nZmx-2 DO BEGIN
1095	alpha_interlay_gcm(k) = (alpha_gcm_lev(k,lt_plotindex_gcm)+alpha_gcm_lev(k+1,lt_plotindex_gcm))/2.
1096	ENDFOR
1097	alpha_interlay_gcm(nZmx-1)=0.
1098
1099	smoothed_alpha1_les = make_array(nz)
1100	smoothed_alpha2_les = make_array(nz)
1101	smoothed_beta1_les = make_array(nz)
1102	FOR t=-ns,ns DO BEGIN
1103	smoothed_alpha1_les = smoothed_alpha1_les + REFORM(alpha1out(*,lt_plotindex_les+t))
1104	smoothed_alpha2_les = smoothed_alpha2_les + REFORM(alpha2out(*,lt_plotindex_les+t))
1105	smoothed_beta1_les = smoothed_beta1_les + REFORM(beta1out(*,lt_plotindex_les+t))
1106	ENDFOR
1107
1108	smoothed_alpha1_les = smoothed_alpha1_les/nstot
1109	smoothed_alpha2_les = smoothed_alpha2_les/nstot
1110	smoothed_beta1_les = smoothed_beta1_les/nstot
1111
1112	ENDIF
1113	; =========================================================================
1114
1115	; * Temperatures *
1116
1117	if (f_offset eq 'false') then begin
1118	print, '........ TEMPERATURES'
1119
1120	xmin = 190
1121	xmax = 250
1122	if (TestCase eq 'Case_Z') then begin
1123	xmin = 170
1124	xmax = 250
1125	endif
1126	if (TestCase eq 'Case_C') then begin
1127	xmin = 180
1128	xmax = 240
1129	endif
1130
1131
1132	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))]]
1133	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)]
1134	title_user = TestCase+SubCase+LayerCase+' Temperatures Comparison'
1135	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Temperature.ps'
1136	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1137	CS, SCALE=28
1138	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=12, TITLE=title_user
1139	cols=INDGEN(6)+2
1140	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1141	AXES, XSTEP = 5, XTITLE='Temperature (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
1142
1143	oplot, temp_les(*,lt_plotindex_les_ini), altitudes_LES/1000., psym=4
1144	oplot, temp_les(*,lt_plotindex_les0), altitudes_LES/1000., psym=4
1145	oplot, temp_les(*,lt_plotindex_les), altitudes_LES/1000., psym=4
1146	oplot, temp_les(*,lt_plotindex_les2), altitudes_LES/1000., psym=4
1147	oplot, temp_les(*,lt_plotindex_les3), altitudes_LES/1000., psym=4
1148	oplot, temp_les(*,lt_plotindex_les4), altitudes_LES/1000., psym=4
1149
1150	if(overplot_convadj EQ 'true') then begin
1151	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj0), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1152	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1153	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj2), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1154	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj3), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1155	oplot, temp_gcm_convadj(*,lt_plotindex_gcm_convadj4), altitudes_GCM/1000., thick=0.1,color=8,linestyle=3
1156	endif
1157
1158
1159	PSCLOSE, /NOVIEW
1160
1161	spawn, 'ps2png '+filename
1162
1163	endif
1164
1165	; * Pressions *
1166
1167	print, '........ PRESSURES'
1168
1169	;what_I_plot = make_array(nZmx)
1170	;labels=['TH P 1d, lt='+string(lt_plot)]
1171	;title_user = TestCase+SubCase+' Pressure Comparisons'
1172	;filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Pressure.ps'
1173	;PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1174	;CS, SCALE=28
1175	;GSET, XMIN=600, XMAX=900, YMIN=0, YMAX=0.5, TITLE=title_user
1176	;cols=INDGEN(1)+2
1177	;GPLOT, X=pplay(,lt_plotindex_gcm), Y=-alog(pplay(,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1178	;AXES, XSTEP = 25, XTITLE='Log P', YSTEP=0.1, YTITLE='Altitude (km)',NDECS=1
1179	;
1180	;oplot, pt(,lt_plotindex_les),-alog(pt(,lt_plotindex_les)/pGround), psym=4
1181	;
1182	;PSCLOSE, /NOVIEW
1183	;
1184	;spawn, 'ps2png '+filename
1185
1186	what_I_plot = make_array(nZmx)
1187	labels=['TH P 1d, lt='+string(lt_plot)]
1188	title_user = TestCase+SubCase+LayerCase+' Pressure Comparisons'
1189	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Pressure.ps'
1190	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1191	CS, SCALE=28
1192	GSET, XMIN=400, XMAX=900, YMIN=0, YMAX=10, TITLE=title_user
1193	cols=INDGEN(1)+2
1194	GPLOT, X=pplay(*,lt_plotindex_gcm), Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1195	AXES, XSTEP = 25, XTITLE='Pression (Pa)', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
1196
1197	oplot, pt(*,lt_plotindex_les),altitudes_LES/1000., psym=4
1198
1199	PSCLOSE, /NOVIEW
1200
1201	spawn, 'ps2png '+filename
1202
1203
1204	; * Temperatures potentielles *
1205	if(full eq 'true') then begin
1206
1207	xmin = 210
1208	xmax = 230
1209	if (TestCase eq 'Case_C') then begin
1210	xmin = 225
1211	xmax = 255
1212	endif
1213	if (TestCase eq 'Case_I') then begin
1214	xmin = 195
1215	xmax = 250
1216	endif
1217	if (TestCase eq 'Case_Z') then begin
1218	xmin = 230
1219	xmax = 290
1220	endif
1221	if (TestCase eq 'ExtremeCase') then begin
1222	xmin = 225
1223	xmax = 255
1224	endif
1225	if (TestCase eq 'Exomars') then begin
1226	xmin = 225
1227	xmax = 255
1228	endif
1229
1230	if (f_offset eq 'false') then begin
1231
1232	print, '........ POTENTIAL TEMPERATURES'
1233	;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))]]
1234	what_I_plot = [[reform(teta_gcm(,lt_plotindex_gcm))],[reform(teta_gcm(,lt_plotindex_gcm2))],[reform(teta_gcm(*,lt_plotindex_gcm3))]]
1235	labels=['TH teta 1d, lt='+string(lt_plot),'TH teta 1d, lt='+string(lt_plot2),'TH teta 1d, lt='+string(lt_plot3)]
1236	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1237	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta.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.6, TITLE=title_user
1241	cols=INDGEN(3)+2
1242	GPLOT, X=what_I_plot, Y=-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.2, YTITLE='-Log(P/P0) ',NDECS=1
1244
1245	;oplot, teta_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4,SYMSIZE=0.5, thick=0.5
1246	;oplot, teta_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4,SYMSIZE=0.5, thick=0.5
1247	oplot, teta_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4,SYMSIZE=0.5, thick=0.5
1248	oplot, teta_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4,SYMSIZE=0.5, thick=0.5
1249	oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4,SYMSIZE=0.5, thick=0.5
1250	;oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4,SYMSIZE=0.5, thick=0.5
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), thick=0.5
1259
1260	PSCLOSE, /NOVIEW
1261
1262	spawn, 'ps2png '+filename
1263
1264
1265	xmin = 210
1266	xmax = 225
1267
1268	if (TestCase eq 'Case_C') then begin
1269	xmin = 230
1270	xmax = 255
1271	endif
1272
1273	if (TestCase eq 'ExtremeCase') then begin
1274	xmin = 225
1275	xmax = 255
1276	endif
1277
1278	;xmin = 280
1279	;xmax = 300
1280	what_I_plot = [[tsurf_gcm(lt_plotindex_gcm),reform(teta_gcm(,lt_plotindex_gcm))],[tsurf_gcm(lt_plotindex_gcm2),reform(teta_gcm(,lt_plotindex_gcm2))]]
1281	labels=['TH teta 1d, lt='+string(lt_plot),'TH teta 1d, lt='+string(lt_plot2)]
1282	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1283	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta_zoom.ps'
1284	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1285	CS, SCALE=28
1286	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=0.05, TITLE=title_user
1287	cols=INDGEN(2)+2
1288	GPLOT, X=what_I_plot, Y=[0.,-alog(pplay(*,lt_plotindex_gcm)/pGround)], /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1289	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.005, YTITLE='-Log(P/P0) ',NDECS=1
1290
1291	;oplot, teta_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1292	;oplot, teta_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1293	oplot, teta_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1294	oplot, teta_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1295	;oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1296	;oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1297	if(overplot_convadj EQ 'true') then begin
1298	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj0), -alog(pplay_convadj(,lt_plotindex_gcm_convadj0)/pGround), thick=0.1,color=8,linestyle=3
1299	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj), -alog(pplay_convadj(,lt_plotindex_gcm_convadj)/pGround), thick=0.1,color=8,linestyle=3
1300	oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj2), -alog(pplay_convadj(,lt_plotindex_gcm_convadj2)/pGround), thick=0.1,color=8,linestyle=3
1301	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj3), -alog(pplay_convadj(,lt_plotindex_gcm_convadj3)/pGround), thick=0.1,color=8,linestyle=3
1302	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj4), -alog(pplay_convadj(,lt_plotindex_gcm_convadj4)/pGround), thick=0.1,color=8,linestyle=3
1303	endif
1304	;oplot, teta_gcm_0(), -alog(pplay(,lt_plotindex_gcm)/pGround)
1305
1306	PSCLOSE, /NOVIEW
1307
1308	spawn, 'ps2png '+filename
1309
1310
1311	xmin = 225
1312	xmax = 270
1313	;xmin = 280
1314	;xmax = 300
1315	if (TestCase eq 'Case_C') then begin
1316	xmin = 215
1317	xmax = 240
1318	endif
1319
1320	what_I_plot = [[tsurf_gcm(lt_plotindex_gcm),reform(temp_gcm(,lt_plotindex_gcm))],[tsurf_gcm(lt_plotindex_gcm2),reform(temp_gcm(,lt_plotindex_gcm2))]]
1321	ymax=0.05
1322	ystep=0.005
1323	print, 'tsurf gcm :',tsurf_gcm(lt_plotindex_gcm)
1324
1325	if (TestCase eq 'Exomars') then begin
1326	xmin = 235
1327	;xmax = 255
1328	xmax=295
1329	what_I_plot = [[tsurf_gcm(lt_plotindex_gcm),reform(temp_gcm(,lt_plotindex_gcm))],[10.tsurf_gcm(lt_plotindex_gcm2),10.reform(temp_gcm(,lt_plotindex_gcm2))]]
1330	;what_I_plot = 5.*what_I_plot
1331	;ymax=0.027
1332	ymax=0.003
1333	ystep=0.0005
1334	endif
1335
1336
1337
1338	labels=['TH T 1d, lt='+string(lt_plot),'TH T 1d, lt='+string(lt_plot2)]
1339	title_user = TestCase+SubCase+LayerCase+' T comparisons'
1340	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_T_zoom.ps'
1341	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1342	CS, SCALE=28
1343	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=ymax, TITLE=title_user
1344	cols=INDGEN(2)+2
1345	GPLOT, X=what_I_plot, Y=[0.,-alog(pplay(*,lt_plotindex_gcm)/pGround)], /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
1346	AXES, XSTEP = 5, XTITLE='Temperature (K)', YSTEP=ystep, YTITLE='-Log(P/P0) ',NDECS=4
1347
1348	;oplot, temp_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1349	;oplot, temp_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1350
1351	oplot, [tsurf_gcm(lt_plotindex_gcm),temp_les(,lt_plotindex_les)], [0.,-alog(pt(,lt_plotindex_les)/pGround)], psym=4
1352	oplot, temp_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), thick=0.1
1353	;oplot, temp_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1354	;oplot, temp_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1355
1356	;oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1357	;oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1358	if(overplot_convadj EQ 'true') then begin
1359	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj0), -alog(pplay_convadj(,lt_plotindex_gcm_convadj0)/pGround), thick=0.1,color=8,linestyle=3
1360
1361	oplot, temp_gcm_convadj(,lt_plotindex_gcm_convadj), -alog(pplay_convadj(,lt_plotindex_gcm_convadj)/pGround), thick=0.1,color=8,psym=2
1362	oplot, temp_gcm_convadj(,lt_plotindex_gcm_convadj), -alog(pplay_convadj(,lt_plotindex_gcm_convadj)/pGround), thick=0.1,color=8,linestyle=3
1363	;oplot, temp_gcm_convadj(,lt_plotindex_gcm_convadj2), -alog(pplay_convadj(,lt_plotindex_gcm_convadj2)/pGround), thick=0.1,color=8,linestyle=3
1364
1365
1366	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj3), -alog(pplay_convadj(,lt_plotindex_gcm_convadj3)/pGround), thick=0.1,color=8,linestyle=3
1367	;oplot, teta_gcm_convadj(,lt_plotindex_gcm_convadj4), -alog(pplay_convadj(,lt_plotindex_gcm_convadj4)/pGround), thick=0.1,color=8,linestyle=3
1368	endif
1369	;oplot, teta_gcm_0(), -alog(pplay(,lt_plotindex_gcm)/pGround)
1370
1371	PSCLOSE, /NOVIEW
1372
1373	spawn, 'ps2png '+filename
1374
1375
1376	endif
1377
1378	if (plot_3d eq 'true') then begin
1379	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))]]
1380	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)]
1381	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1382	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta.ps'
1383	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1384	CS, SCALE=28
1385	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=1.4, TITLE=title_user
1386	cols=INDGEN(6)+2
1387	GPLOT, X=what_I_plot, Y=-alog(pplay_3d(*,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1388	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.2, YTITLE='-Log(P/P0) ',NDECS=1
1389
1390	oplot, teta_les(,lt_plotindex_les_ini), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1391	oplot, teta_les(,lt_plotindex_les0), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1392	oplot, teta_les(,lt_plotindex_les), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1393	oplot, teta_les(,lt_plotindex_les2), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1394	oplot, teta_les(,lt_plotindex_les3), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1395	oplot, teta_les(,lt_plotindex_les4), -alog(pt(,lt_plotindex_les)/pGround), psym=4
1396
1397	PSCLOSE, /NOVIEW
1398
1399	spawn, 'ps2png '+filename
1400	endif
1401
1402	endif else begin
1403
1404	print, '........ POTENTIAL TEMPERATURES'
1405
1406	what_I_plot = reform(teta_gcm(*,lt_plotindex_gcm))
1407	labels=['TH teta 1d, lt='+string(lt_plot)]
1408	title_user = TestCase+SubCase+LayerCase+' Teta comparisons (recomputed from T and P)'
1409	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Teta.ps'
1410	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1411	CS, SCALE=28
1412	GSET, XMIN=210, XMAX=240, YMIN=0, YMAX=2, TITLE=title_user
1413	cols=INDGEN(1)+2
1414	GPLOT, X=what_I_plot, Y=-alog(pplay(*,lt_plotindex_gcm)/pGround), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1415	AXES, XSTEP = 5, XTITLE='Potential temperature (K)', YSTEP=0.2, YTITLE='-Log(P/P0) ',NDECS=1
1416
1417	oplot, teta_les(,lt_plotindex_les),-alog(pt(,lt_plotindex_les)/pGround), psym=4
1418	PSCLOSE, /NOVIEW
1419
1420	spawn, 'ps2png '+filename
1421	endelse
1422
1423	print, '........ SURFACE TEMPERATURES'
1424
1425
1426	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_tsurf_gcmsoil.nc', charvar='TSURF', invar=tsurf_les_tmp
1427	;tsurf_les=make_array(nttot)
1428	;FOR l=0,nttot-1 DO BEGIN
1429	; tsurf_les(l)=TOTAL(TOTAL(tsurf_les_tmp(,,l),1),1)/float(n_elements(reform(tsurf_les_tmp(,0,0)))) /float(n_elements(reform(tsurf_les_tmp(0,,0))))
1430	;ENDFOR
1431	;tsurf_les_tmp=0.
1432
1433	what_I_plot = tsurf_gcm
1434	labels=['TH 1d tsurf']
1435	title_user = TestCase+SubCase+LayerCase+' Surface temperatures (recomputed from T and P)'
1436	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_tsurf.ps'
1437	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1438	CS, SCALE=28
1439	GSET, XMIN=6, XMAX=18, YMIN=150, YMAX=320, TITLE=title_user
1440	cols=INDGEN(1)+2
1441	GPLOT, X=localtime_gcm, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1442	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=5., YTITLE='Surface temperature',NDECS=1
1443
1444	;oplot, localtime, tsurf, psym=4,thick=0.3 ;tsurf les with les soil
1445	;oplot, localtime, tsurf_les, psym=4,thick=0.3 ;tsurf les with gcm soil
1446
1447	PSCLOSE, /NOVIEW
1448
1449	spawn, 'ps2png '+filename
1450
1451
1452	if (got_swdownz eq 'true') then begin
1453
1454	print, '........ DOWNWARD SOLAR FLUX AT SURFACE'
1455	print, '.. and LW FLUX AT SURFACE'
1456
1457	localtime=localtime-history_interval_s/3700.
1458
1459	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_swdownz.nc', charvar='SWDOWNZ', invar=swdownz_les_tmp
1460	getcdf, file=file1, charvar='fluxsurf_sw', invar=swdownz_gcm
1461	;swdownz_les=make_array(nttot)
1462	;FOR l=0,nttot-1 DO BEGIN
1463	; swdownz_les(l)=TOTAL(TOTAL(swdownz_les_tmp(,,l),1),1)/float(nx) /float(ny)
1464	;ENDFOR
1465	;swdownz_les_tmp=0.
1466
1467	swdownz_les=temporary(swdownz)
1468	swdownz_les_Int=INTERPOL(swdownz_les,localtime,localtime_gcm)
1469	swdownz_gcm_Int=INTERPOL(swdownz_gcm,localtime_gcm,localtime,/QUADRATIC)
1470
1471	what_I_plot = swdownz_gcm-swdownz_les_Int
1472	labels=['TH1D SW flux - LES SW flux']
1473	title_user = TestCase+SubCase+LayerCase+' Difference in Solar Fluxes at surface'
1474	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_diffswdownz.ps'
1475	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1476	CS, SCALE=28
1477	GSET, XMIN=6, XMAX=18, YMIN=-20, YMAX=20, TITLE=title_user
1478	cols=INDGEN(1)+2
1479	GPLOT, X=localtime_gcm, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1480	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=2., YTITLE='Solar Flux (W.m-2)',NDECS=1
1481
1482	PSCLOSE, /NOVIEW
1483
1484	spawn, 'ps2png '+filename
1485
1486	what_I_plot = swdownz_les
1487	labels=['LES SW flux']
1488	title_user = TestCase+SubCase+LayerCase+' Solar Fluxes at surface'
1489	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_swdownz.ps'
1490	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1491	CS, SCALE=28
1492	GSET, XMIN=6, XMAX=18, YMIN=0, YMAX=500, TITLE=title_user
1493	cols=INDGEN(1)+2
1494	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1495	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=50., YTITLE='Solar Flux (W.m-2)',NDECS=1
1496
1497	oplot, localtime_gcm, swdownz_gcm, psym=1,thick=0.2
1498
1499	PSCLOSE, /NOVIEW
1500
1501	spawn, 'ps2png '+filename
1502
1503
1504	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_lwdownz.nc', charvar='LWDOWNZ', invar=lwdownz_les_tmp
1505	getcdf, file=file1, charvar='fluxsurf_lw', invar=lwdownz_gcm
1506	;lwdownz_les=make_array(nttot)
1507	;FOR l=0,nttot-1 DO BEGIN
1508	; lwdownz_les(l)=TOTAL(TOTAL(lwdownz_les_tmp(,,l),1),1)/float(nx) /float(ny)
1509	;ENDFOR
1510	;lwdownz_les_tmp=0.
1511
1512	lwdownz_les=temporary(lwdownz)
1513
1514	what_I_plot = lwdownz_les
1515	labels=['LES LW downward flux']
1516	title_user = TestCase+SubCase+LayerCase+' LW downward fluxes at surface'
1517	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_lwdownz.ps'
1518	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1519	CS, SCALE=28
1520	GSET, XMIN=6, XMAX=18, YMIN=-50, YMAX=50, TITLE=title_user
1521	cols=INDGEN(1)+2
1522	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1523	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=10., YTITLE='LW downward surface flux (W.m-2)',NDECS=1
1524
1525	oplot, localtime_gcm, lwdownz_gcm, psym=1,thick=0.2
1526
1527	PSCLOSE, /NOVIEW
1528
1529	spawn, 'ps2png '+filename
1530
1531
1532
1533
1534
1535	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_flxgrd.nc', charvar='FLXGRD', invar=fluxgrd_les_tmp
1536	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_flxgrd_gcmsoil.nc', charvar='FLXGRD', invar=fluxgrd_les_tmp
1537	getcdf, file=file1, charvar='fluxgrd', invar=fluxgrd_gcm
1538	;fluxgrd_les=make_array(nttot)
1539	;FOR l=0,nttot-1 DO BEGIN
1540	; fluxgrd_les(l)=TOTAL(TOTAL(fluxgrd_les_tmp(,,l),1),1)/float(n_elements(reform(fluxgrd_les_tmp(,0,0)))) /float(n_elements(reform(fluxgrd_les_tmp(0,,0))))
1541	;ENDFOR
1542	;fluxgrd_les_tmp=0.
1543
1544	fluxgrd_les=temporary(flxgrd)
1545
1546	what_I_plot = fluxgrd_les
1547	labels=['LES ground flux']
1548	title_user = TestCase+SubCase+LayerCase+' Ground flux at surface'
1549	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fluxgrd.ps'
1550	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1551	CS, SCALE=28
1552	GSET, XMIN=6, XMAX=18, YMIN=-70, YMAX=40, TITLE=title_user
1553	cols=INDGEN(1)+2
1554	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1555	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=10., YTITLE='Ground flux at surface (W.m-2)',NDECS=1
1556
1557	oplot, localtime_gcm, fluxgrd_gcm, psym=1,thick=0.2
1558
1559	PSCLOSE, /NOVIEW
1560
1561	spawn, 'ps2png '+filename
1562
1563
1564
1565
1566
1567	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_flxrad.nc', charvar='FLXRAD', invar=fluxrad_les_tmp
1568	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_flxrad_gcmsoil.nc', charvar='FLXRAD', invar=fluxrad_les_tmp
1569	getcdf, file=file1, charvar='fluxrad', invar=fluxrad_gcm
1570	;fluxrad_les=make_array(nttot)
1571	;FOR l=0,nttot-1 DO BEGIN
1572	; fluxrad_les(l)=TOTAL(TOTAL(fluxrad_les_tmp(,,l),1),1)/float(n_elements(reform(fluxrad_les_tmp(,0,0)))) /float(n_elements(reform(fluxrad_les_tmp(0,,0))))
1573	;ENDFOR
1574	;fluxrad_les_tmp=0.
1575
1576	fluxrad_les=temporary(flxrad)
1577
1578	what_I_plot = fluxrad_les
1579	labels=['LES radiation flux at surface']
1580	title_user = TestCase+SubCase+LayerCase+' Radiation flux at surface'
1581	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fluxrad.ps'
1582	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1583	CS, SCALE=28
1584	GSET, XMIN=6, XMAX=18, YMIN=-40, YMAX=80, TITLE=title_user
1585	cols=INDGEN(1)+2
1586	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1587	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=10., YTITLE='Radiation flux at surface (W.m-2)',NDECS=1
1588
1589	oplot, localtime_gcm, fluxrad_gcm, psym=1,thick=0.2
1590
1591	PSCLOSE, /NOVIEW
1592
1593	spawn, 'ps2png '+filename
1594
1595
1596	localtime=localtime+history_interval_s/3700.
1597
1598	endif
1599
1600	if (got_hfx eq 'true') then begin
1601
1602	print, '........ SENSIBLE HEAT FLUX'
1603
1604	localtime=localtime-history_interval_s/3700.
1605	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_hfx.nc', charvar='HFX', invar=hfx_les_tmp
1606	;getcdf, file=les_path+'/wrfout_d01_9999-01-01_hfx_gcmsoil.nc', charvar='HFX', invar=hfx_les_tmp
1607
1608	hfx_les=temporary(hfx)
1609
1610	getcdf, file=file1, charvar='hfx', invar=hfx_gcm
1611
1612	;hfx_les=make_array(nttot)
1613	;FOR l=0,nttot-1 DO BEGIN
1614	; hfx_les(l)=TOTAL(TOTAL(hfx_les_tmp(,,l),1),1)/float(n_elements(reform(hfx_les_tmp(,0,0)))) /float(n_elements(reform(hfx_les_tmp(0,,0))))
1615	;ENDFOR
1616	;hfx_les_tmp=0.
1617	;hfx_les_Int=INTERPOL(hfx_les,localtime,localtime_gcm)
1618	;hfx_gcm_Int=INTERPOL(hfx_gcm,localtime_gcm,localtime,/QUADRATIC)
1619
1620	what_I_plot = hfx_les
1621	labels=['LES sensible heat flux']
1622	title_user = TestCase+SubCase+LayerCase+' Sensible heat Fluxes at surface'
1623	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_hfx.ps'
1624	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1625	CS, SCALE=28
1626	GSET, XMIN=6, XMAX=18, YMIN=-50, YMAX=50, TITLE=title_user
1627	cols=INDGEN(1)+2
1628	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1629	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=10., YTITLE='Sensible heat Flux (W.m-2)',NDECS=1
1630
1631	oplot, localtime_gcm, hfx_gcm, psym=1,thick=0.2
1632
1633	PSCLOSE, /NOVIEW
1634
1635	spawn, 'ps2png '+filename
1636
1637
1638	if (1 eq 1) then begin
1639
1640	openw, lun, "input_hfx.def", /get_lun
1641	for l=0, nttot-1 do printf, lun, localtime(l),hfx_les(l),fluxrad_les(l)+fluxgrd_les(l)-hfx_les(l), format='((2x,F5.2)(4x,F8.2)(4x,F8.2))'
1642	FREE_LUN, lun
1643	close, lun
1644
1645	endif
1646
1647
1648
1649	what_I_plot = fluxrad_les+fluxgrd_les-hfx_les
1650	labels=['LES flux vdifc']
1651	title_user = TestCase+SubCase+LayerCase+' Flux vdifc at surface'
1652	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_vdifc.ps'
1653	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1654	CS, SCALE=28
1655	GSET, XMIN=6, XMAX=18, YMIN=-50, YMAX=50, TITLE=title_user
1656	cols=INDGEN(1)+2
1657	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1658	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=10., YTITLE='Heat Flux (W.m-2)',NDECS=1
1659
1660	oplot, localtime_gcm, fluxrad_gcm+fluxgrd_gcm-hfx_gcm, psym=1,thick=0.2
1661
1662	PSCLOSE, /NOVIEW
1663
1664	spawn, 'ps2png '+filename
1665
1666
1667
1668
1669
1670	;what_I_plot = [[fluxrad_les-swdownz_les-lwdownz_les],[fluxgrd_les],[hfx_les],[swdownz_les],[lwdownz_les],[fluxrad_les-fluxgrd_les-hfx_les]]
1671	what_I_plot = [[fluxrad_les],[swdownz_les],[lwdownz_les],[swdownz_les-lwdownz_les]]
1672	;labels=['lw up','grd','hfx','sw down','lw down','total']
1673	labels=['rad','sw down','lw down','sw +lw down']
1674	title_user = TestCase+SubCase+LayerCase+' Fluxes at surface'
1675	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_totflux.ps'
1676	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1677	CS, SCALE=28
1678	GSET, XMIN=6, XMAX=18, YMIN=-100, YMAX=100, TITLE=title_user
1679	cols=INDGEN(4)+2
1680	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1681	AXES, XSTEP = 1, XTITLE='local time (h)', YSTEP=10., YTITLE='Fluxes at surface (W.m-2)',NDECS=1
1682
1683	;oplot, localtime_gcm, fluxrad_gcm, psym=1,thick=0.2
1684
1685	PSCLOSE, /NOVIEW
1686
1687	spawn, 'ps2png '+filename
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698	localtime=localtime+history_interval_s/3700.
1699
1700
1701
1702
1703
1704
1705
1706
1707	endif
1708
1709	; ---------------------- * Vitesses verticales * -------------------------------
1710	; ------------ Verification de l'approx terrestre wmax = vmoy dans la couche instable
1711
1712	print, '........ CHECKING wmax = vmoy in unstable layer'
1713
1714	what_I_plot = uv_moy(*,lt_plotindex_les)
1715	labels=['LES uv_moy']
1716	title_user = TestCase+SubCase+LayerCase+' LES mean UV comp to max W in plume trac1'
1717	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_UV.ps'
1718	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1719	CS, SCALE=28
1720	GSET, XMIN=0, XMAX=10, YMIN=0, YMAX=10, TITLE=title_user
1721	cols=INDGEN(1)+2
1722	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1723	AXES, YSTEP = 1, YTITLE='Altitude (km)', XSTEP=1, XTITLE='Mean horizontal velocity inside domain (m/s)',NDECS=1
1724
1725	oplot, make_array(nz,value=wmax(lt_plotindex_les)), altitudes_LES/1000., psym=4
1726
1727	PSCLOSE, /NOVIEW
1728
1729	spawn, 'ps2png '+filename
1730
1731	; ------------ Profil de vitesse
1732
1733	print, '........ VERTICAL VELOCITY'
1734
1735	what_I_plot = make_array(nZmx)
1736	FOR k=0, nZmx-2 DO BEGIN
1737	what_I_plot(k) = 0.5*(zw2_lev(k,lt_plotindex_gcm) + zw2_lev(k+1,lt_plotindex_gcm))
1738	ENDFOR
1739	what_I_plot(nZmx-1) = 0.
1740
1741	smoothed_w_mean1_les = make_array(nz)
1742	smoothed_w_mean2_les = make_array(nz)
1743	smoothed_w_mean1_down_les = make_array(nz)
1744	FOR t=-ns,ns DO BEGIN
1745	smoothed_w_mean1_les = smoothed_w_mean1_les + REFORM(w_mean1(*,lt_plotindex_les+t))
1746	smoothed_w_mean2_les = smoothed_w_mean2_les + REFORM(w_mean2(*,lt_plotindex_les+t))
1747	smoothed_w_mean1_down_les = smoothed_w_mean1_down_les + REFORM(w_mean1_down(*,lt_plotindex_les+t))
1748	ENDFOR
1749
1750	smoothed_w_mean1_les = smoothed_w_mean1_les/nstot
1751	smoothed_w_mean2_les = smoothed_w_mean2_les/nstot
1752	smoothed_w_mean1_down_les = smoothed_w_mean1_down_les/nstot
1753
1754	ratio = make_array(nz)
1755	FOR k=0, nz-1 DO BEGIN
1756	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.
1757	ENDFOR
1758
1759	labels=['TH 1d w, lt='+string(lt_plot)]
1760	title_user = TestCase+SubCase+LayerCase+' Vertical velocity comparisons (inside thermal)'
1761	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Wprofile.ps'
1762	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1763	CS, SCALE=28
1764	GSET, XMIN=-6, XMAX=8, YMIN=0, YMAX=10, TITLE=title_user
1765	cols=INDGEN(1)+2
1766	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1767	AXES, YSTEP = 2, YTITLE='Altitude (km)', XSTEP=1, XTITLE='Vertical velocity inside thermal (m/s)',NDECS=1
1768
1769	oplot, smoothed_w_mean1_les, altitudes_LES/1000., psym=4
1770	oplot, smoothed_w_mean2_les, altitudes_LES/1000., psym=5
1771	oplot, smoothed_w_mean1_down_les, altitudes_LES/1000., psym=6
1772
1773	PSCLOSE, /NOVIEW
1774
1775	spawn, 'ps2png '+filename
1776
1777
1778	; * Static stability *
1779
1780	print, '........ STATIC STABILITY'
1781
1782	dteta_dz_gcm = deriv(altitudes_GCM,reform(teta_gcm(*,lt_plotindex_gcm)))
1783	dteta_dz_les = deriv(altitudes_LES,reform(teta_les(*,lt_plotindex_les)))
1784
1785	what_I_plot = dteta_dz_gcm
1786	labels=['TH static stability 1d']
1787	title_user = TestCase+SubCase+LayerCase+' Static stability comparison'
1788	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_dTetadz.ps'
1789	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1790	CS, SCALE=28
1791	GSET, XMIN=-0.002, XMAX=0.006, YMIN=0, YMAX=10, TITLE=title_user
1792	cols=INDGEN(1)+2
1793	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1794	AXES, XSTEP = 0.0005 , XTITLE='Static stability (K.m-1)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
1795
1796	oplot, dteta_dz_les, altitudes_LES/1000., psym=4
1797
1798	PSCLOSE, /NOVIEW
1799
1800	spawn, 'ps2png '+filename
1801
1802	print,' -----------------------------------------------------------------------------------------------------------------------'
1803	print,' * LES diagnostics of the PLUME * MUAHAHAHAHAHA'
1804	print,' V2 with E and D computed for a UDE plume'
1805	print,' -----------------------------------------------------------------------------------------------------------------------'
1806
1807	print, '........ EXTRACTING DATA'
1808
1809	; --- Reinterpolation of F
1810
1811	fm_therm_gcm_interlay = make_array(nZmx,nTmx)
1812
1813	FOR k=0, nZmx-2 DO BEGIN
1814	fm_therm_gcm_interlay(k,) = (fm_therm_gcm_lev(k,) + fm_therm_gcm_lev(k+1,*))/2.
1815	ENDFOR
1816	fm_therm_gcm_interlay(nZmx-1,*)=0.
1817
1818	; --- Calculation of gcm df/dz using entrainment and detrainments and NOT F
1819
1820	df_dz_gcm = deriv(altitudes_GCM,reform(fm_therm_gcm_interlay(*,lt_plotindex_gcm)))
1821	; --- Smoothing of the mass flux on a user-defined window
1822
1823	smoothed_fm_trac1_les = make_array(nz)
1824	smoothed_fm_trac2_les = make_array(nz)
1825	smoothed_downward_fm_trac1_les = make_array(nz)
1826	FOR t=-ns,ns DO BEGIN
1827	smoothed_fm_trac1_les = smoothed_fm_trac1_les + REFORM(fm_trac1_les(*,lt_plotindex_les+t))
1828	smoothed_fm_trac2_les = smoothed_fm_trac2_les + REFORM(fm_trac2_les(*,lt_plotindex_les+t))
1829	smoothed_downward_fm_trac1_les = smoothed_downward_fm_trac1_les + REFORM(downward_flux1(*,lt_plotindex_les+t))
1830	ENDFOR
1831
1832	smoothed_fm_trac1_les = smoothed_fm_trac1_les/nstot
1833	smoothed_fm_trac2_les = smoothed_fm_trac2_les/nstot
1834	smoothed_downward_fm_trac1_les = smoothed_downward_fm_trac1_les/nstot
1835
1836	; --- Calculation of the entrainement rate according to Rio(2010)
1837	; done in the heavy part at the begining (reeaaaally heavy)
1838
1839	; --- Smoothing of the entrainment on a ~30min window
1840
1841	; term 1
1842
1843
1844	smoothed_e_term1_trac1_les = make_array(nz)
1845	smoothed_e_term1_trac2_les = make_array(nz)
1846	FOR t=-ns,ns DO BEGIN
1847	smoothed_e_term1_trac1_les = smoothed_e_term1_trac1_les + REFORM(e_trac1_les(*,lt_plotindex_les+t))
1848	smoothed_e_term1_trac2_les = smoothed_e_term1_trac2_les + REFORM(e_trac2_les(*,lt_plotindex_les+t))
1849	ENDFOR
1850
1851	smoothed_e_term1_trac1_les = smoothed_e_term1_trac1_les/nstot
1852	smoothed_e_term1_trac2_les = smoothed_e_term1_trac2_les/nstot
1853
1854	smoothed_e_rate_term1_trac1_les = make_array(nz)
1855	smoothed_e_rate_trac2_les = smoothed_e_term1_trac2_les
1856
1857	; it already is an entrainment rate ! KIND OF : it is E/Mc, and Mc is not F !! NOW it is Mc/deltaTeta * dchi/dz
1858	FOR k=0, nz-1 DO BEGIN
1859	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.
1860	; 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.
1861	ENDFOR
1862
1863	; term 2 & 3
1864
1865
1866	smoothed_e_term2_trac1_les = make_array(nz)
1867	smoothed_e_term3_trac1_les = make_array(nz)
1868	smoothed_e_term1_ude_trac1_les = make_array(nz)
1869	smoothed_e_term2_ude_trac1_les = make_array(nz)
1870	smoothed_e_term3_ude_trac1_les = make_array(nz)
1871
1872	FOR t=-ns,ns DO BEGIN
1873	smoothed_e_term2_trac1_les = smoothed_e_term2_trac1_les + REFORM(e1_term2(*,lt_plotindex_les+t))
1874	smoothed_e_term2_trac1_les = smoothed_e_term2_trac1_les + REFORM(e1_term2(*,lt_plotindex_les+t))
1875	smoothed_e_term1_ude_trac1_les = smoothed_e_term1_ude_trac1_les + REFORM(e1_term1_ude(*,lt_plotindex_les+t))
1876	smoothed_e_term2_ude_trac1_les = smoothed_e_term2_ude_trac1_les + REFORM(e1_term2_ude(*,lt_plotindex_les+t))
1877	smoothed_e_term3_ude_trac1_les = smoothed_e_term3_ude_trac1_les + REFORM(e1_term3_ude(*,lt_plotindex_les+t))
1878	ENDFOR
1879
1880	smoothed_e_term2_trac1_les = smoothed_e_term2_trac1_les/nstot
1881	smoothed_e_term3_trac1_les = smoothed_e_term3_trac1_les/nstot
1882	smoothed_e_term1_ude_trac1_les = smoothed_e_term1_ude_trac1_les/nstot
1883	smoothed_e_term2_ude_trac1_les = smoothed_e_term2_ude_trac1_les/nstot
1884	smoothed_e_term3_ude_trac1_les = smoothed_e_term3_ude_trac1_les/nstot
1885
1886	smoothed_e_rate_term2_trac1_les = make_array(nz)
1887	smoothed_e_rate_term3_trac1_les = make_array(nz)
1888	smoothed_e_rate_term1_ude_trac1_les = make_array(nz)
1889	smoothed_e_rate_term2_ude_trac1_les = make_array(nz)
1890	smoothed_e_rate_term3_ude_trac1_les = make_array(nz)
1891
1892	FOR k=0, nz-1 DO BEGIN
1893	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.
1894	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.
1895	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.
1896	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.
1897	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.
1898
1899	ENDFOR
1900
1901	; --- Summing...
1902
1903	smoothed_e_rate_trac1_les = smoothed_e_rate_term1_trac1_les + smoothed_e_rate_term2_trac1_les + smoothed_e_rate_term3_trac1_les
1904	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
1905
1906	;print, 'ommiting term3'
1907	;smoothed_e_rate_trac1_les = smoothed_e_rate_term1_trac1_les + smoothed_e_rate_term2_trac1_les
1908
1909	; --- Smoothing of the detrainment rate
1910
1911	smoothed_d_term1_trac1_les = make_array(nz)
1912	smoothed_d_term2_trac1_les = make_array(nz)
1913	smoothed_d_term3_trac1_les = make_array(nz)
1914	smoothed_d_term1_ude_trac1_les = make_array(nz)
1915	smoothed_d_term2_ude_trac1_les = make_array(nz)
1916	smoothed_d_term3_ude_trac1_les = make_array(nz)
1917
1918	FOR t=-ns,ns DO BEGIN
1919	smoothed_d_term1_trac1_les = smoothed_d_term1_trac1_les + REFORM(d1_term1(*,lt_plotindex_les+t))
1920	smoothed_d_term2_trac1_les = smoothed_d_term2_trac1_les + REFORM(d1_term2(*,lt_plotindex_les+t))
1921	smoothed_d_term3_trac1_les = smoothed_d_term3_trac1_les + REFORM(d1_term3(*,lt_plotindex_les+t))
1922	smoothed_d_term1_ude_trac1_les = smoothed_d_term1_ude_trac1_les + REFORM(d1_term1_ude(*,lt_plotindex_les+t))
1923	smoothed_d_term2_ude_trac1_les = smoothed_d_term2_ude_trac1_les + REFORM(d1_term2_ude(*,lt_plotindex_les+t))
1924	smoothed_d_term3_ude_trac1_les = smoothed_d_term3_ude_trac1_les + REFORM(d1_term3_ude(*,lt_plotindex_les+t))
1925	ENDFOR
1926
1927	smoothed_d_term1_trac1_les = smoothed_d_term1_trac1_les/nstot
1928	smoothed_d_term2_trac1_les = smoothed_d_term2_trac1_les/nstot
1929	smoothed_d_term3_trac1_les = smoothed_d_term3_trac1_les/nstot
1930	smoothed_d_term1_ude_trac1_les = smoothed_d_term1_ude_trac1_les/nstot
1931	smoothed_d_term2_ude_trac1_les = smoothed_d_term2_ude_trac1_les/nstot
1932	smoothed_d_term3_ude_trac1_les = smoothed_d_term3_ude_trac1_les/nstot
1933
1934	smoothed_d_rate_term1_trac1_les = make_array(nz)
1935	smoothed_d_rate_term2_trac1_les = make_array(nz)
1936	smoothed_d_rate_term3_trac1_les = make_array(nz)
1937	smoothed_d_rate_term1_ude_trac1_les = make_array(nz)
1938	smoothed_d_rate_term2_ude_trac1_les = make_array(nz)
1939	smoothed_d_rate_term3_ude_trac1_les = make_array(nz)
1940
1941	FOR k=0, nz-1 DO BEGIN
1942	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.
1943	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.
1944	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.
1945	IF (smoothed_fm_trac1_les(k) ne 0.) THEN BEGIN
1946	smoothed_d_rate_term1_ude_trac1_les(k) = smoothed_d_term1_ude_trac1_les(k)/smoothed_fm_trac1_les(k)
1947	smoothed_d_rate_term2_ude_trac1_les(k) = smoothed_d_term2_ude_trac1_les(k)/smoothed_fm_trac1_les(k)
1948	smoothed_d_rate_term3_ude_trac1_les(k) = smoothed_d_term3_ude_trac1_les(k)/smoothed_fm_trac1_les(k)
1949	ENDIF ELSE BEGIN
1950	smoothed_d_rate_term1_ude_trac1_les(k)=0.
1951	smoothed_d_rate_term2_ude_trac1_les(k)=0.
1952	smoothed_d_rate_term3_ude_trac1_les(k)=0.
1953	ENDELSE
1954	ENDFOR
1955
1956	; --- Summing...
1957
1958	full_d_rate_ude = d1_term1_ude + d1_term2_ude + d1_term3_ude
1959
1960	smoothed_d_rate_trac1_les = smoothed_d_rate_term1_trac1_les+smoothed_d_rate_term2_trac1_les+smoothed_d_rate_term3_trac1_les
1961	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
1962	;print, 'ommiting term3'
1963	;smoothed_d_rate_trac1_les = smoothed_d_rate_term1_trac1_les+smoothed_d_rate_term2_trac1_les
1964
1965	; --- PLOTTING : BUOYANCY TERM
1966
1967	smoothed_buoyancy_trac1_les = make_array(nz)
1968	smoothed_buoyancy_ude_trac1_les = make_array(nz)
1969	smoothed_buoyancy_trac2_les = make_array(nz)
1970	smoothed_buoyancy_downdraft1_les_ude = make_array(nz)
1971
1972	FOR t=-ns,ns DO BEGIN
1973	smoothed_buoyancy_trac1_les = smoothed_buoyancy_trac1_les + REFORM(buoyancy1_les(*,lt_plotindex_les+t))
1974	smoothed_buoyancy_ude_trac1_les = smoothed_buoyancy_ude_trac1_les + REFORM(grav(tplume1moy(,lt_plotindex_les+t)/tenv1moy_ude(*,lt_plotindex_les+t)-1.))
1975	smoothed_buoyancy_trac2_les = smoothed_buoyancy_trac2_les + REFORM(buoyancy2_les(*,lt_plotindex_les+t))
1976	smoothed_buoyancy_downdraft1_les_ude = smoothed_buoyancy_downdraft1_les_ude + REFORM(grav(tdown1moy(,lt_plotindex_les+t)/tenv1moy_ude(*,lt_plotindex_les+t)-1.))
1977	ENDFOR
1978
1979	smoothed_buoyancy_trac1_les = smoothed_buoyancy_trac1_les/nstot
1980	smoothed_buoyancy_ude_trac1_les = smoothed_buoyancy_ude_trac1_les/nstot
1981	smoothed_buoyancy_trac2_les = smoothed_buoyancy_trac2_les/nstot
1982	smoothed_buoyancy_downdraft1_les_ude = smoothed_buoyancy_downdraft1_les_ude/nstot
1983
1984	print, '........ BUOYANCY'
1985
1986	what_I_plot = [[buoyancy_gcm(,lt_plotindex_gcm)],[buoyancy_est_gcm(,lt_plotindex_gcm)]]
1987	labels=['TH buoyancy term','TH estimated buoyancy in plume']
1988	title_user = TestCase+SubCase+LayerCase+' UDE plume buoyancy'
1989	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_B.ps'
1990	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
1991	CS, SCALE=28
1992	GSET, XMIN=-0.06, XMAX=0.06, YMIN=0, YMAX=10, TITLE=title_user
1993	cols=INDGEN(2)+2
1994	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
1995	AXES, XSTEP = 0.01 , XTITLE='N.m-1', YSTEP=1, YTITLE=' Altitude (km)',NDECS=3
1996
1997	;oplot, smoothed_buoyancy_trac1_les, altitudes_LES/1000., psym=4
1998	;oplot, smoothed_buoyancy_trac2_les, altitudes_LES/1000., psym=5
1999	print, smoothed_buoyancy_ude_trac1_les
2000	oplot, smoothed_buoyancy_ude_trac1_les, altitudes_LES/1000., psym=4
2001	oplot, smoothed_buoyancy_downdraft1_les_ude, altitudes_LES/1000., psym=7
2002
2003	PSCLOSE, /NOVIEW
2004
2005	spawn, 'ps2png '+filename
2006
2007
2008	; --- PLOTTING : MASS FLUX
2009
2010	print, '........ MASS FLUX'
2011
2012	f_gcm = fm_therm_gcm_interlay(*,lt_plotindex_gcm)
2013	what_I_plot = f_gcm
2014	labels=['TH mass flux']
2015	title_user = TestCase+SubCase+LayerCase+' mass flux comparison, average over '+taverage+' mn'
2016	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_f.ps'
2017	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2018	CS, SCALE=28
2019	GSET, XMIN=-0.008, XMAX=0.008, YMIN=0, YMAX=10, TITLE=title_user
2020	cols=INDGEN(1)+2
2021	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2022	AXES, XSTEP = 0.002 , XTITLE='Kg.m-2.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2023
2024	oplot, smoothed_fm_trac1_les, altitudes_LES/1000., psym=4
2025	oplot, smoothed_fm_trac2_les, altitudes_LES/1000., psym=5
2026	oplot, smoothed_downward_fm_trac1_les, altitudes_LES/1000., psym=6
2027
2028	PSCLOSE, /NOVIEW
2029
2030	spawn, 'ps2png '+filename
2031
2032
2033	; --- PLOTTING : MASS FLUX DERIVATIVE
2034
2035	print, '........ MASS FLUX DERIVATIVE'
2036
2037
2038	B_w2_trac1 = make_array(nz)
2039
2040	FOR k=0, nz-1 DO BEGIN
2041	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.
2042	; 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.
2043	ENDFOR
2044
2045	df_dz_les1 = deriv(altitudes_LES,reform(smoothed_fm_trac1_les))
2046	df_dz_les2 = deriv(altitudes_LES,reform(smoothed_fm_trac2_les))
2047	df_dz_param = make_array(nz)
2048	df_dz_param2 = make_array(nz)
2049
2050	;dlow=0.0013 ;baseline from continuity equation
2051	dlow=0.0003 ; svn baseline
2052	;dlow=0.0005
2053	dcoeff=-0.3
2054	;dcoeff=-0.4
2055
2056	;aaa1=2.5
2057	;bbb1=0.0015
2058	;Ae1=0.045
2059	;Be1=0.6
2060
2061	aaa1=1.60226
2062	bbb1=0.0006
2063	Ae1=0.0454
2064	Be1=0.57
2065
2066	FOR k=0, nz-1 DO BEGIN
2067	IF (2.5*B_w2_trac1(k) GE 0.) THEN BEGIN
2068	IF ((aaa1*B_w2_trac1(k)-bbb1) GE 0.) THEN BEGIN
2069	; df_dz_param(k)=smoothed_fm_trac1_les(k)(Ae1(aaa1B_w2_trac1(k)-bbb1)^(Be1) - 0.06(aaa1*B_w2_trac1(k))^(0.75))
2070	df_dz_param2(k)=smoothed_fm_trac1_les(k)(Ae1(aaa1B_w2_trac1(k)-bbb1)^(Be1) - MAX([(dcoeffB_w2_trac1(k) + dlow),0.]))
2071	ENDIF ELSE BEGIN
2072	; df_dz_param(k)=smoothed_fm_trac1_les(k)(-0.06(aaa1*B_w2_trac1(k))^(0.75))
2073	df_dz_param2(k)=smoothed_fm_trac1_les(k)(-MAX([(dcoeffB_w2_trac1(k) + dlow),0.]))
2074	ENDELSE
2075	ENDIF ELSE BEGIN
2076	; df_dz_param(k)=smoothed_fm_trac1_les(k)(-0.06(-aaa1*B_w2_trac1(k))^(0.75))
2077	df_dz_param2(k)=smoothed_fm_trac1_les(k)(-MAX([(dcoeffB_w2_trac1(k) + dlow),0.]))
2078	ENDELSE
2079	ENDFOR
2080
2081	what_I_plot = df_dz_gcm
2082	labels=['TH mass flux vertical derivative']
2083	title_user = TestCase+SubCase+LayerCase+' mass flux vertical derivative comparison, average over '+taverage+' mn'
2084	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_dfdz.ps'
2085	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2086	CS, SCALE=28
2087	GSET, XMIN=-0.00002, XMAX=0.00002, YMIN=0, YMAX=10, TITLE=title_user
2088	cols=INDGEN(1)+2
2089	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2090	AXES, XSTEP = 0.000005 , XTITLE='Kg.m-3.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=6
2091	oplot, df_dz_les1, altitudes_LES/1000., psym=4
2092	;oplot, df_dz_les2, altitudes_LES/1000., psym=5
2093	oplot, df_dz_param, altitudes_LES/1000., psym=6, color=5
2094	oplot, df_dz_param2, altitudes_LES/1000., psym=6, color=6
2095	;print, "fm*(e-d)"
2096	;print, smoothed_fm_trac1_les(0.045(2.5B_w2_trac1-0.0015)^(0.6) - 0.06(-2.5*B_w2_trac1)^(0.75))
2097	;print, smoothed_fm_trac1_les
2098	print, B_w2_trac1
2099
2100	PSCLOSE, /NOVIEW
2101
2102	spawn, 'ps2png '+filename
2103
2104	; --- PLOTTING : ENTRAINMENT RATE e = E/f
2105
2106	print, '........ ENTRAINMENT RATE'
2107
2108	e_gcm = make_array(nZmx)
2109
2110	FOR k=0, nZmx-1 DO BEGIN
2111	IF (fm_therm_gcm_interlay(k,lt_plotindex_gcm) ne 0.) THEN BEGIN
2112	e_gcm(k) = zdz_entr_therm_gcm(k,lt_plotindex_gcm)/(approx_zdz_gcm(k)*fm_therm_gcm_interlay(k,lt_plotindex_gcm))
2113	ENDIF ELSE BEGIN
2114	e_gcm(k) = 0.
2115	ENDELSE
2116	ENDFOR
2117
2118
2119	what_I_plot = e_gcm
2120	labels=['TH entrainment rate']
2121	title_user = TestCase+SubCase+LayerCase+' UDE entrainment rate comparison, average over '+taverage+' mn'
2122	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e.ps'
2123	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2124	CS, SCALE=28
2125	GSET, XMIN=-0.003, XMAX=0.003, YMIN=0, YMAX=10, TITLE=title_user
2126	cols=INDGEN(1)+2
2127	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2128	AXES, XSTEP = 0.0006 , XTITLE='entrainment rate m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2129
2130	oplot, smoothed_e_rate_ude_trac1_les, altitudes_LES/1000., psym=4
2131	;oplot, smoothed_e_rate_trac1_les, altitudes_LES/1000., psym=4
2132	;oplot, smoothed_e_rate_trac2_les, altitudes_LES/1000., psym=5
2133
2134	PSCLOSE, /NOVIEW
2135
2136	spawn, 'ps2png '+filename
2137
2138	print, '........ EXTENDED ENTRAINMENT RATE'
2139
2140	;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]]
2141	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]]
2142	labels=['LES base entrainment rate','LES term2 e rate','LES term3 e rate','LES total e rate']
2143	title_user = TestCase+SubCase+LayerCase+' UDE entrainment rate comparison, average over '+taverage+' mn'
2144	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e_terms.ps'
2145	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2146	CS, SCALE=28
2147	GSET, XMIN=-0.01, XMAX=0.01, YMIN=0, YMAX=10, TITLE=title_user
2148	cols=INDGEN(4)+2
2149	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2150	AXES, XSTEP = 0.005 , XTITLE='entrainment rate m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=3
2151
2152	PSCLOSE, /NOVIEW
2153
2154	spawn, 'ps2png '+filename
2155
2156
2157	; --- PLOTTING : EXTENDED DETRAINMENT RATE
2158
2159	print, '........ EXTENDED DETRAINMENT RATE'
2160
2161	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]]
2162	labels=['LES term 1 detrainment rate','LES term2 d rate','LES term3 d rate','LES Total d rate']
2163	title_user = TestCase+SubCase+LayerCase+' UDE detrainment rate comparison, average over '+taverage+' mn'
2164	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d_terms.ps'
2165	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2166	CS, SCALE=28
2167	GSET, XMIN=-0.01, XMAX=0.01, YMIN=0, YMAX=10, TITLE=title_user
2168	cols=INDGEN(4)+2
2169	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2170	AXES, XSTEP = 0.002 , XTITLE='detrainment rate m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2171
2172	PSCLOSE, /NOVIEW
2173
2174	spawn, 'ps2png '+filename
2175
2176	; --- PLOTTING : DETRAINMENT RATE d = D/f
2177
2178	print, '........ DETRAINMENT RATE'
2179
2180	;smoothed_d_rate_trac2_les = make_array(nz)
2181	;
2182	;FOR k=0, nz-1 DO BEGIN
2183	; 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.
2184	; 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.
2185	;ENDFOR
2186	;
2187	d_gcm = make_array(nZmx)
2188	FOR k=0, nZmx-1 DO BEGIN
2189	IF (fm_therm_gcm_interlay(k,lt_plotindex_gcm) ne 0.) THEN BEGIN
2190	d_gcm(k) = zdz_detr_therm_gcm(k,lt_plotindex_gcm)/(approx_zdz_gcm(k)*fm_therm_gcm_interlay(k,lt_plotindex_gcm))
2191	ENDIF ELSE BEGIN
2192	d_gcm(k) = 0.
2193	ENDELSE
2194	ENDFOR
2195
2196	what_I_plot = d_gcm
2197	labels=['TH detrainment rate']
2198	title_user = TestCase+SubCase+LayerCase+' UDE detrainment rate comparison, average over '+taverage+' mn'
2199	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d.ps'
2200	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2201	CS, SCALE=28
2202	GSET, XMIN=0., XMAX=0.03, YMIN=0, YMAX=10, TITLE=title_user
2203	cols=INDGEN(1)+2
2204	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2205	AXES, XSTEP = 0.005 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=2
2206
2207	oplot, smoothed_d_rate_ude_trac1_les, altitudes_LES/1000., psym=4
2208	;oplot, smoothed_d_rate_trac1_les, altitudes_LES/1000., psym=4
2209	;oplot, smoothed_d_rate_trac2_les, altitudes_LES/1000., psym=5
2210
2211	PSCLOSE, /NOVIEW
2212	spawn, 'ps2png '+filename
2213
2214	; --- PLOTTING : FRACTION COVERAGE
2215
2216	print, '........ EXTENDED ALPHA'
2217
2218	what_I_plot = alpha_interlay_gcm
2219	labels=['TH alpha']
2220	title_user = TestCase+SubCase+LayerCase+' fraction coverage comparison, average over '+taverage+' mn'
2221	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_alpha.ps'
2222	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2223	CS, SCALE=28
2224	GSET, XMIN=0., XMAX=1., YMIN=0, YMAX=10, TITLE=title_user
2225	cols=INDGEN(1)+2
2226	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2227	AXES, XSTEP = 0.1 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
2228
2229	oplot, smoothed_alpha1_les, altitudes_LES/1000., psym=4
2230	oplot, smoothed_alpha2_les, altitudes_LES/1000., psym=5
2231	oplot, smoothed_beta1_les, altitudes_LES/1000., psym=6
2232
2233	PSCLOSE, /NOVIEW
2234	spawn, 'ps2png '+filename
2235
2236	; --- PLOTTING : THEORETICAL ENTRAINMENT RATE FROM LES DATA
2237
2238	print, '........ PARAMETRIZED RATES'
2239
2240	approx_zdz_les = make_array(nz)
2241
2242	approx_zdz_les(0)=altitudes_LES(1)
2243	FOR k=1, nz-2 DO BEGIN
2244	approx_zdz_les(k) = altitudes_LES(k+1) - altitudes_LES(k)
2245	ENDFOR
2246	approx_zdz_les(nz-1)=approx_zdz_les(nz-2)
2247
2248
2249	theoretical_e_trac1_les = make_array(nz)
2250	theoretical_e_trac2_les = make_array(nz)
2251
2252
2253	FOR k=0, nz-1 DO BEGIN
2254	theoretical_e_trac1_les(k) = MAX([0.,(betalpha/(1.+betalpha))((afactsmoothed_buoyancy_trac1_les(k)/((smoothed_w_mean1_les(k))^2.)) - fact_epsilon)])
2255	theoretical_e_trac2_les(k) = MAX([0.,(betalpha/(1.+betalpha))((afactsmoothed_buoyancy_trac2_les(k)/((smoothed_w_mean2_les(k))^2.)) - fact_epsilon)])
2256	ENDFOR
2257
2258
2259	what_I_plot = [[theoretical_e_trac1_les],[theoretical_e_trac2_les]]
2260	labels=['LES TH theo e rate trac1','LES TH theo e rate trac2']
2261	title_user = TestCase+SubCase+LayerCase+' comp. theor. entr. rate comparison, average over '+taverage+' mn'
2262	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e_theoretical.ps'
2263	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2264	CS, SCALE=28
2265	GSET, XMIN=-0.015, XMAX=0.03, YMIN=0, YMAX=10, TITLE=title_user
2266	cols=INDGEN(2)+2
2267	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2268	AXES, XSTEP = 0.003 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=3
2269
2270	oplot, smoothed_e_rate_trac1_les, altitudes_LES/1000., psym=4
2271	oplot, smoothed_e_rate_trac2_les, altitudes_LES/1000., psym=5
2272
2273	PSCLOSE, /NOVIEW
2274
2275	spawn, 'ps2png '+filename
2276
2277	; --- PLOTTING : THEORETICAL DETRAINMENT RATE FROM LES DATA
2278	; ZDZ STUFF REMOVED
2279
2280	print, '........ PARAMETRIZED DETRAINMENT'
2281
2282	theoretical_d_trac1_les = make_array(nz)
2283	theoretical_d_trac2_les = make_array(nz)
2284
2285	FOR k=0, nz-1 DO BEGIN
2286	theoretical_d_trac1_les(k) = MAX([detr_min,-afact(betalpha/(1.+betalpha))(smoothed_buoyancy_trac1_les(k)/((smoothed_w_mean1_les(k))^2.))])
2287	theoretical_d_trac2_les(k) = MAX([detr_min,-afact(betalpha/(1.+betalpha))(smoothed_buoyancy_trac2_les(k)/((smoothed_w_mean2_les(k))^2.))])
2288	ENDFOR
2289
2290	what_I_plot = [[theoretical_d_trac1_les],[theoretical_d_trac2_les]]
2291	labels=['LES TH theo d rate trac1','LES TH theo d rate trac2']
2292	title_user = TestCase+SubCase+LayerCase+' comp. theor. detr. rate comparison, average over '+taverage+' mn'
2293	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d_theoretical.ps'
2294	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2295	CS, SCALE=28
2296	GSET, XMIN=-0.1, XMAX=0.1, YMIN=0, YMAX=10, TITLE=title_user
2297	cols=INDGEN(2)+2
2298	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2299	AXES, XSTEP = 0.01 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=2
2300
2301	oplot, smoothed_d_rate_trac1_les, altitudes_LES/1000., psym=4
2302	;oplot, smoothed_d_rate_trac2_les, altitudes_LES/1000., psym=5
2303
2304	PSCLOSE, /NOVIEW
2305
2306	spawn, 'ps2png '+filename
2307
2308	; --- PLOTTING : THEORETICAL E-D RATE FROM LES DATA
2309
2310	print, '........ PARAMETRIZED MASS FLUX DERIVATIVE'
2311
2312	theoretical_dfdz_f_trac1_les = make_array(nz)
2313	theoretical_dfdz_f_trac2_les = make_array(nz)
2314
2315	theoretical_dfdz_f_trac1_les = theoretical_e_trac1_les - theoretical_d_trac1_les
2316	theoretical_dfdz_f_trac2_les = theoretical_e_trac2_les - theoretical_d_trac2_les
2317
2318	df_dz_f_les1 = make_array(nz)
2319	df_dz_f_les2 = make_array(nz)
2320
2321	FOR k=0, nz-1 DO BEGIN
2322	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.
2323	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.
2324	ENDFOR
2325
2326	what_I_plot = [[theoretical_dfdz_f_trac1_les],[theoretical_dfdz_f_trac2_les]]
2327	labels=['LES TH theo 1/f df/dz trac1','LES TH theo 1/f df/dz trac2']
2328	title_user = TestCase+SubCase+LayerCase+' comp. theor. entr. - detr. rate comparison, average over '+taverage+' mn'
2329	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_dfdzf_theoretical.ps'
2330	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2331	CS, SCALE=28
2332	GSET, XMIN=-0.02, XMAX=0.02, YMIN=0, YMAX=10, TITLE=title_user
2333	cols=INDGEN(2)+2
2334	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2335	AXES, XSTEP = 0.01 , XTITLE='entr - detr (rates) m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2336
2337	oplot, df_dz_f_les1, altitudes_LES/1000., psym=4
2338	oplot, df_dz_f_les2, altitudes_LES/1000., psym=5
2339
2340	PSCLOSE, /NOVIEW
2341
2342	spawn, 'ps2png '+filename
2343
2344	; --- PLOTTING : e versus B/w2
2345
2346	print, '........ EXTENDED TURBULENCE'
2347
2348	buoyancy1_les_ude = grav*(tplume1moy/tenv1moy_ude -1.)
2349	;Gamma_full = buoyancy1_les_ude + Gamma_1 + Gamma_2 + Gamma_3
2350	Gamma_full = buoyancy1_les + Gamma_1 + Gamma_2 + Gamma_3
2351	;Gamma_full = buoyancy1_les + Gamma_1
2352
2353	time_indices = [[lt_plotindex_les0],[lt_plotindex_les0a],[lt_plotindex_les],[lt_plotindex_lesa],[lt_plotindex_les2],[lt_plotindex_les2a],[lt_plotindex_les3]]
2354	ctime = [[lt_plot0],[lt_plot0a],[lt_plot],[lt_plota],[lt_plot2],[lt_plot2a],[lt_plot3]]
2355
2356	openw, lun, "fit_ab_simple_thermiques_"+TestCase+SubCase+les_special, /get_lun
2357	printf, lun, " a1 ", " b1 "," LT"
2358	FREE_LUN, lun
2359	close, lun
2360
2361	openw, lun, "fit_ab_double_thermiques_"+TestCase+SubCase+les_special, /get_lun
2362	printf, lun, " a1a ", " b1a ", " a1b ", " b1b "," LT"
2363	FREE_LUN, lun
2364	close, lun
2365
2366	FOR ttt=0, n_elements(time_indices)-1 DO BEGIN
2367
2368	smooth_t,input=Gamma_full,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_Gamma_full
2369	smooth_t,input=buoyancy1_les,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_buoyancy1_les
2370	smooth_t,input=Gamma_1,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_Gamma_1
2371	smooth_t,input=Gamma_2,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_Gamma_2
2372	smooth_t,input=Gamma_3,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_Gamma_3
2373	smooth_t,input=w_mean1,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_w_mean1
2374	smooth_t,input=full_d_rate_ude,nz=nz,ndt=6,t0=time_indices(ttt),output=sm_full_d_rate_ude
2375
2376	www = where(w_mean1(*,time_indices(ttt))^2 ne 0.)
2377	nw=n_elements(www)
2378
2379	Y=make_array(nw-4)
2380	X=make_array(nw-4)
2381	FOR zz=0, nw-5 DO BEGIN
2382	; Y(zz) = Gamma_full(www(zz+2),lt_plotindex_les)/w_mean1(www(zz+2),lt_plotindex_les)^2
2383	; X(zz) = buoyancy1_les(www(zz+2),lt_plotindex_les)/w_mean1(www(zz+2),lt_plotindex_les)^2
2384
2385	; Approche Rio et al (2010)
2386	Y(zz) = sm_Gamma_full(www(zz+2))/sm_w_mean1(www(zz+2))^2
2387	X(zz) = sm_buoyancy1_les(www(zz+2))/sm_w_mean1(www(zz+2))^2
2388
2389	; Approche Gregory et al (2001)
2390	; Y(zz) = sm_Gamma_full(www(zz+2))/(sm_full_d_rate_ude(www(zz+2))*sm_w_mean1(www(zz+2))^2)
2391	; X(zz) = sm_buoyancy1_les(www(zz+2))/(sm_full_d_rate_ude(www(zz+2))*sm_w_mean1(www(zz+2))^2)
2392	ENDFOR
2393
2394	A = [2.5,-0.0015]
2395
2396	err_gamma=make_array(nw-4,value=0.001)
2397	FOR zz=floor(nw/4.), nw-5 DO BEGIN
2398	err_gamma(zz)=0.1
2399	ENDFOR
2400
2401	coefs = lmfit(X,Y,A,/DOUBLE,function_name = 'myfunct', itmax = 500, measure_error = err_gamma)
2402
2403	B = [2.5,-0.0015]
2404
2405	err_gamma=make_array(nw-4,value=0.1)
2406	FOR zz=floor(nw/4.), floor(nw*3./4.-5) DO BEGIN
2407	err_gamma(zz)=0.001
2408	ENDFOR
2409
2410	coefs = lmfit(X,Y,B,/DOUBLE,function_name = 'myfunct', itmax = 500, measure_error = err_gamma)
2411
2412	C = [1.5,-0.0010]
2413
2414	err_gamma=make_array(nw-4,value=0.001)
2415	FOR zz=floor(nw*3./4.-5), nw-5 DO BEGIN
2416	err_gamma(zz)=0.1
2417	ENDFOR
2418
2419	coefs = lmfit(X,Y,C,/DOUBLE,function_name = 'myfunct', itmax = 1000, measure_error = err_gamma)
2420
2421	openw, lun, 'fit_ab_simple_thermiques_'+TestCase+SubCase+les_special, /append
2422	printf, lun, C[0],C[1],ctime(ttt), format='((2x,F6.3)(4x,F9.6)(4x,I0))'
2423	close, lun
2424
2425	openw, lun, 'fit_ab_double_thermiques_'+TestCase+SubCase+les_special, /append
2426	printf, lun, A[0],A[1],B[0],B[1],ctime(ttt), format='((2x,F6.3)(4x,F9.6)(4x,F6.3)(4x,F9.6)(4x,I0))'
2427	close, lun
2428
2429
2430	print, '~~~~~> LT: '+string(ctime(ttt))+' <~~~~~~~~~'
2431	print, 'suggested coefs for fit, a,b in alim layer:'
2432	print, A
2433	print, 'suggested coefs for fit, a,b above alim layer:'
2434	print, B
2435	print, 'suggested coefs for uniform fit, a,b:'
2436	print, C
2437
2438	;what_I_plot = [[Gamma_full(,lt_plotindex_les)],[buoyancy1_les_ude(,lt_plotindex_les)],[Gamma_1(,lt_plotindex_les)],[Gamma_2(,lt_plotindex_les)],[Gamma_3(*,lt_plotindex_les)]]
2439
2440	;what_I_plot = [[Gamma_full(,lt_plotindex_les)],[buoyancy1_les(,lt_plotindex_les)],[Gamma_1(,lt_plotindex_les)],[Gamma_2(,lt_plotindex_les)],[Gamma_3(*,lt_plotindex_les)]]
2441
2442	what_I_plot = [[sm_Gamma_full()],[sm_buoyancy1_les()],[sm_Gamma_1()],[sm_Gamma_2()],[sm_Gamma_3(*)]]
2443	labels=['Tot','B','G1','G2','G3']
2444	title_user = TestCase+SubCase+LayerCase+' UDE turbulence term, average over '+taverage+' mn'
2445	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Gamma'+string(ctime(ttt),format='(I0)')+'.ps'
2446	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2447	CS, SCALE=28
2448	GSET, XMIN=-0.12, XMAX=0.12, YMIN=0, YMAX=6, TITLE=title_user
2449	cols=INDGEN(5)+2
2450	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2451	AXES, XSTEP = 0.03 , XTITLE='Gamma term m.s-2', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
2452
2453	oplot, -0.0015w_mean1(,lt_plotindex_les)^2, altitudes_LES/1000., psym=5, thick=0.05
2454	;oplot, 2.5buoyancy1_les_ude(,lt_plotindex_les), altitudes_LES/1000., psym=5, thick=0.05
2455	;oplot, 2.5buoyancy1_les_ude(,lt_plotindex_les)-0.0015w_mean1(,lt_plotindex_les)^2, altitudes_LES/1000., psym=4, thick=0.3
2456	;oplot, 2.5buoyancy1_les(,lt_plotindex_les), altitudes_LES/1000., psym=5, thick=0.05
2457	;oplot, 2.buoyancy1_les(,lt_plotindex_les)-0.0012w_mean1(,lt_plotindex_les)^2, altitudes_LES/1000., psym=4, thick=0.3
2458	oplot, 2.5buoyancy1_les(,lt_plotindex_les)-0.0015w_mean1(,lt_plotindex_les)^2, altitudes_LES/1000., psym=4, thick=0.3
2459	;oplot, A[0]buoyancy1_les(,lt_plotindex_les)+A[1]w_mean1(,lt_plotindex_les)^2, altitudes_LES/1000., thick=0.3
2460	;oplot, B[0]buoyancy1_les(,lt_plotindex_les)+B[1]w_mean1(,lt_plotindex_les)^2, altitudes_LES/1000., thick=0.3
2461
2462	; Rio et al 2010 :
2463	oplot, A[0]sm_buoyancy1_les()+A[1]sm_w_mean1()^2, altitudes_LES/1000., thick=0.3
2464	oplot, B[0]sm_buoyancy1_les()+B[1]sm_w_mean1()^2, altitudes_LES/1000., thick=0.3
2465	oplot, C[0]sm_buoyancy1_les()+C[1]sm_w_mean1()^2, altitudes_LES/1000., thick=0.3, linestyle=2
2466
2467	; Gregory et al 2001 :
2468	;oplot, A[0]sm_buoyancy1_les()+A[1]sm_full_d_rate_ude()sm_w_mean1()^2, altitudes_LES/1000., thick=0.3
2469	;oplot, B[0]sm_buoyancy1_les()+B[1]sm_full_d_rate_ude()sm_w_mean1()^2, altitudes_LES/1000., thick=0.3
2470	;oplot, C[0]sm_buoyancy1_les()+C[1]sm_full_d_rate_ude()sm_w_mean1()^2, altitudes_LES/1000., thick=0.3, linestyle=2
2471
2472
2473	PSCLOSE, /NOVIEW
2474
2475	spawn, 'ps2png '+filename
2476
2477	ENDFOR
2478
2479	print, ' ~~ comparing fit approaches :'
2480
2481	a1_simple=[1.848,1.723,1.455,1.280,1.993,1.789,1.582,1.050]
2482	b1_simple=[-0.000842,-0.000511,-0.000251,-0.000329,-0.001473,-0.001081,-0.000129,-0.000453]
2483	a1a_double=[1.820,1.716,1.454,1.274,1.918,1.767,1.568,1.029]
2484	b1a_double=[0.000015,-0.000331,-0.000202,-0.000158, 0.001490,-0.000450, 0.000330, 0.000500]
2485	a1b_double=[3.145,2.067,1.222,1.792,4.506,3.071,2.581,0.679]
2486	b1b_double=[-0.001913,-0.000716,-0.000227,-0.000519,-0.005111,-0.001961,-0.000667,-0.000788]
2487
2488	a1_simple=[1.993,1.870,1.789,1.811,1.582,1.565,1.050,1.607,1.494,1.415,1.508,1.523,1.267,1.266,2.008,1.840,1.683,1.732,1.676,1.577,1.348,1.849,1.933,1.757,1.767,1.659,1.608,1.529,1.992,1.918,1.748,1.737,1.605,1.384,0.083,1.848,1.745,1.723,1.598,1.455,1.473,1.280]
2489
2490	b1_simple=[-0.001473,-0.001211,-0.001081,-0.000492,-0.000129,-0.000343,-0.000453,-0.001194,-0.000856,-0.000386,-0.000072,-0.000187,-0.000292,-0.000806,-0.001287,-0.000947,-0.000587,-0.000326,-0.000759,-0.000560, 0.000108,-0.001072,-0.001034,-0.000792,-0.000685,-0.000468,-0.000724,-0.000208,-0.000995,-0.000851,-0.000374,-0.000493,-0.000140,-0.000291,-0.000450,-0.000842,-0.000569,-0.000511,-0.000496,-0.000251,-0.000276,-0.000329]
2491
2492	print, 'simple approach'
2493	print, 'mean a1 and mean b1'
2494	print, MEAN(a1_simple),MEAN(b1_simple)
2495	print, 'sigma'
2496	print, STDDEV(a1_simple),STDDEV(b1_simple)
2497
2498	print, 'double approach'
2499	print, 'mean a1a,b1a,a1b,b1b'
2500	print, MEAN(a1a_double),MEAN(b1a_double),MEAN(a1b_double),MEAN(b1b_double)
2501	print, 'sigma'
2502	print, STDDEV(a1a_double),STDDEV(b1a_double),STDDEV(a1b_double),STDDEV(b1b_double)
2503
2504	;aa1=2.5 & bb1=0.0015
2505	aa1=MEAN(a1_simple) & bb1=abs(MEAN(b1_simple))
2506
2507
2508	print, '........ BUOYANCY AND VERTICAL VELOCITY ENTRAINMENT RATE DEPENDENCY'
2509
2510	B_w2_trac2 = make_array(nz)
2511
2512	dwdz_trac1 = deriv(altitudes_LES,smoothed_w_mean1_les)
2513	;dwdz_trac2 = deriv(altitudes_LES,smoothed_w_mean2_les)
2514	full_dwdz_trac1 = make_array(nz,nttot)
2515	full_dadz_trac1 = make_array(nz,nttot)
2516	FOR l=0, nttot -1 DO BEGIN
2517	full_dwdz_trac1(,l) = deriv(altitudes_LES,w_mean1(,l))
2518	full_dadz_trac1(,l) = deriv(altitudes_LES,alpha1out(,l))
2519	ENDFOR
2520	;alpha = 0.
2521
2522	;FOR zzz=0.,30 DO BEGIN
2523
2524	;alpha = zzz/10.
2525	;
2526	;FOR k=0, nz-1 DO BEGIN
2527	; 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.
2528	; 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.
2529	;ENDFOR
2530
2531	;print, smoothed_buoyancy_trac1_les()/(smoothed_w_mean1_les())^2
2532	;print, (1./smoothed_w_mean1_les())dwdz_trac1(*)
2533
2534	full_e1 = make_array(nz,nttot)
2535	full_bw2 = make_array(nz,nttot)
2536	FOR k=0, nz-1 DO BEGIN
2537	FOR l=0, nttot-1 DO BEGIN
2538	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.
2539	; if(fm_trac1_les(k,l) ne 0.) then full_e1(k,l)=(e_trac1_les(k,l)+e1_term2(k,l)+e1_term3(k,l))/fm_trac1_les(k,l) else full_e1(k,l)=0.
2540	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.
2541	; if(w_mean1(k,l) ne 0.) then full_bw2(k,l)=grav*(tplume1moy(k,l)/tenv1moy(k,l) -1.)/(w_mean1(k,l)^2) else full_bw2(k,l)=0.
2542
2543	; 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.
2544	ENDFOR
2545	ENDFOR
2546
2547	;~~~~~~~~~~~~~~~ fit
2548
2549	lt_plotindex_les3=fix(lt_plotindex_les3(0))
2550	lt_plotindex_les0=fix(lt_plotindex_les0(0))
2551
2552	offset_fits=0
2553	IF (TestCase eq 'Case_Z') THEN BEGIN
2554	offset_fits =-10
2555	ENDIF
2556
2557	nttotfit = lt_plotindex_les3+offset_fits - lt_plotindex_les0
2558
2559	Gamma_w2=make_array(nz,nttotfit,value=0.)
2560
2561	FOR ttt=0,nttotfit-1 DO BEGIN
2562	www = where(w_mean1(*,ttt+lt_plotindex_les0) ne 0.)
2563	if (www(0) eq -1 ) then print, 'AIE AIE AIE!'
2564	Gamma_w2(www,ttt)=aa1*buoyancy1_les(www,ttt+lt_plotindex_les0)/(w_mean1(www,ttt+lt_plotindex_les0))^2 - bb1
2565	ENDFOR
2566
2567	D_out=make_array(2,nttotfit)
2568	;D_out=make_array(3,nttotfit)
2569
2570	print, 'begining fits :'
2571	print, lt_plotindex_les0, nttotfit, lt_plotindex_les3+offset_fits
2572
2573	FOR ttt=0,nttotfit-1 DO BEGIN
2574	www = where(Gamma_w2(*,ttt) gt 0.)
2575	eee = where(full_e1(*,ttt+lt_plotindex_les0) gt 0.)
2576	nw = n_elements(www)
2577	nee = n_elements(eee)
2578	if (nw gt nee) then begin
2579	Y = make_array(nee)
2580	X = make_array(nee)
2581	smooth_t,input=full_e1(eee,*),nz=nee,ndt=6,t0=ttt+lt_plotindex_les0,output=Y
2582	Y = reform(full_e1(eee,ttt+lt_plotindex_les0))
2583	X = reform(Gamma_w2(eee,ttt))
2584	endif else begin
2585	Y = make_array(nw)
2586	X = make_array(nw)
2587	smooth_t,input=full_e1(www,*),nz=nw,ndt=6,t0=ttt+lt_plotindex_les0,output=Y
2588	Y = reform(full_e1(www,ttt+lt_plotindex_les0))
2589	X = reform(Gamma_w2(www,ttt))
2590	endelse
2591	D = [0.08,0.6]
2592	;D = [0.08,0.6,0.1]
2593	coefs = lmfit(X,Y,D,/DOUBLE,function_name = 'myfunct2', itmax = 500)
2594	;coefs = lmfit(X,Y,D,/DOUBLE,function_name = 'myfunct3', itmax = 500)
2595	D_out(*,ttt)=D
2596	ENDFOR
2597
2598	openw, lun, "fit_power_epsilon_AB_thermiques_"+TestCase+SubCase+les_special, /get_lun
2599	printf, lun, " Ae "," Be"
2600	close, lun
2601
2602	openw, lun, "fit_power_epsilon_AB_thermiques_"+TestCase+SubCase+les_special, /append
2603	for l=0, nttotfit-1 do printf, lun, D_out(0,l),D_out(1,l), format='((2x,F8.6)(4x,F8.6))'
2604	FREE_LUN, lun
2605	close, lun
2606
2607	print, 'fits complete, output in fit_power_epsilon_AB_thermiques_'+TestCase+SubCase+les_special
2608
2609	print, 'power approach'
2610	print, 'mean Ae and mean Be'
2611	print, MEAN(D_out(0,)),MEAN(D_out(1,))
2612	print, 'sigma'
2613	print, STDDEV(D_out(0,)),STDDEV(D_out(1,))
2614
2615	;what_I_plot = smoothed_e_rate_ude_trac1_les
2616	what_I_plot = smoothed_e_rate_trac1_les
2617	labels=['e_rate trac1']
2618	title_user = TestCase+SubCase+LayerCase+' LES UDE entrainment rate dep with B/w2, average over '+taverage+' mn,'
2619	;filename = TestCase+SubCase+LayerCase+string(alpha,format='(F3.1)')+'Gcm_Les_Comp_e_Bw2.ps'
2620	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_e_Bw2.ps'
2621	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2622	CS, SCALE=28
2623	GSET, XMIN=0., XMAX=0.2, YMIN=0., YMAX=0.1, TITLE=title_user
2624	cols=INDGEN(1)+2
2625	GPLOT, Y=what_I_plot, X=B_w2_trac1, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30, SYM=5, /NOLINES
2626	AXES, XSTEP = 0.05 , XTITLE='B/w2', YSTEP=0.05, YTITLE='Entrainement rate m-1',NDECS=4
2627
2628	;oplot, smoothed_e_rate_trac2_les, B_w2_trac2, psym=5
2629	FOR l=lt_plotindex_les0, lt_plotindex_les3-1 DO BEGIN
2630	oplot, full_bw2(,l),full_e1(,l),thick=0.05,psym=1
2631	ENDFOR
2632	;mean_full_e1 = make_array(nz) & mean_full_bw2 = make_array(nz)
2633	;FOR k=0, nz-1 DO BEGIN
2634	; mean_full_e1(k) = MEAN(reform(full_e1(k,*)))
2635	; mean_full_bw2(k) = MEAN(reform(full_bw2(k,*)))
2636	;ENDFOR
2637	;oplot, mean_full_e1, mean_full_bw2, thick=0.3, psym = 2,color=5
2638	;oplot, theoretical_e_trac1_les, B_w2_trac1,psym=2,thick=0.8,color=7
2639	oplot,B_w2_trac1,(B_w2_trac1)/2.2222 + 0.0005,thick=0.3,color=7
2640	;oplot, 0.0118*(B_w2_trac1/0.043)^(1./1.65),B_w2_trac1,thick=0.3,color=7
2641	FOR l=lt_plotindex_les0, lt_plotindex_les3+offset_fits-1 DO BEGIN
2642	;oplot, full_bw2(,l),0.012(full_bw2(*,l)/0.048)^(1./1.60),thick=0.1,color=7,psym=1
2643
2644	;oplot, full_bw2(,l),0.04(2.5full_bw2(,l))^(0.5)-0.0015,thick=0.1,color=7,psym=1
2645	;oplot, full_bw2(,l),0.045(aa1full_bw2(,l)-bb1)^(0.6),thick=0.1,color=7,psym=1 ;entrainment formulation,baseline for a=2.5 & b=0.0015
2646
2647	;oplot, full_bw2(,l),0.06(aa1full_bw2(,l)-bb1)^(0.6),thick=0.1,color=7,psym=1
2648	oplot, full_bw2(,l),0.06(aa1full_bw2(,l))^(0.75),thick=0.1,color=6,psym=1 ;detrainment formulation
2649
2650	oplot, full_bw2(,l),MEAN(D_out(0,))(aa1full_bw2(,l)-bb1)^(MEAN(D_out(1,))),thick=0.1,color=7,psym=1
2651
2652
2653	ENDFOR
2654	beta1=0.15
2655
2656	;FOR l=0, nttot-1 DO BEGIN
2657	;;oplot, full_bw2(,l),beta1(2.5full_bw2(,l) - 0.0015)/(1.+beta1(1.-w_mean1_env(,l)/w_mean1(*,l))),thick=0.1,color=5,psym=1
2658	;oplot, full_bw2(,l),beta1(2.5full_bw2(,l) - 0.0015)/(1.+beta1),thick=0.1,color=5,psym=1 ; earth formulation
2659	;ENDFOR
2660
2661	;oplot, alog((B_w2_trac1 - 0.000942361)/0.0444855) - 3.85453, B_w2_trac1, thick=0.3,color=7
2662
2663	;print, alog((B_w2_trac1)/0.0444855) - 3.85453
2664
2665	PSCLOSE, /NOVIEW
2666
2667	spawn, 'ps2png '+filename
2668
2669	;ENDFOR
2670
2671	what_I_plot = full_bw2(*,lt_plotindex_les)
2672	labels=['B/w2']
2673	title_user = TestCase+SubCase+LayerCase+' LES UDE B/w2, average over '+taverage+' mn,'
2674	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Bw2.ps'
2675	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2676	CS, SCALE=28
2677	GSET, XMIN=-0.01, XMAX=0.01, YMIN=0., YMAX=6., TITLE=title_user
2678	cols=INDGEN(1)+2
2679	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2680	AXES, XSTEP = 0.002 , XTITLE='B/w2 term in LES UDE', YSTEP=0.5, YTITLE='Altitude (km)',NDECS=4
2681
2682	;FOR l=0, nttot-1 DO BEGIN
2683	; oplot, full_bw2(*,l),altitudes_LES/1000.,thick=0.05,psym=1
2684	;ENDFOR
2685
2686	PSCLOSE, /NOVIEW
2687
2688	spawn, 'ps2png '+filename
2689
2690
2691	print, '........ BUOYANCY AND VERTICAL VELOCITY DETRAINMENT RATE DEPENDENCY'
2692
2693	full_d1 = make_array(nz,nttot)
2694	full_dSiebesma = make_array(nz,nttot)
2695
2696
2697	FOR k=0, nz-1 DO BEGIN
2698	FOR l=0, nttot-1 DO BEGIN
2699	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.
2700	; if(fm_trac1_les(k,l) ne 0.) then full_d1(k,l)=(d1_term1(k,l)+d1_term2(k,l)+d1_term3(k,l))/fm_trac1_les(k,l) else full_d1(k,l)=-0.
2701	if(w_mean1(k,l) ne 0.) then full_dSiebesma(k,l)=0.750.5buoyancy1_les(k,l)/(w_mean1(k,l)^2) -1.5*full_dwdz_trac1(k,l)/w_mean1(k,l) - full_dadz_trac1(k,l)/alpha1out(k,l) else full_dSiebesma(k,l)=-0.
2702	ENDFOR
2703	ENDFOR
2704
2705	;what_I_plot = smoothed_d_rate_ude_trac1_les
2706	what_I_plot = smoothed_d_rate_trac1_les
2707	labels=['d_rate trac1']
2708	title_user = TestCase+SubCase+LayerCase+' LES UDE detrainment rate dep with B/w2, average over '+taverage+' mn'
2709	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d_Bw2.ps'
2710	;filename = TestCase+SubCase+LayerCase+string(alpha,format='(F3.1)')+'Gcm_Les_Comp_d_Bw2.ps'
2711	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2712	CS, SCALE=28
2713	GSET, XMIN=-0.1, XMAX=0.1, YMIN=0., YMAX=0.1, TITLE=title_user
2714	cols=INDGEN(1)+2
2715	GPLOT, Y=what_I_plot, X=full_bw2(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30, SYM=5, /NOLINES
2716	AXES, XSTEP = 0.05 , XTITLE='B/w2', YSTEP=0.05, YTITLE='Detrainment rate',NDECS=4
2717
2718	FOR l=0, nttot-1 DO BEGIN
2719	oplot, full_bw2(,l),full_d1(,l),thick=0.05,psym=1
2720	ENDFOR
2721	;oplot, theoretical_d_trac1_les, full_bw2(*,lt_plotindex_les),psym=2,thick=0.8,color=7
2722	;oplot,B_w2_trac1/2.7 + 0.0002,B_w2_trac1,thick=0.3,color=7
2723	oplot,B_w2_trac1,B_w2_trac1/2.222 + 0.0002,thick=0.3,color=7
2724	FOR l=0, nttot-1 DO BEGIN
2725	oplot, full_bw2(,l),0.06(2.5full_bw2(,l))^(0.75),thick=0.1,color=7,psym=1 ;detrainment formulation
2726	;oplot, full_bw2(,l),0.045(2.5full_bw2(,l)-0.0015)^(0.6),thick=0.1,color=6,psym=1 ;entrainment formulation
2727	ENDFOR
2728	FOR l=0, nttot-1 DO BEGIN
2729	oplot, full_bw2(,l),0.06(-2.5full_bw2(,l))^(0.75),thick=0.1,color=7,psym=1 ;detrainment formulation
2730	ENDFOR
2731	;a1=3.
2732
2733	;FOR l=0, nttot-1 DO BEGIN
2734	;oplot, full_bw2(,l),a1(beta1/(1.+beta1))full_bw2(,l),thick=0.1,color=5,psym=1 ;earth formulation
2735	;ENDFOR
2736
2737	;oplot, 0.0118*(B_w2_trac1/0.043)^(1./1.65),B_w2_trac1,thick=0.3,color=6
2738	oplot, B_w2_trac1,0.0105*(B_w2_trac1/0.048)^(1./1.7),thick=0.3,color=6
2739
2740	PSCLOSE, /NOVIEW
2741
2742	spawn, 'ps2png '+filename
2743
2744	print, 'Detrainment : new approach '
2745
2746	full_d1_v2 = make_array(nz,nttot)
2747	df_dz_les_full = make_array(nz,nttot)
2748	da_dt_les_full = make_array(nz,nttot)
2749
2750	FOR k=0, nz-1 DO BEGIN
2751	da_dt_les_full(k,) = deriv(localtime,reform(alpha1out(k,)))/3700.
2752	ENDFOR
2753
2754	FOR l=0, nttot-1 DO BEGIN
2755	df_dz_les_full(,l) = deriv(altitudes_LES,reform(fm_trac1_les(,l)))
2756	FOR k=0, nz-1 DO BEGIN
2757	if(fm_trac1_les(k,l) ne 0.) then full_d1_v2(k,l)=full_e1(k,l) - df_dz_les_full(k,l)/fm_trac1_les(k,l) -(rho(k,l)/fm_trac1_les(k,l))*da_dt_les_full(k,l) else full_d1_v2(k,l)=0.
2758	; if(fm_trac1_les(k,l) ne 0.) then full_d1_v2(k,l)=full_e1(k,l) - df_dz_les_full(k,l)/fm_trac1_les(k,l) else full_d1_v2(k,l)=0.
2759	ENDFOR
2760	ENDFOR
2761
2762
2763	;; ~~~~~~~~~~~~~~~~~~~~ DETRAINMENT FIT
2764
2765	;lt_plotindex_les3=fix(lt_plotindex_les3(0))
2766	;lt_plotindex_les0=fix(lt_plotindex_les0(0))
2767
2768	offset_fits=0
2769	IF (TestCase eq 'Case_Z') THEN BEGIN
2770	offset_fits =-10
2771	ENDIF
2772
2773	nttotfit = lt_plotindex_les3+offset_fits - lt_plotindex_les0
2774
2775	B_w2_fits=make_array(nz,nttotfit,value=0.)
2776
2777	FOR ttt=0,nttotfit-1 DO BEGIN
2778	www = where(w_mean1(*,ttt+lt_plotindex_les0) ne 0.)
2779	if (www(0) eq -1 ) then print, 'AIE AIE AIE!'
2780	; B_w2_fits(www,ttt)=aa1*buoyancy1_les(www,ttt+lt_plotindex_les0)/(w_mean1(www,ttt+lt_plotindex_les0))^2
2781	B_w2_fits(www,ttt)=buoyancy1_les(www,ttt+lt_plotindex_les0)/(w_mean1(www,ttt+lt_plotindex_les0))^2
2782	ENDFOR
2783
2784	E_out=make_array(2,nttotfit)
2785	print, 'detrainment: begining fits :'
2786	print, lt_plotindex_les0, nttotfit, lt_plotindex_les3+offset_fits
2787
2788	FOR ttt=0,nttotfit-1 DO BEGIN
2789	eee = where((full_d1_v2(,ttt+lt_plotindex_les0) gt 0.) and (abs(B_w2_fits(,ttt)) gt 0.001))
2790	nee = n_elements(eee)
2791
2792	eee=eee(4:nee-1)
2793	nee = n_elements(eee)
2794
2795	Y = make_array(nee)
2796	X = make_array(nee)
2797
2798	smooth_t,input=full_d1_v2(eee,*),nz=nee,ndt=6,t0=ttt+lt_plotindex_les0,output=Y
2799
2800	;Y = reform(full_d1_v2(eee,ttt+lt_plotindex_les0))
2801	X = reform(B_w2_fits(eee,ttt))
2802
2803	E = [-0.38,0.0005]
2804	coefs = lmfit(X,Y,E,/DOUBLE,function_name = 'myfunct', itmax = 1000)
2805	E_out(*,ttt)=E
2806	ENDFOR
2807
2808	openw, lun, "fit_lin_delta_AB_thermiques_"+TestCase+SubCase+les_special, /get_lun
2809	printf, lun, " Ad "," Bd"
2810	close, lun
2811
2812	openw, lun, "fit_lin_delta_AB_thermiques_"+TestCase+SubCase+les_special, /append
2813	for l=0, nttotfit-1 do printf, lun, E_out(0,l),E_out(1,l), format='((2x,F9.6)(4x,F9.6))'
2814	FREE_LUN, lun
2815	close, lun
2816
2817	print, 'fits complete, output in fit_lin_delta_AB_thermiques_'+TestCase+SubCase+les_special
2818
2819	print, 'delta: lin approach'
2820	print, 'mean Ad and mean Bd'
2821	print, MEAN(E_out(0,)),MEAN(E_out(1,))
2822	print, 'sigma'
2823	print, STDDEV(E_out(0,)),STDDEV(E_out(1,))
2824
2825	;~~~~~~~~~~~
2826
2827
2828	;what_I_plot = smoothed_d_rate_ude_trac1_les
2829	what_I_plot = full_d1_v2(*,lt_plotindex_les)
2830	labels=['d_rate trac1']
2831	title_user = TestCase+SubCase+LayerCase+' LES UDE detrainment rate dep with B/w2, average over '+taverage+' mn'
2832	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_d2_Bw2.ps'
2833	;filename = TestCase+SubCase+LayerCase+string(alpha,format='(F3.1)')+'Gcm_Les_Comp_d_Bw2.ps'
2834	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2835	CS, SCALE=28
2836	GSET, XMIN=-0.1, XMAX=0.1, YMIN=0., YMAX=0.1, TITLE=title_user
2837	cols=INDGEN(1)+2
2838	GPLOT, Y=what_I_plot, X=full_bw2(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30, SYM=5, /NOLINES
2839	AXES, XSTEP = 0.01 , XTITLE='B/w2', YSTEP=0.01, YTITLE='Detrainment rate',NDECS=2
2840
2841	FOR l=lt_plotindex_les0, lt_plotindex_les3+offset_fits-1 DO BEGIN
2842	oplot, full_bw2(,l),full_d1_v2(,l),thick=0.05,psym=1
2843	ENDFOR
2844	;oplot, theoretical_d_trac1_les, full_bw2(*,lt_plotindex_les),psym=2,thick=0.8,color=7
2845	;oplot,B_w2_trac1/2.7 + 0.0002,B_w2_trac1,thick=0.3,color=7
2846	oplot,B_w2_trac1,B_w2_trac1/2.222 + 0.0002,thick=0.3,color=7
2847	FOR l=lt_plotindex_les0, lt_plotindex_les3+offset_fits-1 DO BEGIN
2848
2849	oplot, full_bw2(,l), MEAN(E_out(0,))full_bw2(,l)+MEAN(E_out(1,*)),thick=0.1,color=7,psym=1 ;new detrainment formulation
2850
2851	oplot, full_bw2(,l), -0.45full_bw2(*,l)+0.0005,thick=0.1,color=8,psym=1 ;new detrainment formulation
2852
2853	;oplot, full_bw2(,l),0.06(2.5full_bw2(,l))^(0.75),thick=0.1,color=7,psym=1 ;detrainment formulation (classical)
2854	;oplot, full_bw2(,l),0.045(2.5full_bw2(,l)-0.0015)^(0.6),thick=0.1,color=6,psym=1 ;entrainment formulation
2855	ENDFOR
2856	;FOR l=0, nttot-1 DO BEGIN
2857	;oplot, full_bw2(,l),0.06(-2.5full_bw2(,l))^(0.75),thick=0.1,color=7,psym=1 ;detrainment formulation
2858	;ENDFOR
2859	;a1=3.
2860
2861	;FOR l=0, nttot-1 DO BEGIN
2862	;oplot, full_bw2(,l),a1(beta1/(1.+beta1))full_bw2(,l),thick=0.1,color=5,psym=1 ;earth formulation
2863	;ENDFOR
2864
2865	;oplot, 0.0118*(B_w2_trac1/0.043)^(1./1.65),B_w2_trac1,thick=0.3,color=6
2866	oplot, B_w2_trac1,0.0105*(B_w2_trac1/0.048)^(1./1.7),thick=0.3,color=6
2867
2868	PSCLOSE, /NOVIEW
2869
2870	spawn, 'ps2png '+filename
2871
2872
2873
2874	;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ mass flux derivative fits
2875
2876
2877	;offset_fits=0
2878	;IF (TestCase eq 'Case_Z') THEN BEGIN
2879	;offset_fits =-10
2880	;ENDIF
2881
2882	;nttotfit = lt_plotindex_les3+offset_fits - lt_plotindex_les0
2883
2884	full_d1_fits=make_array(nz,nttotfit,value=0.)
2885
2886
2887	full_1_f_dfdz = make_array(nz,nttot)
2888
2889	FOR l=0, nttot-1 DO BEGIN
2890	FOR k=0,nz-1 DO BEGIN
2891	IF (fm_trac1_les(k,l) GT 0.) THEN full_1_f_dfdz(k,l)=df_dz_les_full(k,l)/fm_trac1_les(k,l) ELSE full_1_f_dfdz(k,l)=0.
2892	ENDFOR
2893	ENDFOR
2894
2895	e_from_fit = make_array(nz,nttotfit)
2896	FOR l=0, nttotfit-1 DO BEGIN
2897	FOR k=0,nz-1 DO BEGIN
2898	e_from_fit(k,l) = MEAN(D_out(0,))(MAX([0.,aa1full_bw2(k,l+lt_plotindex_les0)-bb1]))^(MEAN(D_out(1,)))
2899	ENDFOR
2900	ENDFOR
2901
2902
2903	FOR ttt=0,nttotfit-1 DO BEGIN
2904	www = where(e_from_fit(,ttt) gt full_1_f_dfdz(,ttt+lt_plotindex_les0) )
2905	if (www(0) eq -1 ) then print, 'AIE AIE AIE!'
2906	full_d1_fits(www,ttt)=e_from_fit(www,ttt) - full_1_f_dfdz(www,ttt+lt_plotindex_les0)
2907	ENDFOR
2908
2909	F_out=make_array(3,nttotfit)
2910	print, 'mass flux derivative: begining fits :'
2911	print, lt_plotindex_les0, nttotfit, lt_plotindex_les3+offset_fits
2912
2913	FOR ttt=0,nttotfit-1 DO BEGIN
2914	;FOR ttt=3,nttotfit-4 DO BEGIN
2915
2916	eee = where((full_d1_fits(,ttt) gt 0.) and (altitudes_LES() gt 500.) and (B_w2_fits(*,ttt) gt 0.))
2917	fff = where((full_d1_fits(,ttt) gt 0.) and (altitudes_LES() gt 500.) and (B_w2_fits(*,ttt) lt 0.))
2918	nee = n_elements(eee)
2919	nff = n_elements(fff)
2920
2921	;eee=eee(floor(nee/4.):nee-1)
2922	;nee = n_elements(eee)
2923
2924	Y1 = make_array(nee)
2925	X1 = make_array(nee)
2926
2927	Y2 = make_array(nff)
2928	X2 = make_array(nff)
2929
2930	;smooth_t,input=full_d1_v2(eee,*),nz=nee,ndt=6,t0=ttt+lt_plotindex_les0,output=Y
2931
2932	Y1 = reform(full_d1_fits(eee,ttt))
2933	X1 = reform(B_w2_fits(eee,ttt))
2934
2935
2936	Y2 = reform(full_d1_fits(fff,ttt))
2937	X2 = reform(B_w2_fits(fff,ttt))
2938
2939	;smooth_t,input=full_d1_fits(eee,*),nz=nee,ndt=4,t0=ttt,output=Y
2940	;smooth_t,input=B_w2_fits(eee,*),nz=nee,ndt=4,t0=ttt,output=X
2941
2942	;F = [-0.38,0.0001]
2943	F = -0.5
2944	coefs = lmfit(X2,Y2,F,/DOUBLE,function_name = 'myfunct4', itmax = 1000)
2945
2946	F_out(0,ttt)=MEAN(reform(Y1))
2947	;F_out(1:2,ttt)=F
2948	F_out(1,ttt)=F
2949
2950	ENDFOR
2951
2952	;openw, lun, "fit_lin_delta_on_f_AB_thermiques_"+TestCase+SubCase+les_special, /get_lun
2953	;printf, lun, " Ada "," Ad "," Bd"
2954	;close, lun
2955	;
2956	;openw, lun, "fit_lin_delta_on_f_AB_thermiques_"+TestCase+SubCase+les_special, /append
2957	;for l=0, nttotfit-1 do printf, lun, F_out(0,l),F_out(1,l),F_out(2,l), format='((2x,F9.6)(4x,F9.6)(4x,F9.6))'
2958	;FREE_LUN, lun
2959	;close, lun
2960
2961	openw, lun, "fit_lin_delta_on_f_AB_thermiques_"+TestCase+SubCase+les_special, /get_lun
2962	printf, lun, " Ad "," Bd"
2963	close, lun
2964
2965	openw, lun, "fit_lin_delta_on_f_AB_thermiques_"+TestCase+SubCase+les_special, /append
2966	for l=0, nttotfit-1 do printf, lun, F_out(0,l),F_out(1,l), format='((2x,F9.6)(4x,F9.6))'
2967	FREE_LUN, lun
2968	close, lun
2969
2970
2971	print, 'delta:lin approach'
2972	print, 'mean Ada and mean Ad, Bd'
2973	print, MEAN(F_out(0,)),MEAN(F_out(1,)) ;,MEAN(F_out(2,*))
2974	print, 'sigma'
2975	print, STDDEV(F_out(0,)),STDDEV(F_out(1,)) ;,STDDEV(F_out(2,*))
2976
2977	what_I_plot = full_1_f_dfdz(*,lt_plotindex_les)
2978	labels=['TH mass flux vertical derivative']
2979	title_user = TestCase+SubCase+LayerCase+' mass flux vertical derivative comparison, average over '+taverage+' mn'
2980	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_1_f_dfdz.ps'
2981	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
2982	CS, SCALE=28
2983	GSET, XMIN=-0.02, XMAX=0.02, YMIN=0, YMAX=10, TITLE=title_user
2984	cols=INDGEN(1)+2
2985	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
2986	AXES, XSTEP = 0.005 , XTITLE='m-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=6
2987
2988	df_dz_from_param = make_array(nz)
2989	;FOR k=0, nz-1 DO BEGIN
2990	; df_dz_from_param(k) = MEAN(D_out(0,))(MAX([0.,aa1full_bw2(k,lt_plotindex_les)-bb1]))^(MEAN(D_out(1,))) - MEAN(F_out(0,))full_bw2(k,lt_plotindex_les)+MEAN(F_out(1,*))
2991	;ENDFOR
2992	FOR k=0, nz-1 DO BEGIN
2993	IF (full_bw2(k,lt_plotindex_les) gt 0.) THEN BEGIN
2994	df_dz_from_param(k) = MEAN(D_out(0,))(MAX([0.,aa1full_bw2(k,lt_plotindex_les)-bb1]))^(MEAN(D_out(1,))) - MEAN(F_out(0,*))
2995	ENDIF ELSE BEGIN
2996	; df_dz_from_param(k) = MEAN(D_out(0,))(MAX([0.,aa1full_bw2(k,lt_plotindex_les)-bb1]))^(MEAN(D_out(1,))) - (MEAN(F_out(1,))full_bw2(k,lt_plotindex_les)+MEAN(F_out(2,*)))
2997	df_dz_from_param(k) = MEAN(D_out(0,))(MAX([0.,aa1full_bw2(k,lt_plotindex_les)-bb1]))^(MEAN(D_out(1,))) - (MEAN(F_out(1,))full_bw2(k,lt_plotindex_les))
2998	ENDELSE
2999	ENDFOR
3000
3001	oplot, df_dz_from_param, altitudes_LES/1000., psym=5
3002	;oplot, df_dz_les1/smoothed_fm_trac1_les, altitudes_LES/1000., psym=4
3003	;oplot, df_dz_les2, altitudes_LES/1000., psym=5
3004	;oplot, df_dz_param, altitudes_LES/1000., psym=6, color=5
3005	;oplot, df_dz_param2, altitudes_LES/1000., psym=6, color=6
3006	;print, "fm*(e-d)"
3007	;print, smoothed_fm_trac1_les(0.045(2.5B_w2_trac1-0.0015)^(0.6) - 0.06(-2.5*B_w2_trac1)^(0.75))
3008	;print, smoothed_fm_trac1_les
3009	;print, B_w2_trac1
3010
3011	PSCLOSE, /NOVIEW
3012
3013	spawn, 'ps2png '+filename
3014
3015
3016
3017
3018	print, '........ BUOYANCY AND VERTICAL VELOCITY VERTICAL VELOCITY DERIVATIVE DEPENDANCY'
3019
3020	full_dw2_w2 = make_array(nz,nttot)
3021	w2_mean1 = make_array(nz,nttot)
3022	dw2_dz = make_array(nz,nttot)
3023	w2_param = make_array(nz,nttot)
3024	FOR l=0, nttot-1 DO BEGIN
3025	w2_mean1(,l) = w_mean1(,l)^2.
3026	dw2_dz(,l) = deriv(altitudes_LES,w2_mean1(,l))
3027	ENDFOR
3028	FOR k=0, nz-1 DO BEGIN
3029	FOR l=0, nttot-1 DO BEGIN
3030	if(w2_mean1(k,l) ne 0.) then begin
3031	; if(full_bw2(k,l) ge 0.) then full_dw2_w2(k,l) = 0.5dw2_dz(k,l)/w2_mean1(k,l) + (1.-w_mean1_env_ude(k,l)/w_mean1(k,l))0.012(full_bw2(k,l)/0.048)^(1./1.6) else full_dw2_w2(k,l)=0.5dw2_dz(k,l)/w2_mean1(k,l)
3032	; if(full_bw2(k,l) ge 0.) then full_dw2_w2(k,l) = 0.5dw2_dz(k,l)/w2_mean1(k,l) + (1.-w_mean1_env(k,l)/w_mean1(k,l))0.009(full_bw2(k,l)/0.048)^(1./1.9) else full_dw2_w2(k,l)=0.5dw2_dz(k,l)/w2_mean1(k,l)
3033	full_dw2_w2(k,l) = 0.5dw2_dz(k,l)/w2_mean1(k,l); + (1.-w_mean1_env_ude(k,l)/w_mean1(k,l))full_e1(k,l)
3034	endif else begin
3035	full_dw2_w2(k,l)= 0.
3036	endelse
3037
3038	if(w2_mean1(k,l) ne 0.) then begin
3039	; if(full_bw2(k,l) ge 0.) then w2_param(k,l) = full_bw2(k,l)/6. - (0.014*(full_bw2(k,l)/0.05)^(1./1.35))/2. else w2_param(k,l) = full_bw2(k,l)/6. - full_d1(k,l)/2.
3040	w2_param(k,l) = 2.5*full_bw2(k,l) - 0.0015
3041	endif else begin
3042	w2_param(k,l)=0.
3043	endelse
3044	; if(w2_mean1(k,l) ne 0.) then w2_param(k,l) = full_bw2(k,l) else w2_param(k,l)=0.
3045	ENDFOR
3046	ENDFOR
3047
3048	what_I_plot = full_dw2_w2(*,lt_plotindex_les)
3049
3050	labels=['0.5*(dw2/dz)/w2']
3051	title_user = TestCase+SubCase+LayerCase+' LES w2 equation vs B/w2, average over '+taverage+' mn'
3052	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_w2_Bw2.ps'
3053	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3054	CS, SCALE=28
3055	GSET, YMIN=-0.1, YMAX=0.1, XMIN=-0.1, XMAX=0.1, TITLE=title_user
3056	cols=INDGEN(1)+2
3057	GPLOT,Y=what_I_plot, X=w2_param(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30, SYM=5, /NOLINES
3058	AXES, XSTEP = 0.1 , XTITLE='aB/w2 - b', YSTEP=0.1, YTITLE='0.5*(dw2/dz)/w2 + E/fm',NDECS=4
3059
3060	FOR l=0, nttot-1 DO BEGIN
3061	oplot, w2_param(,l), full_dw2_w2(,l),thick=0.05,psym=1
3062	ENDFOR
3063	oplot, w2_param(,lt_plotindex_les),w2_param(,lt_plotindex_les),thick=0.05,color=7
3064
3065	PSCLOSE, /NOVIEW
3066
3067	spawn, 'ps2png '+filename
3068
3069	; --- PLOTTING : 0.5*dwu2/dz
3070	print, '........ ///////////// starting local thermal model ///////////'
3071
3072	print, ' -> alimentation'
3073	a_star = make_array(nz, value=0.)
3074	f_star = make_array(nz, value=0.)
3075	f_star(0) = 1.
3076	teta_est = make_array(nz, value=0.)
3077	teta_p = make_array(nz, value=0.)
3078	zw2 = make_array(nz,value=0.)
3079	w_est = make_array(nz,value=0.)
3080	entr_star = make_array(nz,value=0.)
3081	detr_star = make_array(nz,value=0.)
3082	zbuoy_est = make_array(nz,value=0.)
3083	zbuoy = make_array(nz,value=0.)
3084	a_star_tot = 0.
3085
3086	zw2(1)= 0.40.38115522grav(teta_les(0,lt_plotindex_les)/teta_les(1,lt_plotindex_les) -1.)*approx_zdz_les(0)
3087	w_est(1) = zw2(1)
3088	FOR k=0, nz-2 DO BEGIN
3089	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
3090	a_star(k) = MAX([(teta_les(k,lt_plotindex_les)-teta_les(k+1,lt_plotindex_les)),0.])*sqrt(altitudes_LES(k))
3091	a_star_tot = a_star_tot + a_star(k)
3092	lalim = k+1.
3093	endif
3094	ENDFOR
3095	FOR k=0, nz-1 DO BEGIN
3096	a_star(k) = a_star(k)/a_star_tot
3097	ENDFOR
3098	print, 'alimentation :'
3099	;print, a_star
3100	f_star(0)=0.
3101	f_star(1) = a_star(0)
3102	teta_p(0) = teta_les(0,lt_plotindex_les)
3103	teta_est(0) = teta_les(0,lt_plotindex_les)
3104	print, ' -> plume'
3105	FOR k=1, nz-2 DO BEGIN
3106	if (k LT lalim) then begin
3107	teta_est(k) = (f_star(k)teta_p(k-1)+a_star(k)0.25*(teta_les(k,lt_plotindex_les) + teta_p(k-1)))/(f_star(k) + a_star(k))
3108	endif else begin
3109	teta_est(k) = teta_p(k-1)
3110	endelse
3111	zbuoy_est(k) = grav*(teta_est(k)/teta_les(k,lt_plotindex_les) -1.)
3112	zw2fact=fact_epsilon2.approx_zdz_les(k)/(1.+betalpha)
3113	zdw2=afact*zbuoy_est(k)/fact_epsilon
3114	; w_est(k+1) = MAX([0.0001,exp(-zw2fact)*(w_est(k)-zdw2)+zdw2])
3115	w_est(k+1) = MAX([0.0001,w_est(k)(1.-2.approx_zdz_les(k)0.5)+2.1.approx_zdz_les(k)zbuoy_est(k)])
3116	if (w_est(k+1) lt 0.) then begin
3117	w_est(k+1)=zw2(k)
3118	endif
3119	if (w_est(k+1) gt 0.001) then begin
3120	entr_star(k)=f_star(k)approx_zdz_les(k)(betalpha/(1.+betalpha))MAX([0.,afactzbuoy_est(k)/w_est(k+1)])
3121	detr_star(k)=f_star(k)approx_zdz_les(k)MAX([detr_min,-afact(betalpha/(1.+betalpha))zbuoy_est(k)/w_est(k+1)])
3122	endif
3123	if (k lt lalim) then begin
3124	a_star(k)=max([a_star(k),entr_star(k)])
3125	entr_star(k)=0.
3126	endif
3127	if (w_est(k+1) gt 0.001) then begin
3128	f_star(k+1)=f_star(k)+a_star(k)+entr_star(k)-detr_star(k)
3129	if (k lt lalim) then begin
3130	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))
3131	endif else begin
3132	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))
3133	endelse
3134	zbuoy(k) = grav*(teta_p(k)/teta_les(k,lt_plotindex_les) -1.)
3135	zw2fact=fact_epsilon2.approx_zdz_les(k)/(1.+betalpha)
3136	zdw2=afact*zbuoy(k)/fact_epsilon
3137	; zw2(k+1) = MAX([0.0001,exp(-zw2fact)*(zw2(k)-zdw2)+zdw2])
3138	zw2(k+1) = MAX([0.0001,zw2(k)(1.-2.approx_zdz_les(k)0.5)+2.1.approx_zdz_les(k)zbuoy_est(k)])
3139	endif
3140	ENDFOR
3141	print, ' -> done'
3142
3143	print, '........ CHECKING VERTICAL VELOCITY FORMULATION'
3144	what_I_plot = make_array(nz,value=0.)
3145	what_I_overplot = make_array(nz,value=0.)
3146	FOR k=0, nz-2 DO BEGIN
3147	what_I_plot(k) = 0.5*(sqrt(zw2(k)) + sqrt(zw2(k+1)))
3148	ENDFOR
3149	FOR k=0, nZmx-2 DO BEGIN
3150	what_I_overplot(k) = 0.5*(zw2_lev(k,lt_plotindex_gcm) + zw2_lev(k+1,lt_plotindex_gcm))
3151	ENDFOR
3152	labels=['zw2 in les calc as in TH']
3153	title_user = TestCase+SubCase+LayerCase+' LES vertical velocity formulation check'
3154	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_w2_check.ps'
3155	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3156	CS, SCALE=28
3157	GSET, XMIN=0., XMAX=8., YMIN=0., YMAX=7, TITLE=title_user
3158	cols=INDGEN(1)+2
3159	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3160	AXES, XSTEP = 1 , XTITLE='w2 in LES from TH calc m/s', YSTEP=1., YTITLE='Altitude (km)',NDECS=3
3161
3162	oplot, w_mean1(*,lt_plotindex_les), altitudes_LES/1000.
3163	oplot, what_I_overplot, altitudes_GCM/1000.,psym =4
3164
3165	PSCLOSE, /NOVIEW
3166
3167	spawn, 'ps2png '+filename
3168
3169	print, '........ CHECKING TETA ESTIMATIONS FORMULATION'
3170
3171	what_I_plot = [[teta_est],[teta_p],[tplume1moy(*,lt_plotindex_les)]]
3172	labels=['LES estimated teta','LES teta plume calc as in TH','LES teta plume']
3173	title_user = TestCase+SubCase+LayerCase+' LES estimated teta formulation check'
3174	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_teta_check.ps'
3175	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3176	CS, SCALE=28
3177	GSET, XMIN=214., XMAX=220., YMIN=0., YMAX=7, TITLE=title_user
3178	cols=INDGEN(3)+2
3179	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3180	AXES, XSTEP = 1 , XTITLE='Teta plume and Est in LES from TH calc K', YSTEP=1., YTITLE='Altitude (km)',NDECS=3
3181
3182	oplot, teta_gcm(,lt_plotindex_gcm)(buoyancy_gcm(*,lt_plotindex_gcm)/grav +1.), altitudes_GCM/1000.
3183	oplot, teta_gcm(,lt_plotindex_gcm)(buoyancy_est_gcm(*,lt_plotindex_gcm)/grav +1.), altitudes_GCM/1000.
3184
3185	PSCLOSE, /NOVIEW
3186
3187	spawn, 'ps2png '+filename
3188
3189	print, '........ CHECKING MASS FLUX FORMULATION'
3190
3191	what_I_plot = [[f_star/MAX(f_star)],[smoothed_fm_trac1_les(,lt_plotindex_les)/MAX(smoothed_fm_trac1_les(,lt_plotindex_les))]]
3192	labels=['LES normalized f_star ','LES normalized updraft mass flux']
3193	title_user = TestCase+SubCase+LayerCase+' LES normalized f_star formulation check'
3194	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fm_check.ps'
3195	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3196	CS, SCALE=28
3197	GSET, XMIN=0., XMAX=1., YMIN=0., YMAX=7, TITLE=title_user
3198	cols=INDGEN(2)+2
3199	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3200	AXES, XSTEP = 0.1 , XTITLE='f/max(f) in LES from TH calc', YSTEP=1., YTITLE='Altitude (km)',NDECS=3
3201
3202	oplot,fm_therm_gcm_interlay(,lt_plotindex_gcm)/MAX(fm_therm_gcm_interlay(,lt_plotindex_gcm)), altitudes_GCM/1000.,psym=4
3203
3204	PSCLOSE, /NOVIEW
3205
3206	spawn, 'ps2png '+filename
3207
3208	; COMPUTING THE CONTINUITY EQUATION IN THE QUASI-BOUSSINESQ APPROX
3209
3210	da_dt = make_array(nz,n_elements(localtime))
3211	smoothed_da_dt = make_array(nz)
3212	FOR k=0, nz-1 DO BEGIN
3213	da_dt(k,) = deriv(localtime,reform(alpha1out(k,)))/3700.
3214	ENDFOR
3215	FOR t=-ns,ns DO BEGIN
3216	smoothed_da_dt = smoothed_da_dt + REFORM(da_dt(*,lt_plotindex_les+t))
3217	ENDFOR
3218	smoothed_da_dt = smoothed_da_dt/nstot
3219
3220	smoothed_rho = make_array(nz)
3221	FOR t=-ns,ns DO BEGIN
3222	smoothed_rho = smoothed_rho + REFORM(rho(*,lt_plotindex_les+t))
3223	ENDFOR
3224	smoothed_rho = smoothed_rho/nstot
3225
3226	;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
3227	continuity1 = smoothed_rhosmoothed_da_dt + df_dz_les1 - smoothed_e_rate_trac1_lessmoothed_fm_trac1_les + smoothed_d_rate_trac1_les*smoothed_fm_trac1_les
3228
3229	print, '........ CONTINUITY CHECK'
3230
3231	;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]]
3232	what_I_plot = [[continuity1],[smoothed_rhosmoothed_da_dt],[df_dz_les1],[-smoothed_e_rate_trac1_lessmoothed_fm_trac1_les],[smoothed_d_rate_trac1_les*smoothed_fm_trac1_les]]
3233	labels=['total continuity','rho*da/dt','df/dz','-E','D']
3234	title_user = TestCase+SubCase+LayerCase+' LES UDE continuity check, average over '+taverage+' mn'
3235	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_continuity.ps'
3236	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3237	CS, SCALE=28
3238	GSET, XMIN=-0.0001, XMAX=0.0001, YMIN=0, YMAX=10, TITLE=title_user
3239	cols=INDGEN(5)+2
3240	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3241	AXES, XSTEP = 0.00001 , XTITLE='kg.m-2.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
3242	PSCLOSE, /NOVIEW
3243
3244	spawn, 'ps2png '+filename
3245
3246	; COMPUTING THE E-D TERM FROM THE CONTINUITY EQUATION
3247
3248	eminusd1=make_array(nz)
3249	FOR k=0, nz-1 DO BEGIN
3250	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.
3251	ENDFOR
3252	what_I_plot = eminusd1
3253	labels=['e-d']
3254	title_user = TestCase+SubCase+LayerCase+' LES e-d, average over '+taverage+' mn'
3255	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_EminusD.ps'
3256	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3257	CS, SCALE=28
3258	GSET, XMIN=-0.002, XMAX=0.002, YMIN=0, YMAX=10, TITLE=title_user
3259	cols=INDGEN(1)+2
3260	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3261	AXES, XSTEP = 0.0005 , XTITLE='kg.m-2.s-1', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
3262	PSCLOSE, /NOVIEW
3263
3264	spawn, 'ps2png '+filename
3265
3266	; COMPUTING THE TURBULENT FLUX DECOMPOSITION IN PASSIVE ENV AND ACTIVE PLUME
3267	; TO CHECK CONSISTENCY
3268
3269	print, '........ CHECKING CONSISTENCY OF UPDRAFT/ENV DECOMPOSITION'
3270
3271	smoothed_hf1_term1 = make_array(nz)
3272	smoothed_hf1_term2 = make_array(nz)
3273	smoothed_hf1_term3 = make_array(nz)
3274	smoothed_wt = make_array(nz)
3275	FOR t=-ns,ns DO BEGIN
3276	smoothed_hf1_term1 = smoothed_hf1_term1 + REFORM(hf1_term1(*,lt_plotindex_les+t))
3277	smoothed_hf1_term2 = smoothed_hf1_term2 + REFORM(hf1_term2(*,lt_plotindex_les+t))
3278	smoothed_hf1_term3 = smoothed_hf1_term3 + REFORM(hf1_term3(*,lt_plotindex_les+t))
3279	smoothed_wt = smoothed_wt + REFORM(wt(*,lt_plotindex_les+t))
3280	ENDFOR
3281	smoothed_hf1_term1 = smoothed_hf1_term1/nstot
3282	smoothed_hf1_term2 = smoothed_hf1_term2/nstot
3283	smoothed_hf1_term3 = smoothed_hf1_term3/nstot
3284	smoothed_wt = smoothed_wt/nstot
3285
3286	what_I_plot = [[smoothed_hf1_term1],[smoothed_hf1_term2],[smoothed_hf1_term3],[smoothed_hf1_term1+smoothed_hf1_term2+smoothed_hf1_term3]]
3287	labels=['within plume turbulence','within env. turbulence','organized turbulence','TOTAL']
3288	title_user = TestCase+SubCase+LayerCase+' LES turbulence decomposition check, average over '+taverage+' mn'
3289	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_turbu.ps'
3290	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3291	CS, SCALE=28
3292	GSET, XMIN=-1, XMAX=1.5, YMIN=0, YMAX=6, TITLE=title_user
3293	cols=INDGEN(4)+2
3294	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
3295	AXES, XSTEP = 0.5 , XTITLE='m.K/s', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
3296	oplot, smoothed_wt, altitudes_LES/1000.,psym=3
3297	PSCLOSE, /NOVIEW
3298	spawn, 'ps2png '+filename
3299
3300
3301	; COMPUTING THE TURBULENT FLUX DECOMPOSITION IN PASSIVE ENV ,ACTIVE PLUME and ACTIVE DOWNDRAFT
3302	; TO CHECK CONSISTENCY
3303
3304	print, '........ CHECKING CONSISTENCY OF UPDRAFT/DOWNDRAFT/ENV DECOMPOSITION'
3305
3306	smoothed_hf1_ude_term1 = make_array(nz)
3307	smoothed_hf1_ude_term2 = make_array(nz)
3308	smoothed_hf1_ude_term3 = make_array(nz)
3309	smoothed_hf1_ude_term4 = make_array(nz)
3310	FOR t=-ns,ns DO BEGIN
3311	smoothed_hf1_ude_term1 = smoothed_hf1_ude_term1 + REFORM(hf1_ude_term1(*,lt_plotindex_les+t))
3312	smoothed_hf1_ude_term2 = smoothed_hf1_ude_term2 + REFORM(hf1_ude_term2(*,lt_plotindex_les+t))
3313	smoothed_hf1_ude_term3 = smoothed_hf1_ude_term3 + REFORM(hf1_ude_term3(*,lt_plotindex_les+t))
3314	smoothed_hf1_ude_term4 = smoothed_hf1_ude_term4 + REFORM(hf1_ude_term4(*,lt_plotindex_les+t))
3315	ENDFOR
3316	smoothed_hf1_ude_term1 = smoothed_hf1_ude_term1/nstot
3317	smoothed_hf1_ude_term2 = smoothed_hf1_ude_term2/nstot
3318	smoothed_hf1_ude_term3 = smoothed_hf1_ude_term3/nstot
3319	smoothed_hf1_ude_term4 = smoothed_hf1_ude_term4/nstot
3320
3321	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]]
3322	labels=['within plume turbulence','within downdraft turbulence','within env. turbulence','organized turbulence','TOTAL']
3323	title_user = TestCase+SubCase+LayerCase+' LES UDE turbulence decomposition check, average over '+taverage+' mn'
3324	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_turbu_ude.ps'
3325	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3326	CS, SCALE=28
3327	GSET, XMIN=-1, XMAX=1.5, YMIN=0, YMAX=6, TITLE=title_user
3328	cols=INDGEN(5)+2
3329	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
3330	AXES, XSTEP = 0.25 , XTITLE='m.K/s', YSTEP=0.5, YTITLE='Altitude (km)',NDECS=3
3331	oplot, smoothed_wt, altitudes_LES/1000.,psym=3
3332	PSCLOSE, /NOVIEW
3333
3334	spawn, 'ps2png '+filename
3335
3336	print, '........ UDE TRANSPORT : Q'
3337
3338	what_I_plot = [[wq(,lt_plotindex_les)],[wq_updraft(,lt_plotindex_les)],[wq_downdraft(*,lt_plotindex_les)]]
3339
3340	labels=['total turbulent transport','within updraft turbulent transport','within downdraft turbulent transport']
3341	title_user = TestCase+SubCase+LayerCase+' LES UDE turbulence Q transport, average over '+taverage+' mn'
3342	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Qturbu_ude.ps'
3343	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3344	CS, SCALE=28
3345	GSET, XMIN=-1, XMAX=1.5, YMIN=0, YMAX=6, TITLE=title_user
3346	cols=INDGEN(3)+2
3347	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=9, COL=cols, LABELS=labels, THICK = 30
3348	AXES, XSTEP = 0.25 , XTITLE='m.K/s', YSTEP=0.5, YTITLE='Altitude (km)',NDECS=3
3349	PSCLOSE, /NOVIEW
3350
3351	spawn, 'ps2png '+filename
3352
3353
3354
3355	; CHECK CONSISTENCY OF We = W and THETA e = THETA approximation
3356
3357	print, '........ CHECKING CONSISTENCY OF env variable (w_e,theta_e) = mean variable (w_overbar,theta_overbar)'
3358	print, 'as well as mean(w) = 0, in the UDE decomposition'
3359
3360	smoothed_delta_theta_ude = make_array(nz)
3361	smoothed_delta_w_ude = make_array(nz)
3362	smoothed_w_mean1_full = make_array(nz)
3363
3364	FOR t=-ns,ns DO BEGIN
3365	smoothed_delta_theta_ude = smoothed_delta_theta_ude + REFORM(w_mean1_env_ude(,lt_plotindex_les+t)-w_mean1_full(,lt_plotindex_les+t))
3366	smoothed_delta_w_ude = smoothed_delta_w_ude + REFORM(tenv1moy_ude(,lt_plotindex_les+t)-tmoy_full(,lt_plotindex_les+t))
3367	smoothed_w_mean1_full = smoothed_w_mean1_full + REFORM(w_mean1_full(*,lt_plotindex_les+t))
3368	ENDFOR
3369
3370	smoothed_delta_theta_ude = smoothed_delta_theta_ude/nstot
3371	smoothed_delta_w_ude = smoothed_delta_w_ude/nstot
3372	smoothed_w_mean1_full = smoothed_w_mean1_full/nstot
3373
3374	what_I_plot = [[smoothed_delta_theta_ude],[smoothed_delta_w_ude],[smoothed_w_mean1_full]]
3375	labels=['theta env_ude - theta moy','w env_ude - w moy','mean w over domain (,,k)']
3376	title_user = TestCase+SubCase+LayerCase+' LES UDE env/mean approximation check, average over '+taverage+' mn'
3377	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_approx_ude.ps'
3378	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3379	CS, SCALE=28
3380	GSET, XMIN=-2, XMAX=2, YMIN=0, YMAX=10, TITLE=title_user
3381	cols=INDGEN(3)+2
3382	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3383	AXES, XSTEP = 0.5 , XTITLE='(m/s) and (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
3384	PSCLOSE, /NOVIEW
3385
3386	spawn, 'ps2png '+filename
3387
3388	; GETTING SOME INSIGHT ON PLUME'S INSIDE TEMPERATURES
3389
3390	print, '........ STRUCTURE POTENTIAL TEMPERATURES'
3391	xmin = 210
3392	xmax = 220
3393	if (TestCase eq 'Case_Z') then begin
3394	xmin = 260
3395	xmax = 270
3396	endif
3397	ztva = teta_gcm*(buoyancy_gcm/grav +1.)
3398	ztva_est = teta_gcm*(buoyancy_est_gcm/grav +1.)
3399	what_I_plot = [[tplume1moy(,lt_plotindex_les)],[tenv1moy(,lt_plotindex_les)],[teta_les(*,lt_plotindex_les)]]
3400	labels=['Teta updraft','Teta env ','Teta moy']
3401	title_user = TestCase+SubCase+LayerCase+' LES Teta in the structures, no average'
3402	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fullTeta.ps'
3403	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3404	CS, SCALE=28
3405	GSET, XMIN=xmin, XMAX=xmax, YMIN=0, YMAX=6, TITLE=title_user
3406	cols=INDGEN(3)+2
3407	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3408	AXES, XSTEP = 1 , XTITLE='Potential Temperature (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
3409	oplot, ztva(*,lt_plotindex_gcm), altitudes_GCM/1000., thick=0.3
3410	oplot, ztva_est(*,lt_plotindex_gcm), altitudes_GCM/1000., thick=0.3
3411	oplot, teta_gcm(*,lt_plotindex_gcm), altitudes_GCM/1000., thick=0.3
3412	PSCLOSE, /NOVIEW
3413	spawn, 'ps2png '+filename
3414
3415
3416	print, '........ UDE STRUCTURE POTENTIAL TEMPERATURES'
3417
3418	what_I_plot = [[tplume1moy(,lt_plotindex_les)],[tdown1moy(,lt_plotindex_les)],[tenv1moy_ude(,lt_plotindex_les)],[teta_les(,lt_plotindex_les)]]
3419	labels=['Teta updraft','Teta downdraft','Teta env (UDE)','Teta moy']
3420	title_user = TestCase+SubCase+LayerCase+' LES UDE Teta in the structures, no average'
3421	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_fullTeta_ude.ps'
3422	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3423	CS, SCALE=28
3424	GSET, XMIN=210, XMAX=220, YMIN=0, YMAX=6, TITLE=title_user
3425	cols=INDGEN(4)+2
3426	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3427	AXES, XSTEP = 1 , XTITLE='Potential Temperature (K)', YSTEP=1, YTITLE='Altitude (km)',NDECS=4
3428	PSCLOSE, /NOVIEW
3429
3430	spawn, 'ps2png '+filename
3431
3432	; -----------------------------------------------------------------------------------------------------------------------
3433	; End of PLUME diagnostics
3434	; -----------------------------------------------------------------------------------------------------------------------
3435
3436	; * TKE *
3437
3438	print, '........ TKE'
3439
3440	what_I_plot = reform(tke_gcm(*,lt_plotindex_gcm))
3441	labels=['TH tke 1d']
3442	title_user = TestCase+SubCase+LayerCase+' TKE comparison'
3443	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_tke.ps'
3444	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3445	CS, SCALE=28
3446	GSET, XMIN=-1, XMAX=8, YMIN=0, YMAX=10, TITLE=title_user
3447	cols=INDGEN(1)+2
3448	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3449	AXES, XSTEP = 1, XTITLE='Turbulent kinetic energy (kg.m-3)', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
3450
3451	oplot, tke_les(*,lt_plotindex_les), altitudes_LES/1000., psym=4
3452
3453	PSCLOSE, /NOVIEW
3454
3455	spawn, 'ps2png '+filename
3456
3457	; * HEAT FLUX *
3458
3459	print, '........ BL TOP'
3460
3461	dteta_dz_les_full=make_array(nz,nttot)
3462
3463	FOR l=0, nttot-1 DO BEGIN
3464	dteta_dz_les_full(,l)= deriv(altitudes_LES,reform(teta_les(,l)))
3465	ENDFOR
3466
3467	altindex=make_array(nttot,value=0)
3468	altboundary=make_array(nttot)
3469
3470	altindex2=make_array(nttot,value=0)
3471	altboundary2=make_array(nttot)
3472
3473	print, nttot
3474	FOR l=0, nttot-1 DO BEGIN
3475
3476	if (n_elements(where(wt(,l) eq min(wt(,l)))) eq 1) then begin
3477	altindex(l)=where(wt(,l) eq min(wt(,l)))
3478	endif
3479	if (altindex(l) ne -1) then altboundary(l)=altitudes_LES(altindex(l))/1000.
3480	if (l ne 0) then begin
3481	if ((altboundary(l) gt 1.5*altboundary(l-1)) and (localtime(l) gt 10.)) then altboundary(l)=altboundary(l-1)
3482	endif
3483
3484	if (n_elements(where(dteta_dz_les_full(,l) eq max(dteta_dz_les_full(,l)))) eq 1) then begin
3485	altindex2(l)=where(dteta_dz_les_full(,l) eq max(dteta_dz_les_full(,l)))
3486	endif
3487	if (altindex2(l) ne -1) then altboundary2(l)=altitudes_LES(altindex2(l))/1000.
3488
3489	ENDFOR
3490	FOR l=1, nttot-2 DO BEGIN
3491	if (localtime(l) gt 10.) then begin
3492	if ((altboundary2(l) gt 1.5altboundary2(l-1)) and (altboundary2(l+1) lt 1.5altboundary2(l-1))) then altboundary2(l)=(altboundary2(l-1)+altboundary2(l+1))/2.
3493	if (altboundary2(l) gt 1.5*altboundary2(l-1)) then altboundary2(l)=altboundary2(l-1)
3494	endif
3495	if (localtime(l) gt 10.) then begin
3496	if ((altindex2(l) gt 1.5altindex2(l-1)) and (altindex2(l+1) lt 1.5altindex2(l-1))) then altindex2(l)=(altindex2(l-1)+altindex2(l+1))/2.
3497	if (altindex2(l) gt 1.5*altindex2(l-1)) then altindex2(l)=altindex2(l-1)
3498	endif
3499
3500	ENDFOR
3501
3502
3503	altboundary2 = SMOOTH(altboundary2,5,/EDGE_TRUNCATE)
3504	altboundary2 = SMOOTH(altboundary2,10,/EDGE_TRUNCATE)
3505	altindex2 = SMOOTH(altindex2,5,/EDGE_TRUNCATE)
3506	altindex2 = SMOOTH(altindex2,10,/EDGE_TRUNCATE)
3507
3508	FOR l=0, nttot-1 DO BEGIN
3509	if (localtime(l) gt 16.) then altboundary2(l)=0.
3510	if ((localtime(l) lt 9.) or (localtime(l) gt 16.)) then altindex2(l)=0.
3511	ENDFOR
3512
3513	what_I_plot = [[altboundary],[altboundary2]]
3514
3515	labels=['LES Boundary Layer top from Heat flux minimum','LES Boundary Layer top grom Teta gradient']
3516
3517	title_user = TestCase+SubCase+LayerCase+' LES boundary layer top'
3518	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Zi.ps'
3519	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3520	CS, SCALE=28
3521	GSET, XMIN=6, XMAX=18, YMIN=0, YMAX=9, TITLE=title_user
3522	cols=INDGEN(2)+2
3523	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3524	AXES, XSTEP = 1, XTITLE='local time LES', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
3525
3526	;oplot,localtime_gcm,zi_gcm/1000.,thick=0.4,color=5
3527	oplot,localtime_gcm,altitudes_GCM(lmax_gcm)/1000.,thick=0.4,color=6
3528
3529
3530	PSCLOSE, /NOVIEW
3531	spawn, 'ps2png '+filename
3532
3533	what_I_plot = [[altindex],[altindex2]]
3534	labels=['LES Boundary Layer top from Heat flux minimum','LES Boundary Layer top grom Teta gradient']
3535	title_user = TestCase+SubCase+LayerCase+' LES boundary layer top'
3536	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_Ziindex.ps'
3537	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3538	CS, SCALE=28
3539	GSET, XMIN=6, XMAX=18, YMIN=0, YMAX=400, TITLE=title_user
3540	cols=INDGEN(2)+2
3541	GPLOT, X=localtime, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3542	AXES, XSTEP = 1, XTITLE='local time LES', YSTEP=10, YTITLE='Altitude (km)',NDECS=1
3543	PSCLOSE, /NOVIEW
3544	spawn, 'ps2png '+filename
3545
3546	if (1 eq 1) then begin
3547
3548	openw, lun, "input_zipbl", /get_lun
3549	for l=0, nttot-1 do printf, lun, altindex2(l),localtime(l), format='((2x,I0)(4x,F8.2))'
3550	FREE_LUN, lun
3551	close, lun
3552
3553	endif
3554
3555	print, '........ HEAT FLUX'
3556
3557	lay_heatFlux_up = make_array(nZmx,nTmx)
3558	lay_heatFlux_down = make_array(nZmx,nTmx)
3559
3560	FOR k=1, nZmx-1 DO BEGIN
3561	lay_heatFlux_up(k,) = 0.5(heatFlux_up(k,) + heatFlux_up(k-1,))
3562	lay_heatFlux_down(k,) = 0.5(heatFlux_down(k,) + heatFlux_down(k-1,))
3563	ENDFOR
3564	lay_heatFlux_up(0,)=0.5(heatFlux_up(0,*))
3565	lay_heatFlux_down(0,)=0.5(heatFlux_down(0,*))
3566	zkh_gcm_int = make_array(nZmx)
3567
3568	FOR k=0, nZmx-2 DO BEGIN
3569	zkh_gcm_int(k) = 0.5*(zkh(k,lt_plotindex_gcm) + zkh(k+1,lt_plotindex_gcm))
3570	ENDFOR
3571	MY_gcm = -zkh_gcm_int(deriv(altitudes_GCM,reform(zh(,lt_plotindex_gcm)))); - 0.025max(lay_heatFlux_up(,lt_plotindex_gcm)))
3572
3573
3574	;what_I_plot = [[lay_heatFlux_up(,lt_plotindex_gcm)],[MY_gcm],[lay_heatFlux_up(,lt_plotindex_gcm)+MY_gcm]]
3575	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]]
3576
3577	;labels=['TH updraft heat flux','Mellor and Yamada gcm heat flux','Total']
3578	labels=['TH updraft heat flux','Mellor and Yamada gcm heat flux','TH downdraft heat flux','Total']
3579
3580	title_user = TestCase+SubCase+LayerCase+' TH vertical heat flux'
3581	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_WT.ps'
3582	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3583	CS, SCALE=28
3584	GSET, XMIN=-2, XMAX=3, YMIN=0, YMAX=10, TITLE=title_user
3585	;cols=INDGEN(3)+2
3586	cols=INDGEN(4)+2
3587	GPLOT, X=what_I_plot, Y=altitudes_GCM/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3588	AXES, XSTEP = 1, XTITLE='vertical turbulent heat flux', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
3589
3590	;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)
3591
3592	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
3593
3594	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
3595
3596	oplot, smoothed_hf1_ude_term4, altitudes_LES/1000., color = 5
3597
3598	;oplot, smoothed_hf1_term1, altitudes_LES/1000.,thick=0.1,LINESTYLE = 5
3599	;oplot, smoothed_hf1_term2, altitudes_LES/1000.,color=2,thick=0.1,LINESTYLE = 5
3600	;oplot, smoothed_hf1_term3, altitudes_LES/1000.,color=8,thick=0.1,LINESTYLE = 5
3601
3602	oplot, wt(*,lt_plotindex_les), altitudes_LES/1000.
3603
3604	;oplot, smoothed_hf1_term3(*,lt_plotindex_les), altitudes_LES/1000., psym=7
3605
3606	PSCLOSE, /NOVIEW
3607
3608	spawn, 'ps2png '+filename
3609
3610	; * TRACERS *
3611
3612	print, '........ TRACER FLUX'
3613
3614	if (got_tracer_flux eq 'true') then begin
3615
3616	what_I_plot = [[alpha1out(,lt_plotindex_les)wq_updraft(,lt_plotindex_les)],[beta1out(,lt_plotindex_les)wq_downdraft(,lt_plotindex_les)],[(1.-alpha1out(,lt_plotindex_les)-beta1out(,lt_plotindex_les))wq_env_ude(,lt_plotindex_les)],[alpha1out(,lt_plotindex_les)w_mean1(,lt_plotindex_les)(q_mean_up(,lt_plotindex_les) - q_mean(,lt_plotindex_les))+beta1out(,lt_plotindex_les)w_mean1_down(,lt_plotindex_les)(q_mean_down(,lt_plotindex_les) - q_mean(,lt_plotindex_les))],[wq(*,lt_plotindex_les)]]
3617
3618	;labels=['TH updraft heat flux','Mellor and Yamada gcm heat flux','Total']
3619	labels=['LES updraft tracer flux','LES downdraft tracer flux','LES env tracer flux','LES organized structures tracer flux','Total']
3620
3621	title_user = TestCase+SubCase+LayerCase+' LES vertical tracer flux'
3622	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_WQ.ps'
3623	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3624	CS, SCALE=28
3625	GSET, XMIN=-1, XMAX=4, YMIN=0, YMAX=5, TITLE=title_user
3626	;cols=INDGEN(3)+2
3627	cols=INDGEN(5)+2
3628	GPLOT, X=what_I_plot, Y=altitudes_LES/1000., /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3629	AXES, XSTEP = 1, XTITLE='vertical turbulent tracer flux', YSTEP=1, YTITLE='Altitude (km)',NDECS=1
3630
3631	oplot, alpha1out(,lt_plotindex_les)w_mean1(,lt_plotindex_les)(q_mean_up(,lt_plotindex_les) - q_mean(,lt_plotindex_les)), altitudes_LES/1000., color=2,psym=2,thick=0.5
3632	oplot, beta1out(,lt_plotindex_les)w_mean1_down(,lt_plotindex_les)(q_mean_down(,lt_plotindex_les) - q_mean(,lt_plotindex_les)), altitudes_LES/1000., color=3,psym=2,thick=0.5
3633
3634
3635	;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
3636
3637	;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
3638
3639	;oplot, smoothed_hf1_ude_term4, altitudes_LES/1000., color = 5
3640
3641	;oplot, smoothed_hf1_term1, altitudes_LES/1000.,thick=0.1,LINESTYLE = 5
3642	;oplot, smoothed_hf1_term2, altitudes_LES/1000.,color=2,thick=0.1,LINESTYLE = 5
3643	;oplot, smoothed_hf1_term3, altitudes_LES/1000.,color=8,thick=0.1,LINESTYLE = 5
3644
3645	;oplot, wt(*,lt_plotindex_les), altitudes_LES/1000.
3646
3647	;oplot, smoothed_hf1_term3(*,lt_plotindex_les), altitudes_LES/1000., psym=7
3648
3649	PSCLOSE, /NOVIEW
3650
3651	spawn, 'ps2png '+filename
3652
3653	endif
3654
3655	print, '........ TRACER PLOTS DEACTIVATED'
3656
3657	; trying stuff
3658
3659	buoyancy_downdraft = grav*(tdown1moy/tenv1moy_ude-1.)
3660	lmix = make_array(nttot,value=-1.)
3661	altitudes_rel_LES = make_array(nz,nttot)
3662	FOR l=0, nttot-1 DO BEGIN
3663	; kmax = where(w_mean1(,l) eq max(w_mean1(,l)))
3664	; if (kmax(0) ne -1) then lmix(l) = altitudes_LES(kmax(0)) else lmix(l) = -1.
3665	;FOR k=nz-2, 1,-1 DO BEGIN
3666	; 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))
3667	;ENDFOR
3668	FOR k=nz-2, 1,-1 DO BEGIN
3669	if (tdown1moy(k,l) eq 0.) then lmix(l) = altitudes_LES(k)
3670	ENDFOR
3671	ENDFOR
3672
3673	FOR l=0, nttot-1 DO BEGIN
3674	FOR k=0, nz-1 DO BEGIN
3675	altitudes_rel_LES(k,l) = altitudes_LES(k)/lmix(l)
3676	ENDFOR
3677	ENDFOR
3678
3679	print, '........ Teta down / Teta up in UDE'
3680
3681	stuff2=tdown1moy/tenv1moy_ude
3682
3683	what_I_plot = stuff2(*,lt_plotindex_les)
3684	labels=['Teta d/Teta env 12h']
3685	title_user = TestCase+SubCase+LayerCase+' TH trying stuff'
3686	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuff2.ps'
3687	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3688	CS, SCALE=28
3689	GSET, XMIN=0.992, XMAX=1.004, YMIN=0, YMAX=1, TITLE=title_user
3690	cols=INDGEN(1)+2
3691	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3692	AXES, XSTEP = 0.002, XTITLE='Teta d/ Teta env', YSTEP=0.2, YTITLE='Altitude/zi ',NDECS=3
3693
3694	FOR i=0,nttot-1 DO BEGIN
3695	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
3696	ENDFOR
3697	;oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.075)/187.931 + 0.9977, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3698	oplot, (altitudes_rel_LES(,lt_plotindex_les))/19.231 + 0.9938, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3699	oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.075)/187.931 + 0.9982, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3700	oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.60)/(-1333) + 1.00025, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3701	PSCLOSE, /NOVIEW
3702	spawn, 'ps2png '+filename
3703
3704	print, '........ Teta down / Teta up in UDE'
3705
3706	stuff2=tdown1moy/tplume1moy
3707
3708	what_I_plot = stuff2(*,lt_plotindex_les)
3709	labels=['Teta d/Teta u 12h']
3710	title_user = TestCase+SubCase+LayerCase+' TH trying stuff'
3711	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuff2.5.ps'
3712	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3713	CS, SCALE=28
3714	GSET, XMIN=0.95, XMAX=1.1, YMIN=0, YMAX=1, TITLE=title_user
3715	cols=INDGEN(1)+2
3716	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3717	AXES, XSTEP = 0.01, XTITLE='Teta d/ Teta u', YSTEP=0.2, YTITLE='Altitude/zi ',NDECS=3
3718
3719	FOR i=0,nttot-1 DO BEGIN
3720	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
3721	ENDFOR
3722	;oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.075)/187.931 + 0.9977, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3723	;oplot, (altitudes_rel_LES(,lt_plotindex_les))/19.231 + 0.9938, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3724	;oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.60)/(-1333) + 1.00025, altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3725	PSCLOSE, /NOVIEW
3726	spawn, 'ps2png '+filename
3727
3728	print, '........ B down / B up in UDE'
3729
3730	stuff2 = (tdown1moy/tenv1moy_ude -1.)/(tplume1moy/tenv1moy_ude -1.)
3731	what_I_plot = stuff2(*,lt_plotindex_les)
3732	labels=['B down/B up 12h']
3733	title_user = TestCase+SubCase+LayerCase+' TH trying stuff Bd/Bu'
3734	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuffBuBd.ps'
3735	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3736	CS, SCALE=28
3737	GSET, XMIN=-1, XMAX=1., YMIN=0, YMAX=1, TITLE=title_user
3738	cols=INDGEN(1)+2
3739	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3740	AXES, XSTEP = 0.1, XTITLE='B down/ B up', YSTEP=0.1, YTITLE='Altitude/zi ',NDECS=1
3741
3742	FOR i=0,nttot-1 DO BEGIN
3743	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
3744	ENDFOR
3745	;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
3746	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
3747	oplot, (altitudes_rel_LES(,lt_plotindex_les)-0.7)/1., altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3748	;oplot, (altitudes_rel_LES(,lt_plotindex_les))/0.08333333-1., altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3749	oplot, sqrt(altitudes_rel_LES(,lt_plotindex_les)/0.122449)-1., altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3750	PSCLOSE, /NOVIEW
3751	spawn, 'ps2png '+filename
3752
3753	print, '........ F down / F up in UDE'
3754
3755	stuff2 = downward_flux1/fm_trac1_les
3756	what_I_plot = stuff2(*,lt_plotindex_les)
3757	labels=['f down/f up 12h']
3758	title_user = TestCase+SubCase+LayerCase+' TH trying stuff f down/f up'
3759	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stufffufd.ps'
3760	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3761	CS, SCALE=28
3762	GSET, XMIN=-8, XMAX=0.5, YMIN=0, YMAX=1, TITLE=title_user
3763	cols=INDGEN(1)+2
3764	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3765	AXES, XSTEP = 0.5, XTITLE='f down/ f up', YSTEP=0.1, YTITLE='Altitude/zi ',NDECS=1
3766
3767	FOR i=0,nttot-1 DO BEGIN
3768	if(lmix(i) ne -1) then oplot, stuff2(,i), altitudes_rel_LES(,i), thick=0.1
3769	ENDFOR
3770	oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.0149259)/0.00333)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3771	;oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.02)/0.006)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3772	PSCLOSE, /NOVIEW
3773	spawn, 'ps2png '+filename
3774
3775	print, '........ dFdz down / dzFdz up in UDE'
3776	stuff3=make_array(nz,nttot)
3777	FOR l=0, nttot-1 DO BEGIN
3778	stuff3(,l) = deriv(altitudes_LES,downward_flux1(,l))/deriv(altitudes_LES,fm_trac1_les(*,l))
3779	ENDFOR
3780	what_I_plot = stuff3(*,lt_plotindex_les)
3781	labels=['dfdz down/dfdz up 12h']
3782	title_user = TestCase+SubCase+LayerCase+' TH trying stuff f down/f up'
3783	filename = TestCase+SubCase+LayerCase+'Gcm_Les_Comp_stuffdfudfd.ps'
3784	PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3785	CS, SCALE=28
3786	GSET, XMIN=-30, XMAX=30, YMIN=0, YMAX=1, TITLE=title_user
3787	cols=INDGEN(1)+2
3788	GPLOT, X=what_I_plot, Y=altitudes_rel_LES(*,lt_plotindex_les), /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3789	AXES, XSTEP = 5, XTITLE='dfdz down/ dfdz up', YSTEP=0.1, YTITLE='Altitude/zi ',NDECS=1
3790
3791	FOR i=0,nttot-1 DO BEGIN
3792	if(lmix(i) ne -1) then oplot, stuff3(,i), altitudes_rel_LES(,i), thick=0.1
3793	ENDFOR
3794	;oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.0149259)/0.00333)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3795	;oplot, -alog(((altitudes_rel_LES(,lt_plotindex_les)+0.02)/0.006)), altitudes_rel_LES(,lt_plotindex_les),thick=0.3,color=7
3796	PSCLOSE, /NOVIEW
3797	spawn, 'ps2png '+filename
3798	;
3799	;what_I_plot = [[ar_col],[co2_col],[tke_col]]
3800	;labels=['Ar deviation','Co2 deviation','TKE deviation']
3801	;title_user = TestCase+SubCase+' TH 1d tracer conservation'
3802	;filename = TestCase+SubCase+'Gcm_Les_Comp_tracer.ps'
3803	;PSOPEN, THICK=200, CHARSIZE=120, FILE = filename, FONT = 5, TFONT = 5
3804	;CS, SCALE=28
3805	;GSET, XMIN=localtime_gcm(0), XMAX=localtime_gcm(nTmx-1), YMIN=0, YMAX=2, TITLE=title_user
3806	;cols=INDGEN(3)+2
3807	;GPLOT, X=localtime_gcm, Y=what_I_plot, /LEGEND, LEGPOS=1, COL=cols, LABELS=labels, THICK = 30
3808	;AXES, XSTEP = 2, XTITLE='Local time (h)', YSTEP=0.1, YTITLE='Tracer integrated column mass deviation from origin (%)',NDECS=1
3809	;
3810	;PSCLOSE, /NOVIEW
3811	;
3812	;spawn, 'ps2png '+filename
3813	;
3814	;title_user = TestCase+SubCase+' TH argon propagation from first layer'
3815	;PS_START, file = TestCase+SubCase+'Gcm_Les_Comp_Argon.ps'
3816	;
3817	;what_I_plot = transpose(ar(0:6,*))
3818	;
3819	;maxfield_init = 0.05
3820	;minfield_init = 0
3821	;pal=33
3822	;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
3823	;lim_min = minfield_init & w=where(what_I_plot le lim_min) & if (w[0] ne -1) then what_I_plot[w]=lim_min
3824	;
3825	;section, $
3826	; what_I_plot, $ ; 2D field
3827	; localtime_gcm, $ ; horizontal coordinate
3828	; altitudes_gcm(0:6), $ ; altitude coordinate
3829	; minfield=minfield_init, $ ; minimum value of plotted field (=0: calculate)
3830	; maxfield=maxfield_init, $ ; maximum value of plotted field (=0: calculate)
3831	;; minspace=minspace, $ ; minimum value of space window (=0: calculate)
3832	;; maxspace=maxspace, $ ; maximum value of space window (=0: calculate)
3833	;; overcontour=overcontour, $ ; another 2D field to overplot with contour lines (=0: no)
3834	;; overvector_x=overvector_x, $ ; wind vector - x component (=0: no)
3835	;; overvector_y=overvector_y, $ ; wind vector - y component (=0: no)
3836	;; colors=colors, $ ; number of colors/levels (32 is default)
3837	; title_plot=title_user, $ ; title of the plot ('Profile' is default)
3838	; title_axis=['Martian hour (h)','Height above ground (m)'], $ ; title of the [x,y] axis (['Field','Altitude'] is default)
3839	; ct=pal, $ ; color table (33-rainbow is default)
3840	;; topo=topography, $
3841	; format=format ; format of colorbar annotations ('(F6.2)' is default)
3842	;
3843	;PS_END, /PNG
3844	;
3845	;INTERVAL_VOLUME, supermask1, 0.5, 1.,verts, conn
3846	;conn = TETRA_SURFACE(verts, conn)
3847	;oRain = OBJ_NEW('IDLgrPolygon', verts, POLYGONS=conn, $
3848	; COLOR=[255,255,255], SHADING=1)
3849	;XOBJVIEW, oRain, BACKGROUND=[150,200,255]
3850
3851	;INTERVAL_VOLUME, supermask2, 0.5, 1.5,verts, conn
3852	;conn = TETRA_SURFACE(verts, conn)
3853	;oRain = OBJ_NEW('IDLgrPolygon', verts, POLYGONS=conn, $
3854	; COLOR=[255,255,255], SHADING=1)
3855	;XOBJVIEW, oRain, BACKGROUND=[150,200,255]
3856
3857	ENDELSE
3858
3859	print, ''
3860	print, '........ ALL DONE'
3861	print, ''
3862
3863	END

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

Download in other formats:

Original Format