Context Navigation

gettherm.pro @ 638

Last change on this file since 638 was 211, checked in by aslmd, 14 years ago
MESOSCALE_DEV: thermal graphical stuff was lost in the clean-up. this is fixed. plus checked if no file is missing anymore, OK.
File size: 118.2 KB

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

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

Download in other formats:

Original Format