pro legende , n_color,color,niveau,position_x_fin,position_y_debut,position_y_fin,unit

;*************************************************************************************
;This function is for ploting the color scale in the left part of the windows.
;*************************************************************************************

        x0 = position_x_fin+0.03
        x1 = x0+0.03
        y1 = position_y_debut
        largeur = (position_y_fin-position_y_debut)/float((n_color-1))
        epaisseur = 0.003

        while largeur*(n_color+1) ge 1. do largeur = largeur - 0.001

        for i=0,(n_color-1) do begin
                y0 = y1
                y1 = y0+largeur
                POLYFILL,[x0-epaisseur,x0-epaisseur,x1+epaisseur,x1+epaisseur],$
                        [y0-epaisseur,y1+epaisseur,y1+epaisseur,y0-epaisseur],color = 0,/normal
                POLYFILL, [x0,x0,x1,x1],[y0,y1,y1,y0],color = color(i),/normal
                if i eq 0 then begin
                        xyouts,x1+0.01,y0-0.005, STRING(FORMAT='(i4)',niveau(i)),SIZE=0.6,/normal,charthick=1.
                endif else begin
                        xyouts,x1+0.01,y0-0.005,STRING(format='(e10.4)',niveau(i)),SIZE=0.6,$
                                /normal,charthick=1.
                endelse
        endfor
        y0 = y1
        y1 = y0+largeur
        xyouts,x1+0.01,y0-0.005,STRING(format='(e10.4)',niveau(i)),SIZE=0.6,/normal,charthick=1.
        XYOUTS,x1+0.01,y1-largeur*0.4,unit,SIZE=0.6,/normal,charthick=1.
        return
end

;*********************************************************************
;programme soit de tracé, soit de conversion du fichier et de lecture
;*********************************************************************

pro geia_edg

 close, /all

 nomin = ' '
 nomout = ' '
 ligne = ' '
 cpt = 0.
 i = 0.
 j = 0.
 voir = 0

 print, 'Juste visualiser ? (1) oui (0) non'
 read, voir

if not voir then begin
 print, 'Nom du fichier d entree ?'
 read, nomin
 print, 'Nom du fichier de sortie ?'
 read, nomout

 openr, 1, nomin
 for i = 0, 9 do readf, 1, ligne
 while not eof(1) do begin
  readf, 1, ligne
  cpt++
 endwhile
 close, 1
 print, 'cpt : ', cpt
 
 lonlat = fltarr(cpt)
 lon = fltarr(cpt)
 lat = fltarr(cpt)
 lonlattmp = 0.
 val = dblarr(12)
 tab_val = dblarr(cpt,12)

 openr, 1, nomin
 openw, 2, nomout
 for i = 0, 9 do begin
  readf, 1, ligne
  printf, 2, '#'+ligne
 endfor
 printf, 2, '#Converted to kg(C)/m**2/s'
 i = 0.
 while not eof(1) do begin
  readf, 1, lonlattmp, val
  lonlat(i) = lonlattmp
  for j = 0,11 do tab_val(i,j) = val(j)
  i++
 endwhile
 close, 1

 jpm = [31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31.]
 

 for j=0., cpt-1 do begin
   lat(j) = floor(lonlat(j)/1000.)-91.+0.5
   lon(j) = lonlat(j)-floor(lonlat(j)/1000.)*1000-181.+0.5
   for i = 0, 11 do tab_val(j,i) = tab_val(j,i)/(1000000.*24.*3600.*jpm(i))
;Special soil NOx :
;  for i = 0, 11 do tab_val(j,i) = tab_val(j,i)*30.*(1.E-6)/(24.*3600.*jpm(i))
 endfor
 
 lonref = fltarr(360)
 latref = fltarr(180)
 tab = dblarr(360,180)
 for i = 0, 359 do lonref(i) = -179.5 + i
 for i = 0, 179 do latref(i) = -89.5 + i

 for i = 0, 359 do begin
  for j = 0, 179 do begin
   ou = where(lonref(i) eq lon and latref(j) eq lat)
   if ou ne -1 then begin
    if max(tab_val(ou,*)) ne 0. then begin
     ligne = strcompress(string(format='(f7.2)',lonref(i))+' '+string(format='(f7.2)',latref(j))$
                        +' '+string(format='(e13.7)',tab_val(ou,0))+' '+string(format='(e13.7)',tab_val(ou,1))$
                        +' '+string(format='(e13.7)',tab_val(ou,2))+' '+string(format='(e13.7)',tab_val(ou,3))$
                        +' '+string(format='(e13.7)',tab_val(ou,4))+' '+string(format='(e13.7)',tab_val(ou,5))$
                        +' '+string(format='(e13.7)',tab_val(ou,6))+' '+string(format='(e13.7)',tab_val(ou,7))$
                        +' '+string(format='(e13.7)',tab_val(ou,8))+' '+string(format='(e13.7)',tab_val(ou,8))$
                        +' '+string(format='(e13.7)',tab_val(ou,9))+' '+string(format='(e13.7)',tab_val(ou,10))$
                        +' '+string(format='(e13.7)',tab_val(ou,11)))
     printf, 2, ligne
    endif 
   endif 
  endfor
 endfor

 close, 2
endif

;*********************
;Partie visualisation
;*********************

 nomin = ' '
 sortie = 0
 som = 0
 nbniv = 0
;unit = 'kg(C)/m**2/s'
;Special soil NOx :
 unit = 'kg(NO)/m**2/s'
 mois = ['jan','feb','march','apr','may','jun',$
         'jully','aug','sept','oct','nov','dec']

 print, 'Nom du fichier a visualiser ?'
 read, nomin
 espece = nomin

 if voir then begin
  ligne = ' '
  cpt = 0.
  i = 0.
  j = 0.
  jpm = [31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31.]
  openr, 1, nomin
  for i = 0,10 do readf, 1, ligne
  while not eof(1) do begin
   readf, 1, ligne
   cpt++
  endwhile
  close, 1
  print, 'cpt : ', cpt
  lon = fltarr(cpt)
  lat = fltarr(cpt)
  lontmp = 0.
  lattmp = 0.
  val = dblarr(12)
  tab_val = dblarr(cpt,12)
  lonref = fltarr(360)
  latref = fltarr(180)
  tab = dblarr(360,180)
  for i = 0, 359 do lonref(i) = -179.5 + i
  for i = 0, 179 do latref(i) = -89.5 + i
 endif

 openr, 1, nomin
 for i = 0,10 do begin
  readf, 1, ligne
 endfor
 i = 0.
 while not eof(1) do begin
  readf, 1, lontmp, lattmp, val
  lon(i) = lontmp
  lat(i) = lattmp
  for j = 0,11 do tab_val(i,j) = val(j)
  i++
 endwhile
 close, 1
 
;**************************************
;Choix du format de sortie du trace
;et de la sommation ou non sur l annee
;**************************************

 print, 'Format : (0) X (1) ps'
 read, sortie
 print, 'Somme sur l annee ? (1) oui (0) non'
 read, som

 loadct, 33
 tvlct, r, g, b, /get
 !p.region = [0.1,0.1,0.9,0.9]
 mydevice = !d.name

 if sortie  then begin
  set_plot, 'PS'
  device,  /landscape, /color
  !p.position = [0.25,0.2,0.95,0.8]
 endif else begin
  device, decompose=0
  !p.background = 255
  !p.color = 0
  !P.position = [0.2,0.2,0.9,0.9]
  window, xsize = 1200., ysize = 600., retain=2
 endelse

;*********************************************************
;Preparation de la legende et des niveaux pour le contour
;*********************************************************

 temp = 0.
 npog = 0
 nmois = 0
 tab_tmp = dblarr(cpt)
 if som then begin
  for j=0, 11 do tab_tmp = tab_tmp + tab_val(*,j)*jpm(j)/365.
 endif else begin
  print, 'Quel mois ?'
  read, nmois
  tab_tmp = tab_val(*,nmois-1)
 endelse

 ou = where(tab_tmp ne 0)
 nbnnul = n_elements(ou)
 non_nul = dblarr(nbnnul)
 for i=0.,nbnnul-1. do non_nul(i) = tab_tmp(ou(i))
 nmax = max(non_nul)
 nmin = min(non_nul)

 ogmin = floor(alog10(nmin))
 ogmax = floor(alog10(nmax) + 0.5)
 nog = ogmax - ogmin + 1
 npog = nbniv / nog
 if npog lt 1 then begin
  print, 'Pas assez : ', npog
  while npog lt 1 do begin
   nbniv++
   npog = nbniv / nog
  endwhile
 endif
 if npog gt 10 then begin
  print, 'Trop : ', npog
  while npog gt 10 do begin
   nbniv--
   npog = nbniv / nog
  endwhile
 endif
 print, ogmin, ogmax, nog
 print, 'Nombre de niveaux : ', nbniv
 print, 'Nombre de niveaux par ordre de grandeur : ', npog
 kmin = floor(alog10(nmin))
 kmax = floor(alog10(nmax+0.5))

 print, 'De : ', nmin,' a : ', nmax
 tab_col = fltarr(nbniv)
 tab_col2 = fltarr(nbniv+1)

 for i=0.,nbniv-1. do tab_col(i)= 250.*(i+1)/ (float(nbniv))
 tab_col2(0)=255.
 for i=0.,nbniv-1. do tab_col2(i+1)= 250.*(i+1)/ (float(nbniv))
 V = dblarr(nbniv+1)
 V2 = dblarr(nbniv+2)

 V2(0) = 0

 k = kmin
 j = 1
 for i=1,nbniv+1 do begin
  V2(i) = j*(1./npog)*(10.D^k)
  j++
  if j ge npog then begin
   j = 1
   k++
  endif
 endfor

 for i=0.,nbniv do V(i) = V2(i+1)

;********************
;Trace des emissions
;********************

 if som then begin
  for j=0, 11 do begin
   for i=0.,cpt-1 do begin
    tab(where(lonref eq lon(i)),where(latref eq lat(i))) = tab(where(lonref eq lon(i)),where(latref eq lat(i))) + tab_val(i,j)*jpm(j)/365.
   endfor
  endfor
  if sortie then begin
   nom = strcompress('geia_somme_'+espece+'.ps')
   device, filename = nom
   map_set,title = espece+'  Min :'+string(nmin)+'  Max :'+string(nmax)+' '+unit
   contour,tab,lonref,latref,levels=V,/cell_fill,c_colors=tab_col,/overplot
   map_continents, /hires
   map_grid,/label,latdel=10,londel=10,latlab=-180,lonlab=-90,charsize=0.5
   legende, nbniv+1, tab_col2, V2, 0.0, 0.01, 0.9, unit
   device, /close
   print, 'Done !!'
  endif else begin
   map_set,title=espece+'  Min :'+string(nmin)+'  Max :'+string(nmax)+' '+unit
   contour,tab,lonref,latref,levels=V,/cell_fill,c_colors=tab_col,/overplot
   map_continents, /hires
   map_grid,/label,latdel=10,londel=10,latlab=-180,lonlab=-90,charsize=0.5
   legende, nbniv+1, tab_col2, V2, 0.0, 0.01, 0.9, unit
  endelse
 endif else begin
  tot_val = dblarr(360,180,12)
  for j=0, 11 do begin
   for i=0.,cpt-1. do begin
    tot_val(where(lonref eq lon(i)),where(latref eq lat(i)),j) = tab_val(i,j)
   endfor
  endfor
  if sortie then begin
   deb_nom = 'geia_mois_'
   fin_nom = '.ps'
   for i=0,11 do begin
    device, filename=strcompress(deb_nom+string(i+1)+'_'+espece+fin_nom, /remove_all)
    map_set,title=espece+'  Min :'+string(nmin)+'  Max :'+string(nmax)+' '+unit+'  '+mois(i)
    contour,tot_val(*,*,i),lonref,latref,levels=V,/cell_fill,c_colors=tab_col,/overplot
    map_continents, /hires
    map_grid, /label,latdel=10,londel=10,latlab=-180,lonlab=-90,charsize=0.5
    legende, nbniv+1, tab_col2, V2, 0.0, 0.01, 0.9, unit
    device, /close
    print, mois(i), 'Fait !!'
   endfor
  endif else begin
   map_set, title=espece+'  Min :'+string(nmin)+'  Max :'+string(nmax)+$
                  ' '+unit+'  '+mois(nmois-1)
   contour,tot_val(*,*,nmois-1),lonref,latref,levels=V,/cell_fill,c_colors=tab_col,/overplot
   map_continents, /hires
   map_grid, /label,latdel=10,londel=10,latlab=-180,lonlab=-90,charsize=0.5
   legende, nbniv+1, tab_col2, V2, 0.0, 0.01, 0.9, unit
  endelse
 endelse
 set_plot, mydevice

stop
end