;---------------------------  
; Reading Data
;---------------------------  
FUNCTION data_read, cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex

@common
@com_eg
 

; same data as previous read ?

read = 1

IF cmd.var EQ cmd_prev.var AND cmd.grid EQ cmd_prev.grid AND cmd.plt EQ cmd_prev.plt AND cmd.timave EQ cmd_prev.timave AND cmd.date1 EQ cmd_prev.date1 THEN BEGIN 
   IF hotyp NE '-' THEN BEGIN 
      ; time serie
      IF cmd.spec EQ cmd_prev.spec THEN BEGIN 
         fldr = field 
         read = 0
      ENDIF ELSE read = 1
      
   ENDIF ELSE BEGIN 
      read = 1
      IF strpos(cmd.plt, '#') NE -1 THEN BEGIN
         IF  cmd.plt EQ cmd_prev.plt THEN BEGIN ; not same zoom
            fldr = field 
            read = 0
         ENDIF ELSE begin
            read = 1
         ENDELSE
      ENDIF 
   ENDELSE 
ENDIF ELSE read = 1

read = 1

IF read EQ 1 THEN BEGIN 

   cmd_prev = cmd

; difference from cmd.exp ? division from exp ?

   diff_from_exp = 0
   div_from_exp =  0

;   argvar = str_sep(cmd.exp,'-')

   IF strpos(cmd.exp, '-') NE -1 THEN BEGIN
      argvar = strsplit(cmd.exp,'-', /EXTRACT)
   ENDIF ELSE IF strpos(cmd.exp, '/') NE -1 THEN BEGIN
      argvar = strsplit(cmd.exp,'/', /EXTRACT)
   ENDIF ELSE argvar =  cmd.exp
      
   IF n_elements(argvar) EQ 2 THEN BEGIN
      IF strpos(cmd.exp, '-') NE -1 THEN BEGIN
         diff_from_exp = 1
         exp_init = cmd.exp
         cmd.exp = argvar[0]
      ENDIF
      IF strpos(cmd.exp, '/') NE -1 THEN BEGIN
         div_from_exp = 1
         exp_init = cmd.exp
         cmd.exp = argvar[0]
      ENDIF      
   ENDIF 
   varexp = cmd.exp   

; define data base
  
   ncdf_db = def_dbase(cmd.exp)

; define file name

   def_file_name, cmd, ncdf_db, file_name, delta_t1

; get file if remote

   get_file, file_name, ncdf_db

; read grid from file if require (@<grid> option) 

   IF read_grid_from_file EQ 1 THEN BEGIN 
      mesh_from_file, cmd.grid, file_name, ncdf_db, cmd.var
      key_shift_map = key_shift
   ENDIF 

; define data domain
 
   IF hotyp NE '-' THEN BEGIN   ; time serie
      ; get hovmoeller box + stride
      box_plot = def_box(cmd.plt, dimplot, legbox, time_stride)
      IF debug_w THEN print, 'box_plot for domdef in data_read = ', box_plot
      CASE n_elements(box_plot) OF
         4: domdef, box_plot[0:3], /MEMEINDICES
         6: domdef, box_plot[0:5], /MEMEINDICES
      ENDCASE 
      ; time grid for hovmoeller
      time1 = 1+delta_t1
      IF cmd.out EQ 'cdf' THEN nb_cycles = 1
      timearr = compute_time(cmd.timave, cmd.date1, cmd.spec)
      jpt = timearr.count
      time2 = jpt+delta_t1
      time = timearr.scale
      niveau = 1
;      print, '  Data_read, hotyp/jpt = ', hotyp, jpt
;      print, '  Data_read, box = ', box_plot[0:3]
   ENDIF ELSE BEGIN             ; single date
      jpt = 1
      CASE strmid(cmd.timave, strlen(cmd.timave)-1, 1) OF
         'y': time1 = delta_t1+1
         'm': BEGIN
            CASE strmid(cmd.timave, strlen(cmd.timave)-2, 1) OF
               'm': time1 = delta_t1+long(strmid(cmd.date1, 0, 2))  ; mean month
               ELSE: time1 = delta_t1+1
            ENDCASE 
            END 
         'd': time1 = delta_t1+1
         ELSE: time1 = delta_t1+1
      ENDCASE 
      time2 = time1
      time = time1
      print, '    Setting domdef in data_read, plt typ =', plttyp   
      CASE plttyp OF
         'plt': BEGIN
            box_plot = def_box(cmd.plt, 2, legbox, time_stride) 
            IF strmid(cmd.var, 0,2) EQ 'vo' THEN BEGIN ; get right zindex for vertical
               CASE vert_type OF
                  'level':domdef,[ box_plot[0:3],vert_mean],/zindex, /MEMEINDICES
                  'z':domdef,[ box_plot[0:3],vert_mean], /MEMEINDICES
                  ELSE: domdef, box_plot[0:3], /MEMEINDICES
               ENDCASE 
            ENDIF ELSE BEGIN
               domdef, box_plot[0:3], /MEMEINDICES
            ENDELSE
         END 
         'pltz': BEGIN
            box_plot = def_box(cmd.plt, 2, legbox, time_stride) 
            domdef, box_plot[0:3], /MEMEINDICES
         END 
         'plt1d': BEGIN
            box_plot = def_box(cmd.plt, 2, legbox, time_stride) 
            domdef, box_plot[0:3], /MEMEINDICES
         END 
         ELSE: domdef, /MEMEINDICES
      ENDCASE 
   ENDELSE 

;   print, jpt, time1, time2, hotyp

; save time in common fld_att

   time1_r = time1
   time2_r = time2

; call specific read routine

   CASE STRMID(cmd.var, 0, 2) OF
      '@@': BEGIN
         IF debug_w THEN print,  'keyword_set(all_data) : ', keyword_set(all_data)
         fldr = macro_read(file_name, cmd.var, ncdf_db, TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)
         IF stddev_diff EQ 1 THEN fldr.origin =  'diff'
      END 
      ELSE: BEGIN
         CASE plttyp OF
            'yfx': BEGIN        ; scatter plot y=f(x)
               idx = strpos(cmd.var, '=f(')
               var1 = strmid(cmd.var, 0, idx)
               fldr1 = nc_read(file_name, var1, ncdf_db, $
                               TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)
               CASE cmd.grid OF 
                  'U': vargrid1 = 'U'
                  'V': vargrid1 = 'V'
                  ELSE: vargrid1 = 'T'
               ENDCASE 
               IF strpos(cmd.var, '=f(next)') EQ -1 THEN BEGIN  ; type 1 y=f(x) on 1 line
                  var2 = strmid(cmd.var, idx+3, strlen(cmd.var)-idx-4)
                  fldr2 = nc_read(file_name, var2, ncdf_db, $
                                  TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)
                  datyp2 = datyp
                  cmd2 = cmd
                  sw_diffg = 0
                  vargrid2 = vargrid1
               ENDIF ELSE  BEGIN  ; type 2 y=f(next) on 2 lines
                  cmd2 = decode_cmd(cmdline_main, idx_main+1)
                  cmd2.timave = cmd.timave
                  cmd2.date1 = cmd.date1
                  cmd2.spec = cmd.spec
                  datyp2 = def_dptyp(cmd2)
                  ; other grid ?
                  sw_diffg = 0
                  IF cmd2.grid NE cmd.grid THEN BEGIN
                     sw_diffg = 1
                     jptb =  jpt
                     def_grid, cmd2
                     IF debug_w THEN print, 'def_grid, cmd2 dans data_read : '
                     jpt = jptb
                  ENDIF 
                  fldr2 = data_read(cmd2, datyp2.hotyp, datyp2.plttyp, '1', '0', ALL_DATA = all_data, _extra = ex)
                  var2 = cmd2.var
                  print,  'Variable ',  var2,  ' lue sur grille ',  vargrid
                  jptb =  jpt
                  IF sw_diffg EQ 1 THEN BEGIN
                     def_grid, cmd
                     IF debug_w THEN print, 'def_grid, cmd dans data_read : '
                  ENDIF
                  jpt = jptb
                  CASE cmd2.grid OF 
                     'U': vargrid2 = 'U'
                     'V': vargrid2 = 'V'
                     ELSE: vargrid2 = 'T'
                  ENDCASE
               ENDELSE  
                                ; consistency check
               IF fldr1.dim NE fldr2.dim THEN BEGIN
                  print, '  *** dimension inconsistency in y=f(x) : ', cmd.var
                  return, -1
               ENDIF 
                                ; custom stucture for scatter plot
               fldr = {name: var1, data: fldr1.data, legend: fldr1.legend, $
                       units:  fldr1.units, origin: fldr1.origin, dim: fldr1.dim, $
                       name2: var2, data2: fldr2.data, legend2: fldr2.legend, $
                       units2:  fldr2.units, origin2: fldr2.origin, dim2: fldr2.dim}
                                ; diff not possible yet
               IF diff_from_exp EQ 1 OR STRMID(cmd.spec, 0, 2) EQ 'd:' THEN BEGIN
                  diff_from_exp = 0
                  cmd.spec = '-'
                  print, '   data_read : difference in scatter plot not ready'
               ENDIF 
                                ; re-organise cmd.var
               cmd.var = var1
               cmd.var2 = var2
               
            END  
            ELSE: BEGIN
               IF vecplot EQ 1 THEN BEGIN ; vectors case
                  idx = strpos(cmd.var, ',')
                  var1 = strmid(cmd.var, 1, idx-1)
                  var2 = strmid(cmd.var, idx+1, strlen(cmd.var)-idx-2)
;                  print, var1, var2
;old version                  file_nam = strmid(file_name, 0, strlen(file_name)-4)
                  file_nam = strmid(base_file_name+base_suffix, 0, strlen(base_file_name+base_suffix)-1)
                  vargrid = strmid(cmd.grid, 0, 1)
                  CASE STRMID(var1, 0, 2) OF
                     '@@': fldr1 = macro_read(file_name, var1, ncdf_db, TIME_1 = time1,  TIME_2 =  time2, _extra = ex)
                     ELSE: fldr1 = nc_read(file_nam+strmid(cmd.grid, 0, 1)+suff_domain+'.nc', var1, ncdf_db, $
                                             TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)
                  ENDCASE 
                  vargrid = strmid(cmd.grid, 1, 1)
                  CASE STRMID(var2, 0, 2) OF
                     '@@': fldr2 = macro_read(file_name, var2, ncdf_db, TIME_1 = time1,  TIME_2 =  time2, _extra = ex)
                     ELSE: fldr2 = nc_read(file_nam+strmid(cmd.grid, 1, 1)+suff_domain+'.nc', var2, ncdf_db, $
                                  TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)
                  ENDCASE     
                                ; consistency check
                  IF fldr1.dim NE fldr2.dim THEN BEGIN
                     print, '  *** dimension inconsistency in vector plot : ', cmd.var
                     return, -1
                  ENDIF 
                ; custom stucture for vector plot
                  fldr = {name: var1, data: fldr1.data, legend: fldr1.legend, $
                          units:  fldr1.units, origin: fldr1.origin, dim: fldr1.dim, $
                          name2: var2, data2: fldr2.data, legend2: fldr2.legend, $
                          units2:  fldr2.units, origin2: fldr2.origin, dim2: fldr2.dim}
                  IF diff_from_exp EQ 1 OR STRMID(cmd.spec, 0, 2) EQ 'd:' THEN BEGIN
                     diff_from_exp = 0
                     cmd.spec = '-'
                     print, '   data_read : difference in vector plot not ready'
                  ENDIF 
                                ; re-organise cmd.var
                  cmd.var = var1
;                  cmd.grid = strmid(cmd.grid, 2, 1)
                  cmd.var2 = var2
              

               ENDIF ELSE BEGIN  ; general case
                  fldr = nc_read(file_name, cmd.var, ncdf_db, $
                                 TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)

                  ; perform mean if hovmoeller diagram is required
                  ; It is absolutely necessary for the division case with 
                  ; a hovmoeller diagram (before performing the division)
                  IF plttyp EQ 'pltt' AND div_from_exp EQ 1 THEN BEGIN
                     mask_z, fldr.data, cmd, boite_pltt, dimplot, legz
                     z2d = checkfield(fldr.data, plttyp, TYPE = hotyp, BOXZOOM = boite_pltt, DIREC = direc,  _extra = ex)
                     fldr =  { name : fldr.name, data : z2d,  legend : fldr.legend,  units : fldr.units,  $
                               origin : fldr.origin, dim : size(z2d, /N_DIMENSIONS), direc : hotyp}
                     print,  'Averaging made before the call to SAXO plt routines !!!!'
                  ENDIF
; fix
;                  idx = where(fldr.data EQ 1.00000e+20)
;                  IF idx(0) NE -1 THEN fldr.data(idx) = 0.
               ENDELSE 
            END 
         ENDCASE  
      END  
   ENDCASE  

; if tkeavt, take log

;    CASE cmd.var OF 
;       'votkeavt': BEGIN
;          idx = where(fldr.data LE valmask/10)
;          idxm = where(fldr.data EQ valmask)
;          fldr.data(idx) = alog10(fldr.data(idx))
;          ; on T grid
;          vargrid = 'T'
;          fldr.data = (fldr.data+shift(fldr.data, 0, 0, 1))*.5
;          fldr.data(idxm) = valmask
;       END 
;       ELSE: 
;    ENDCASE  


; density projection
   IF splot EQ 1 THEN BEGIN
      IF cmd.timave EQ '1y' AND 1 THEN BEGIN 
         print,  'data_read:  Constructing annual mean from monthly means...'
         cmd.timave = '1m'
         year_s =  cmd.date1
         cumulative = 0.
         FOR i = 1, 12 DO BEGIN
            IF i EQ 1  THEN cmd.date1 = year_s + '01'
            IF i EQ 2  THEN cmd.date1 = year_s + '02'
            IF i EQ 3  THEN cmd.date1 = year_s + '03'
            IF i EQ 4  THEN cmd.date1 = year_s + '04'
            IF i EQ 5  THEN cmd.date1 = year_s + '05'
            IF i EQ 6  THEN cmd.date1 = year_s + '06'
            IF i EQ 7  THEN cmd.date1 = year_s + '07'
            IF i EQ 8  THEN cmd.date1 = year_s + '08'
            IF i EQ 9  THEN cmd.date1 = year_s + '09'
            IF i EQ 10 THEN cmd.date1 = year_s + '10'
            IF i EQ 11 THEN cmd.date1 = year_s + '11'
            IF i EQ 12 THEN cmd.date1 = year_s + '12'
            print,  '   Date1 = ',  cmd.date1
            sfild = fldr
            bin_sigma, cmd, sfild    ;;;;;;;;;;;;;; prob
            cumulative = cumulative + sfild.data
         ENDFOR
         fldr = sfild
         fldr.data = cumulative/12.
         cmd.timave = '1y'
         cmd.date1 = year_s
      ENDIF ELSE BEGIN
         IF cmd.timave EQ '1m' AND strmid(cmd.plt, 0, 2) EQ 'st' AND really_1m_st EQ 1 THEN BEGIN
            sfild = fldr
            timedim = (size(sfild.data))[(size(sfild.data))[0]]
            fldrd = fltarr((size(sfild.data))[1],(size(sfild.data))[2], (sig_max - sig_min)/sig_del + 1, timedim)
            FOR i =  0, timedim-1 DO BEGIN
               print, 'data_read: Performing monthly bining. Indices: first, last, current = 0,',  timedim-1,  i
               sfild2 = {name: sfild.name, data: sfild.data(*, *, *, i), legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim-1}
               bin_sigma, cmd, sfild2
               fldrd(*, *, *, i)= sfild2.data
            ENDFOR 
            fldr = {name: sfild.name, data:fldrd, legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim}
         ENDIF ELSE BEGIN
            sfild = fldr
            bin_sigma, cmd, sfild
            fldr = sfild
         ENDELSE 
      ENDELSE 
   ENDIF 

; modify data info if needed (actual data modification done in; trends.pro called by pltt.pro)

   CASE plttyp OF
      'pltt': BEGIN
         IF run_stddev EQ 0 THEN BEGIN
            CASE strmid(cmd.trend, 0, 1) OF
               '1': fldr.origin = 'diff'
               '2': fldr.origin = 'diff'
               '3': fldr.origin = 'diff'
               '4': fldr.origin = 'diff'
               ELSE:
            ENDCASE
         ENDIF
      END 
      ELSE:
   ENDCASE 
   
   IF n_elements(fldr.data) EQ 1 THEN return, fldr

; Difference 2 cases : from cmd.exp or from cmd.spec


   IF diff_from_exp EQ 1 THEN BEGIN 

      cmd.exp = argvar[1]
      print, '    remove : ', cmd.exp

      ; read field2
      field2 = data_read(cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex)
      IF n_elements(field2.data) EQ 1 THEN return, field2

      ; perform difference

      diff = fldr.data-field2.data
      IF (where(fldr.data EQ valmask))[0] NE -1 THEN $
       diff[where(fldr.data EQ valmask)] = valmask
      fldr.data = diff
      fldr.origin = 'diff'   
      cmd.exp = exp_init
      varexp = cmd.exp 

   ENDIF ELSE IF div_from_exp EQ 1 THEN BEGIN  ; Division : 1 case
      
      cmd.exp = argvar[1]
      print, '    divide by : ', cmd.exp

      div =  fldr.data
      
      ; read field2

      field2 = data_read(cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex)
      IF n_elements(field2.data) EQ 1 THEN return, field2

      IF plttyp EQ 'pltt' THEN BEGIN
         mask_z, field2.data, cmd, boite_pltt, dimplot, legz
         z2d = checkfield(field2.data, plttyp, TYPE = hotyp, BOXZOOM = boite_pltt, DIREC = direc,  _extra = ex)
         field2 =  { name : field2.name, data : z2d,  legend : field2.legend,  units : field2.units,  $
                   origin : field2.origin, dim : size(z2d, /N_DIMENSIONS), direc : hotyp}
         print,  'Averaging made before the call to SAXO plt routines !!!!'
      ENDIF
      
      ; perform division
      IF (where(field2.data EQ 0.0))[0] NE -1 THEN BEGIN
         print,  'Be careful some data are null. They will be masked ! '
         idx0 = where(field2.data EQ 0.0)
         field2.data(idx0) =  valmask
         fldr.data(idx0) =  valmask
         div(idx0) =  valmask
      ENDIF

      idx = where(field2.data NE valmask)
      div(idx) = fldr.data(idx) / field2.data(idx)

      fldr.data = div
      fldr.origin = 'div'   
      cmd.exp = exp_init
      varexp = cmd.exp 
      
   ENDIF ELSE BEGIN

; If difference, decode second field, read and perform difference

      IF STRMID(cmd.spec, 0, 2) EQ 'd:' THEN BEGIN
         
         cmd2 = cmd
         cmdl = STRMID(cmd.spec, 2, strlen(cmd.spec)-2)
         
         print, '    remove : ', cmdl
         
         ; decode field2

         argvar = str_sep(cmdl, '/')
         
         cmd2.exp    = argvar[0]
         cmd2.timave = argvar[1]
         cmd2.date1  = argvar[2]
         cmd2.spec   = '-'
         IF n_elements(argvar) EQ 4 THEN cmd2.spec = argvar[3] 
         
         ; read field2

         field2 = data_read(cmd2, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex)
         IF n_elements(field2.data) EQ 1 THEN return, field2

         ; perform difference

         diff = fldr.data-field2.data
         IF (where(fldr.data EQ valmask))[0] NE -1 THEN $
          diff[where(fldr.data EQ valmask)] = valmask
         fldr.data = diff
         fldr.origin = 'diff'   
         cmd.exp = cmd.exp+' - '+cmd2.exp+ ' ('+cmd2.timave+' '+cmd2.date1+')'
         varexp = cmd.exp 

      ENDIF 

   ENDELSE

ENDIF ELSE BEGIN

   print, ''
   print, '   Data already read (previous field)'
   print, ''
   
ENDELSE

return, fldr

END