;+
;
; @file_comments
; compute delta between *_5d_yyyy_grid_T_orcares.pro Netcdf files with same dimension
;
; @param file1 {in}{required} {type=string}
; name of the first file to be read
;
; @param file2 {in}{required} {type=string}
; name of the second file to be read
;
; @param file3 {in}{required} {type=string}
; file where differences between variables are written
;
; @examples
; test (for header and delta)
; to compute difference between Sigma in Sigma_5d_1993_grid_T_ORCA2.nc and
; itself and write it in ginette.nc :
; $ cd ${GEOMAG}
; $ idl
; IDL> file1=getenv('GEOMAG_OD') + '/Sigma_5d_1993_grid_T_ORCA2.nc'
; IDL> file2=getenv('GEOMAG_OD') + '/Sigma_5d_1993_grid_T_ORCA2.nc'
; IDL> step2_diff, file1, file2, 'ginette.nc'
; values of delta must be 0 everywhere
; $ rm ginette.nc
;
; real life
; to compute difference between Sigma in
; /usr/work/sur/fvi/OPA/ORCA2/Sigma_5d_1994_grid_T.nc and
; ${GEOMGA_OD}/Sigma_5d_1994_grid_T_ORCA2.nc and write it in
; Sigma_5d_1994_grid_T_ORCA2_diff.nc :
; $ cd ${GEOMAG}
; $ idl
; IDL> year=1994
; IDL> file1= '/usr/work/sur/fvi/OPA/ORCA2/Sigma_5d_1994_grid_T.nc'
; IDL> file2=getenv('GEOMAG_OD') + '/Sigma_5d_1994_grid_T_ORCA2.nc'
; IDL> step2_diff, file1, file2, 'Sigma_5d_1994_grid_T_ORCA2_diff.nc'
;
; to see the difference, for example
; IDL> xxx,'Sigma_5d_1994_grid_T_ORCA2_diff.nc'
; select xy on plt wigdet
;
; see also step2_diff.sh for quick start
;
; @history
; fplod 2007-08-22T09:39:45Z aedon.locean-ipsl.upmc.fr (Darwin)
; correction
; fplod 2007-07-30T10:58:12Z aedon.locean-ipsl.upmc.fr (Darwin)
; like <progeomag>step1_diff</progeomag>
;
; @todo
; check difference between headers
;
; @version
; $Id$
;
PRO step2_diff, file1, file2, file3
;
ncverbose=1
;
; check if file1 exists
fullfile1 = isafile(file1, NEW=0,/MUST_EXIST)
IF fullfile1[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfile1 + ' was not found.'
   ras = report(msg)
   RETURN
ENDIF
;
; protection
IF (FILE_TEST(fullfile1[0], /READ) EQ 0) THEN BEGIN
   msg = 'eee : the file ' + fullfile1[0] + ' is not readable.'
   ras = report(msg)
   RETURN
ENDIF
;
; check if file2 exists
fullfile2 = isafile(file2, NEW=0,/MUST_EXIST)
IF fullfile2[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfile2 + ' was not found.'
   ras = report(msg)
   RETURN
ENDIF
;
; protection
IF (FILE_TEST(fullfile2[0], /READ) EQ 0) THEN BEGIN
   msg = 'eee : the file ' + fullfile2[0] + ' is not readable.'
   ras = report(msg)
   RETURN
ENDIF
;
; read the first file
ncdf_read,  fullfile1[0], file1info, file1dinfo, file1vinfo, $
            file1gatts, file1vatts, file1data
;
; read the second file
ncdf_read,  fullfile2[0], file2info, file2dinfo, file2vinfo, $
            file2gatts, file2vatts, file2data
;
; compute delta
delta = *file1data[4] - *file2data[4]
minarr = min(delta, max = maxarr)
;
; check for number of non zero in delta
checknonzero=where(delta NE 0.,count)
IF count EQ 0 THEN BEGIN
   msg = 'iii : delta is zero everywhere'
   ras = report(msg)
ENDIF ELSE BEGIN
   msg = 'iii : delta is not zero ' + STRING(count)  + ' times'
   ras = report(msg)
ENDELSE
;
; write output
; in order to avoid unexpected overwritten
IF (FILE_TEST(file3) EQ 1) THEN BEGIN
   msg = 'eee : the file ' + file3  + ' already exists.'
   ras = report(msg)
   RETURN
ENDIF
;
;---------------------------
; Creation of the NetCdf output file
;---------------------------
;
  netcdf_id = NCDF_CREATE(file3, /clobber)
  NCDF_CONTROL, netcdf_id, /NOFILL
;
; dimension
  dimidx = NCDF_DIMDEF(netcdf_id, file1dinfo[0].name ,  file1dinfo[0].size)
  dimidy = NCDF_DIMDEF(netcdf_id, file1dinfo[1].name ,  file1dinfo[1].size)
  dimidz = NCDF_DIMDEF(netcdf_id, file1dinfo[2].name ,  file1dinfo[2].size)
  dimidt = NCDF_DIMDEF(netcdf_id, file1dinfo[3].name, /UNLIMITED)
  dimidz2 = NCDF_DIMDEF(netcdf_id, file1dinfo[4].name ,  file1dinfo[4].size)
;
; global attributes
  NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite
  NCDF_ATTPUT, netcdf_id, 'file_name'  , file3, /GLOBAL
  NCDF_ATTPUT, netcdf_id, 'Title'      , $
               'delta between '+ file1 + ' and ' + file2, /GLOBAL
;
; declaration of variables
  varid = lonarr(5)
;
  varid[0] = NCDF_VARDEF(netcdf_id, file1vinfo[0].name  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[0].name , $
               STRING(*(*file1vatts[0])[0].values)
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[1].name, $
               STRING(*(*file1vatts[0])[1].values)
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[2].name, $
               STRING(*(*file1vatts[0])[2].values)
;
  varid[1] = NCDF_VARDEF(netcdf_id, file1vinfo[1].name  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[0].name, $
               STRING(*(*file1vatts[1])[0].values)
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[1].name, $
               STRING(*(*file1vatts[1])[1].values)
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[2].name, $
               STRING(*(*file1vatts[1])[2].values)
;
  varid[2] = NCDF_VARDEF(netcdf_id, file1vinfo[2].name, [dimidz2], /FLOAT);
  NCDF_ATTPUT, netcdf_id, varid[2], $
               (*file1vatts[2])[0].name, $
               STRING(*(*file1vatts[2])[0].values)
  NCDF_ATTPUT, netcdf_id, varid[2], $
               (*file1vatts[2])[1].name, $
               STRING(*(*file1vatts[2])[1].values)
  NCDF_ATTPUT, netcdf_id, varid[2], $
               (*file1vatts[2])[2].name, $
               STRING(*(*file1vatts[2])[2].values)
;
  varid[3] = NCDF_VARDEF(netcdf_id, file1vinfo[3].name, [dimidt], /FLOAT);
  NCDF_ATTPUT, netcdf_id, varid[3], $
               (*file1vatts[3])[0].name, $
               STRING(*(*file1vatts[3])[0].values)
  NCDF_ATTPUT, netcdf_id, varid[3], $
               (*file1vatts[3])[1].name, $
               STRING(*(*file1vatts[3])[1].values)
  NCDF_ATTPUT, netcdf_id, varid[3], $
               (*file1vatts[3])[2].name, $
               STRING(*(*file1vatts[3])[2].values)
  NCDF_ATTPUT, netcdf_id, varid[3], $
               (*file1vatts[3])[3].name, $
               STRING(*(*file1vatts[3])[3].values)
  NCDF_ATTPUT, netcdf_id, varid[3], $
               (*file1vatts[3])[4].name, $
               STRING(*(*file1vatts[3])[4].values)
help,delta

  varid[4] =  NCDF_VARDEF(netcdf_id, 'delta',   [dimidx, dimidy,dimidz, dimidt], /DOUBLE)
  NCDF_ATTPUT, netcdf_id, varid[4], 'long_name', 'delta'
  NCDF_ATTPUT, netcdf_id, varid[4], 'valid_min', minarr ,/DOUBLE
  NCDF_ATTPUT, netcdf_id, varid[4], 'valid_max', maxarr ,/DOUBLE
;
;---------------------------
; end of header definition, writing of the NetCdf files
;---------------------------
;
  NCDF_CONTROL, netcdf_id, /ENDEF
;
  NCDF_VARPUT, netcdf_id, file1vinfo[0].name, (*file1data[0])
  NCDF_VARPUT, netcdf_id, file1vinfo[1].name, (*file1data[1])
  NCDF_VARPUT, netcdf_id, file1vinfo[2].name, (*file1data[2])
  NCDF_VARPUT, netcdf_id, file1vinfo[3].name, (*file1data[3])
;
; write the data
  NCDF_VARPUT, netcdf_id, 'delta', delta
;
;---------------------------
; close the netcdf files
  NCDF_CLOSE, netcdf_id
;
  msg = 'iii : ' + file3 + ' created'
  ras = report(msg)
;
END