pro gravitwaveprof
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;
; faire avant api avec kappa 3.9
;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	path='../TESTGW/'
	path='/tmp7/aslmd/'
path='/home/aslmd/GRAVITWAVE/'
	experiment='mcd_lt16_rcpcst'
	;experiment='mcd_lt16'
	;experiment='short_mcd_lt16_rcpcst'
	;experiment='short_mcd_lt16'
	;experiment='short_mcd_lt16_calcprho'
	;experiment='mcd_lt16_calcprho'
	;experiment='short_mcd_lt16_calcprho_tot'
	;experiment='colder_gcm_add_short'
	;experiment='colder_gcm_add'
	experiment=''
	;experiment='wind15'
	;experiment='wind20c'
	;experiment='colder_gcm_add_3.9_fine_rad_mountain'
	experiment='test'
experiment='GW_MARS_highwind_3D.157077'
experiment='GW_MARS_highwind_3D_lower_opacity.161436'
experiment='GW_MARS_highwind_3D_loweropacity_lotspoints.161468'
experiment='GW_MARS_highwind_3D_loweropacity_morepoints.161542'
;experiment='GW_MARS_highwind_3D_loweropacity_morepoints_LT15.161544'
;experiment='GW_MARS_highwind_3D_loweropacity_morepoints_LTcst.168440'
;experiment='GW_MARS_highwind_3D_loweropacity_morepoints_widehill.168468'
	file=path+experiment+'_wrfout_d01_9999-09-09_09:00:00_z'
file=path+experiment+'/wrfout_d01_9999-09-09_09:00:00_z'
;file=path+experiment+'/wrfout_d01_9999-09-09_10:00:00_z'
;file=path+experiment+'/wrfout_d01_9999-09-09_11:00:00_z'
charvar='W' & charvarc='W'
charvar='tk' & charvarc='W'
charvar='tk' & charvarc='W' & cond=1
charvar='tpot' & charvarc='tpot'
;charvar='PTOT' & charvarc='PTOT'
charvar='tk' & charvarc='tk'
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
nts=1  & nte=15   & nxs=0   & nxe=59
nts=3  & nte=nts  & nxs=999 & nxe=999 ;; moyenne ;; attention a la montagne !!!
nts=4  & nte=nts  & nxs=31  & nxe=nxs ;; AU CENTRE ! sinon choisir le point en traçant de 0a 59
nts=1  & nte=nts  & nxs=0   & nxe=59
nts=4  & nte=nts  & nxs=0   & nxe=59
nts=1  & nte=15   & nxs=31  & nxe=31
nts=1  & nte=15   & nxs=45  & nxe=45
nts=9  & nte=9    & nxs=0   & nxe=59
nts=1  & nte=13   & nxs=59  & nxe=59
nts=8  & nte=8    & nxs=0   & nxe=121
nts=1  & nte=9    & nxs=59  & nxe=59
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
nts=1  & nte=10  & nxs=50   & nxe=50
;nts=1  & nte=1  & nxs=6   & nxe=6
nts=1  & nte=100  & nxs=50   & nxe=50
nts=1  & nte=100  & nxs=60   & nxe=60
nts=1  & nte=100  & nxs=20   & nxe=20
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;
;
what_I_plot=0. & overcontour=0.
;SPAWN, '\rm param_plot.idl ; cp gravitwave_inc2.pro param_plot.idl'
;
;
;
getcdf, $
        file=file, $
        charvar=charvar, $
        invar=invar
;getcdf, $
;        file=file, $
;        charvar=charvarc, $
;        invar=invarc
getcdf, $
        file=file, $
        charvar='vert', $
        invar=vert
getcdf, $
        file=file, $
        charvar='PTOT', $
        invar=columnp
;;
;; 
;;
	s = size(invar) & middle = floor(s[2]/2) & print, 'PLOT at y subs ', middle
        nte = min(s[4]-1,nte) & print, 'nte ', nte
        mean_invar = TOTAL(invar(*,*,*,nts:nte),4) / float(nte-nts+1)
	;mean_columnp = TOTAL(columnp(*,*,*,nts:nte),4) / float(nte-nts+1)
        anomal_invar = invar & for i=0,s[4]-1 do anomal_invar(*,*,*,i) = anomal_invar(*,*,*,i) - mean_invar
;;
;;
;;
;	goto, perturb
for nt=nts,nte do begin
for nx=nxs,nxe do begin
;
zefile=experiment+'_'+charvar+'_'+string(nt+100,'(I0)')+'_'+string(nx+100,'(I0)')
;
;
;
  PS_Start, filename=zefile+'.ps'
  print, zefile+'.ps'
  !P.Charsize = 1.2
  !p.charthick = 2.0
  !p.thick = 2.0
  !x.thick = 2.0
  !y.thick = 2.0
;
;
;
;;what_I_plot = reform(invar(*,1,*,nt))
;what_I_plot = reform(invar(*,59,*,nt))
what_I_plot = reform(invar(*,middle,*,nt))

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   file=path+experiment+'/input_sounding' & header='' & nlines_header=1 & ncol = 5
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   nlines = FILE_LINES(file)-nlines_header & data=FLTARR(ncol,nlines)
   OPENR, lun, file, /GET_LUN & READF, lun, header & READF, lun, data & CLOSE, lun & FREE_LUN, lun
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   mcd_column      = reform(data(0,*)) 
   calc_tpot       = reform(data(1,*))

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   file=path+experiment+'/input_therm' & ncol = 5
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   nlines = FILE_LINES(file) & data=FLTARR(ncol,nlines)
   OPENR, lun, file, /GET_LUN & READF, lun, data & CLOSE, lun & FREE_LUN, lun
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   pressure        = reform(data(2,*))
   mcd_tpot        = reform(data(4,*))   

   if (nx lt 999) then begin
     what_I_plot = reform(what_I_plot(nx,*))
   endif else begin
     what_I_plot = TOTAL(what_I_plot,1) / n_elements(reform(what_I_plot(*,0)))
   endelse

     ;;;;;;;;;;;
     ;;;;;;;;;;;
     ;;;;;;;;;;;
     column = vert/1000.
     overplot_column = mcd_column/1000.
     alt = [20.,150.]
     alt = [40.,100.]
     alt = [40.,80.]
     alt = [0.,120.]
     ;;;;;;;;;;;
     ;;;;;;;;;;;
     ;mention = ''
     SPAWN, '\rm param_plot.idl'
     minfield_init = 85.  ;100.
     maxfield_init = 185. ;250.
     title_axis = ['Temperature (K)','Altitude (km)']
     overplot = mcd_tpot
     ;;;;;;;;;;;
     SPAWN, "echo 'intervaly=10.' >> param_plot.idl"
     SPAWN, "echo 'intervalx=5.' >> param_plot.idl"
     ;;;;;;;;;;;

	;overplot = reform(mean_invar(nx,middle,*))
	;overplot_column = column

;;;;;;;;;;;;;;;;;;;;;;;;;;;; PLOT PRESSION mettre profile en log
;column = reform(columnp(nx,middle,*,nt))
;overplot_column = pressure
;	;;overplot = mean_invar
;	;;overplot_column = mean_columnp
;alt = [1.e2,1.e-3]
;title_axis = ['Temperature (K)','Pressure (Pa)']
;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PLOT PRESSION

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; calculer la CAPE ? Tsat - Tenv / Tenv


			;what_I_plot = pressure
			;	what_I_plot = calc_tpot
			;overplot = alog(overplot) / alog(10.)
			;overplot = overplot[where(overplot lt 35.)] ;; evite les valeurs impossibles
			;what_I_plot = alog(what_I_plot) / alog(10.)
			;minfield_init = -5.
			;maxfield_init = 3.
			;	;;tpot
			;	minfield_init = 2.
			;	maxfield_init = 4.
			;what_I_plot = 100. * abs(calc_tpot - overplot) / overplot
			;minfield_init = 0.000001
			;maxfield_init = 100.

   profile, $
        what_I_plot, $                          ; 1D vertical profile
        column, $                               ; altitudes     
        alt=alt, $                              ; altitude range [altmin, altmax]
        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
        inprofile=overplot, $                   ; another vertical profile to overplot
        incolumn=overplot_column, $             ; altitudes of the other vertical profile (in case /= column)
;        discrete=discrete, $                    ; show the profile points (= type of points in !psym)
;        title_plot=title_user, $                ; title of the plot ('Profile' is default)
        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
        mention=mention                         ; add text precision within the plot window (default is nothing or '')

	;;; sponge layer
	;w=where(FINITE(overplot) ne 0)
	;oplot, [minfield_init,maxfield_init], [max(overplot_column[w])-50.,max(overplot_column[w])-50.], linestyle=1

	;column = columnp
	;@tempcond.inc
	;yeye = reform(overplot(nx,1,*,nt)) 
	;;w = where(yeye le 0.) & yeye[w] = !VALUES.F_NAN
	;;oplot, yeye, vert/1000., linestyle=1

@tempcond.inc
overplot = reform(overplot(nx,middle,*,nt))
oplot, overplot, overplot_column, linestyle=1

	;calc_tk = calc_tpot * (pressure/610.)^(1/4.4)
	;oplot, calc_tk, column, linestyle=3

PS_End, /PNG
endfor
endfor

perturb:
PS_Start, filename=experiment+'_'+charvar+'_perturb_'+string(nxs+100,'(I0)')+'.ps'

     what_I_plot_loop = reform(anomal_invar(nxs,middle,*,*))

     what_I_plot = reform(what_I_plot_loop(*,nts))
     column = vert/1000.
     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     alt = [0.,120.]
     minfield_init = -30. 
     maxfield_init = 30. 
     title_axis = ['Temperature anomaly (K)','Altitude (km)']
		;; peut pas changer	
		;SPAWN, "echo 'intervaly=10.' >> param_plot.idl"
     		;SPAWN, "echo 'intervalx=5.' >> param_plot.idl"
     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     profile, $
        what_I_plot, $                          ; 1D vertical profile
        column, $                               ; altitudes     
        alt=alt, $                              ; altitude range [altmin, altmax]
        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
        inprofile=overplot, $                   ; another vertical profile to overplot
        incolumn=overplot_column, $             ; altitudes of the other vertical profile (in case /= column)
;        discrete=discrete, $                    ; show the profile points (= type of points in !psym)
;        title_plot=title_user, $                ; title of the plot ('Profile' is default)
        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
        mention=mention                         ; add text precision within the plot window (default is nothing or '')

     for iii=nts+1, nte do begin
        overplot = reform(what_I_plot_loop(*,iii))
        overplot_column = column
        oplot, overplot, overplot_column
     endfor

	;temppp = reform(mean_invar(nxs,middle,*))
	;hache = 191.*temppp/3.72/1000.
	;ok = column / ( 2. * hache )
	;print, exp(ok)
	;yeah = exp(ok)
	;oplot, yeah, column

PS_End, /PNG

end