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

;
; faire avant api avec kappa 3.9
;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
path='../TESTGW/'
path='/tmp7/aslmd/'
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'
file=path+experiment+'_wrfout_d01_9999-09-09_09: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
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;
;
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
;
;
;
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))

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   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.]
        ;hache = 8.
        ;column = -hache*alog(reform(columnp(nx,1,*,nt))/610.)
        ;overplot_column = -hache*alog(pressure/610.)
;column = reform(columnp(nx,1,*,nt))
;overplot_column = pressure
;alt = [1.e2,1.e-3]
;;;; mettre profile en log
     ;;;;;;;;;;;
     ;;;;;;;;;;;
     mention = ''
     SPAWN, '\rm param_plot.idl'
     minfield_init = 80.  ;100.
     maxfield_init = 190. ;250.
     title_axis = ['Temperature (K)','Altitude (km)']
;title_axis = ['Temperature (K)','Pressure (Pa)']
     overplot = mcd_tpot
     ;;;;;;;;;;;

;; 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,1,*,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
end