;
;  Plot an ecmwf_monthly_summary
;
pro plot_ecmwf,site,searchdate
;
;  common containing the ecmwf data
;
common ecmwf,newpress,newuwind,newvwind,newtemp,newshum,newliq,newice,newcfrac,$
  newrhum,newomega,numtimes,base_time,time_offset,numheights,newheight
;
;  File information
;
curdate=strcompress(string(searchdate),/remove_all)
site=strcompress(site,/remove_all)
;
;  Read in the ecmwf data
;
get_ecmwf_data,site,searchdate
;
;  Select which data to plot
;
dataname=strarr(15)
dataname=['Pressure','Zonal Wind Component','Meridional Wind Component',$
          'Temperature','Specific Humidity',$
          'Specific Cloud Liquid Water Content','Specific Cloud Ice Content',$
          'Cloud Fraction','Relative Humidity',$
          'Omega = Vertical Velocity in Pressure Coordinates',$
          'IWC where Cloud Fraction > 0.1',$
          'IWC where Cloud Fraction > 0.8',$
	  'IWC where Cloud Fraction > 0.5',$
          'IWC where Cloud Fraction is 0.1 - 0.5',$
	  'IWC']
datatag=strarr(15)
datatag=['press','uwind','vwind','temp','shum','liq','ice','cfrac','rhum',$
         'omega','ice_cfrac_10','ice_cfrac_80','ice_cfrac_50',$
         'ice_cfrac_10_50','ice_cfrac_0_100']
cb_title=strarr(15)
cb_title=['mb','m/s','m/s','K','kg/kg','kg/kg','kg/kg','percent/100',$
          'percent/100','Pa/s','log mg/m!E3','log mg/m!E3','log mg/m!E3',$
          'log mg/m!E3','log mg/m!E3']
;
;  Loop through the variables in the file
;
;for datanum=0,14 do begin  ;loop through the different variables
for datanum=14,14 do begin

print, datatag[datanum]
if datanum eq 0 then data=newpress
if datanum eq 1 then data=newuwind 
if datanum eq 2 then data=newvwind  
if datanum eq 3 then data=newtemp  
if datanum eq 4 then data=newshum  
if datanum eq 5 then data=newliq 
if datanum eq 6 then data=newice
if datanum eq 7 then data=newcfrac  
if datanum eq 8 then data=newrhum 
if datanum eq 9 then data=newomega  

if datanum ge 10 and datanum le 14 then begin
  data=make_array(numtimes,numheights,value=-99999)
  data=newice*((newpress*100.)/(287.04*newtemp))
;
;  Select the different restrictions for newice based on cfrac
;
  ;Specific Cloud Ice Content > 0.1
  if datanum eq 10 then data[where(newcfrac le 0.1)]=-99999.  
  ;Specific Cloud Ice Content > 0.8
  if datanum eq 11 then data[where(newcfrac le 0.8)]=-99999.  
  ;Specific Cloud Ice Content > 0.5
  if datanum eq 12 then data[where(newcfrac le 0.5)]=-99999.  
  ;Specific Cloud Ice Content 0.1-0.5
  if datanum eq 13 then data[where(newcfrac lt 0.1 or newcfrac gt 0.5)]=-99999.
  ;Specific Cloud Ice Content 0.0-1.0
  if datanum eq 14 then data[where(newcfrac lt 0.0001 or newcfrac ge 1.0)]=-99999.  
  ;Specific Cloud Ice Content 0.0-1.0
  ;if datanum eq 14 then data[where(newcfrac le 0 or newcfrac gt 1.0)]=-99999.  

  for j=0,numtimes-1 do begin
    for k=0,numheights-1 do begin
      if data[j,k] ne -99999. then begin
        data[j,k]=(data[j,k]/newcfrac[j,k])*1.e6
        if data[j,k] lt 1.e-2 then data[j,k]=-99999.
        if data[j,k] gt 0. then data[j,k]=alog10(data[j,k])
      endif
    endfor
  endfor
endif  ;end of datanum ge 10

;
;  Find the max and min data values
;
dmax=max(data)
result=where(data eq -99999, count)
if count ne 0 then data[result]=1e6  
dmin=min(data)
;if count ne 0 then data[result]=-9999  ;turns no data to min value

;  For cloud fraction and ice, turns 0 cloud fraction and ice to white
if datanum eq 7 or datanum eq 6 then begin
  result=where(data eq 0,count)
  if count gt 0 then data[result]=1e6
endif

;dmin=-4
;dmax=3
print, 'dmax',dmax,'dmin',dmin
if min(data) eq max(data) then begin
	print,'no data in this variable'
	goto, jump1
endif

;
;  Calculate the file starting and ending date and time
;
standard_date_time,base_time+time_offset[0],syear,smonth,sday,shour,smin,ssec
print,'file start',syear,smonth,sday,shour,smin,ssec
standard_date_time,base_time+time_offset[numtimes-1],eyear,emonth,eday,ehour,emin,esec
print,'file end',eyear,emonth,eday,ehour,emin,esec
;
;  Calculate the standard date arrays
;
ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)
;
;  Set up the plot
;
loadct,5  ;load the color table
!y.omargin=[2,5]  ;set aside some room at the top of the page for the title
!x.margin=[8,15]  ;change to accomodate the colorbar
!y.margin=[2,2]  ;change to move the images closer together
!p.multi=[0,1,5]  ;put 5 plots on a page
!p.charsize=2  ;increase the character size
!p.background=!d.n_colors-1  ;make the background white
xwin=800  ;xsize of the plot window
ywin=800  ;ysize of the plot window
window,0,xsize=xwin,ysize=ywin,title='ecmwf',/pixmap
dummy=label_date(date_format=[' %D %M !C %H:%I '])
;
;  Plot 1
;
plot1='yes'
if plot1 eq 'yes' then begin

;week1=num_sec_ep(1999,07,1,0,0,0)
week1=num_sec_ep(syear,smonth,1,0,0,0)
result=where(base_time+time_offset ge week1,count)
if count gt 0 then week1=result[0] else goto, no_plot1
;week2=num_sec_ep(1999,07,8,0,0,0)
week2=num_sec_ep(syear,smonth,8,0,0,0)
result=where(base_time+time_offset lt week2,count)
if count gt 0 then week2=result[count-1] else goto, no_plot1

week=data[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,newheight,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
image=congrid(week,sx,sy)
image=bytscl(image,max=dmax,min=dmin)
tv,image,px[0],py[0]
axis,weektime[0],newheight[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[0],yaxis=0,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],newheight[0],yaxis=1,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
endif
no_plot1:
;
;  Plot 2
;
plot2='yes'
if plot2 eq 'yes' then begin

;week1=num_sec_ep(1999,07,8,0,0,0)
week1=num_sec_ep(syear,smonth,8,0,0,0)
result=where(base_time+time_offset ge week1,count)
if count gt 0 then week1=result[0] else goto, no_plot2
;week2=num_sec_ep(1999,07,15,0,0,0)
week2=num_sec_ep(syear,smonth,15,0,0,0)
result=where(base_time+time_offset lt week2,count)
if count gt 0 then week2=result[count-1] else goto, no_plot2

week=data[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,newheight,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
image=bytscl(week,max=dmax,min=dmin)
tv,congrid(image,sx,sy),px[0],py[0]
axis,weektime[0],newheight[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[0],yaxis=0,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],newheight[0],yaxis=1,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
endif
no_plot2:
;
;  Plot 3
;
plot3='yes'
if plot3 eq 'yes' then begin

;week1=num_sec_ep(1999,07,15,0,0,0)
week1=num_sec_ep(syear,smonth,15,0,0,0)
result=where(base_time+time_offset ge week1,count)
if count gt 0 then week1=result[0] else goto,no_plot3
;week2=num_sec_ep(1999,07,22,0,0,0)
week2=num_sec_ep(syear,smonth,22,0,0,0)
result=where(base_time+time_offset lt week2,count)
if count gt 0 then week2=result[count-1] else goto,no_plot3

week=data[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,newheight,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
image=bytscl(week,max=dmax,min=dmin)
tv,congrid(image,sx,sy),px[0],py[0]
axis,weektime[0],newheight[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[0],yaxis=0,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],newheight[0],yaxis=1,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
endif
no_plot3:
;
;  Plot 4
;
plot4='yes'
if plot4 eq 'yes' then begin

;week1=num_sec_ep(1999,07,22,0,0,0)
week1=num_sec_ep(syear,smonth,22,0,0,0)
result=where(base_time+time_offset ge week1,count)
if count gt 0 then week1=result[0] else goto, no_plot4
;week2=num_sec_ep(1999,07,29,0,0,0)
week2=num_sec_ep(syear,smonth,29,0,0,0)
result=where(base_time+time_offset lt week2,count)
if count gt 0 then week2=result[count-1] else goto, no_plot4 

week=data[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,newheight,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
image=congrid(week,sx,sy)
image=bytscl(image,max=dmax,min=dmin)
tv,image,px[0],py[0]
axis,weektime[0],newheight[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[0],yaxis=0,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],newheight[0],yaxis=1,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
endif
no_plot4:
;
;  Plot 5
;
plot5='yes'
if plot5 eq 'yes' then begin

;week1=num_sec_ep(1999,07,29,0,0,0)
week1=num_sec_ep(syear,smonth,29,0,0,0)
result=where(base_time+time_offset ge week1,count)
if count gt 0 then week1=result[0] else goto,no_plot5
;week2=num_sec_ep(1999,07,31,23,0,0)
week2=num_sec_ep(syear,smonth,31,23,59,59)
result=where(base_time+time_offset le week2,count)
if count gt 0 then week2=result[count-1] else goto,no_plot5

week=data[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,newheight,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
if (eday-sday-28) eq 2 then begin ;determine if 30 or 31 days
	adjust=3./7.   ;31 days
	numticks=3
endif else begin 
	adjust=2./7.  ;30 days
	numticks=2
endelse
image=congrid(bytscl(week,max=dmax,min=dmin),sx*adjust,sy)
tv,image,px[0],py[0]
sz=size(image)

contour,week,weektime,newheight,xstyle=1,ystyle=1,/nodata,/noerase,$
	position=[px[0],py[0],px[0]+sz[1]-1,py[0]+sz[2]-1],/device
tv,image,px[0],py[0]
axis,weektime[0],newheight[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],newheight[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xrange=[weektime[0],weektime[n_elements(weektime)-1]],$
  xtickformat='(A1)'
axis,weektime[0],newheight[0],yaxis=0,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],newheight[0],yaxis=1,color=0,ystyle=1,$
  yrange=[newheight[0],newheight[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
endif
no_plot5:
;
;  Add a title
;
jump1:
if smonth eq 1 then smon='Jan'
if smonth eq 2 then smon='Feb'
if smonth eq 3 then smon='Mar'
if smonth eq 4 then smon='Apr'
if smonth eq 5 then smon='May'
if smonth eq 6 then smon='Jun'
if smonth eq 7 then smon='Jul'
if smonth eq 8 then smon='Aug'
if smonth eq 9 then smon='Sep'
if smonth eq 10 then smon='Oct'
if smonth eq 11 then smon='Nov'
if smonth eq 12 then smon='Dec'
if site eq 'sgp' then begin
	capsite='SGP'
	datastream='sgpecmwf'
endif else if site eq 'nsa' then begin
	capsite='NSA'
	datastream='nsaecmwf'
endif else if site eq 'twpc1' then begin
	capsite='TWP  Manus'
	datastream='twpecmwfman'
endif else if site eq 'twpc2' then begin
	capsite='TWP  Nauru'
	datastream='twpecmwfnau'
endif
syr=strcompress(string(syear),/remove_all)
xyouts,0.5,0.95,align=0.5,charsize=1.7,/normal,color=0,$
	'Monthly ECMWF Summary for '+dataname[datanum]+'!C'+capsite+' - '+smon+' '+syr
;
;  Create gif file
;
strmon=string(smonth)
if smonth lt 10 then strmon='0'+strmon
strday=string(sday)
if sday lt 10 then strday='0'+strday
curdate=strcompress(syr+strmon+strday,/remove_all)

gifdir='/data/mace4/arm_data/temp_image/'
giffile=datastream+'.'+curdate+'.'+datatag[datanum]

write_gif,gifdir+giffile+'.gif',tvrd( )
;
;  Remove png file
; 
unix_command='rm '+gifdir+giffile+'.png'
spawn,unix_command

endfor  ;end loop through variables
end