source: trunk/condmag_on_orca.pro @ 6

;+
;
; @file_comments
; interpolate condmag.nc file (sediment and magnetic fields) on ORCA grid
; and produce cond_sed_<orcares>.nc and Br_<orcares>.nc
;
; @categories
; interpolation, orca grid
;
; @param orcares {in}{required}{type=string}
; code of ORCA grid to use for output results
; must be 'ORCA2'
;
; @param method {in}{required}{type=string}
; code for interpolation method
; must be 'bilinear' or 'imoms3'
; in fact bilinear is used in this geomag project
; ++ est-ce que ce sera tjs vrai ?
;
; @param gridtype {in}{required}{type=string}
; to specify on which grid we do the interpolation T, U, V
; must belong to T,U,V
;
; @restrictions
;  - Requires SAXO tools
;  - must have condmag.nc in ${GEOMAG_ID}/
;  - must have meshmask in ${GEOMAG_ID}/
;  - must not have cond_sed_*.nc in ${GEOMAG_OD}/
;  - must not have Br_*.nc in ${GEOMAG_OD}/
;
; @todo
; add ORCA025
; provide tools to plot output file
; produce a NetCDF GDT or CF compliant
;
; @pre
; be sure to have datafile condmag.nc in the directory defined in
; ${GEOMAG_ID}/ see geomag_env.sh
; be sure to have meshmask of ORCA grid you choos in the directory defined in
; ${GEOMAG_ID}/
; for ORCA2 filename is meshmask_bab.nc
;
; @post
; cond_sed_<I>orcasres</I>.nc is now present in ${GEOMAG_OD}/
; Br_<I>orcasres</I>.nc is now present in ${GEOMAG_OD}/
; see geomag_env.sh
;
; @examples
; IDL> .run condmag_on_orca
; IDL> condmag_on_orca, 'ORCA2', 'bilinear','T'
;
; @history
; fplod 2006-11-23T08:57:10Z aedon.locean-ipsl.upmc.fr (Darwin)
; rename to condmag_on_orca.pro
; fplod 2006-11-22T15:10:31Z aedon.locean-ipsl.upmc.fr (Darwin)
; add creation of Br_<I>orcasres</I>.nc
; ++ not very beautifuly implemented
; fplod 2006-11-22T10:38:51Z aedon.locean-ipsl.upmc.fr (Darwin)
; always use msg = 'iii : ...' or msg = 'eee : ...' 
; information on open io
; fplod 2006-11-20T10:32:02Z aedon.locean-ipsl.upmc.fr (Darwin)
; cleaning before introduction into svn repository
; fplod 2006-03-23T13:05:39Z aedon.lodyc.jussieu.fr (Darwin)
; cond_sed_ORCA2 presque ok manque attributes min/max de cond_sed
; fplod 2006-03-22T09:35:42Z aedon.lodyc.jussieu.fr (Darwin)
; created from
; /Users/fplod/incas/seb/DORAEMON/wind/idl/quikscat_to_orca_scalar.pro written
; by Sebastien Masson
; to reproduce  /usr/work/sur/fvi/OPA/geomag/cond_sed_ORCA2.nc
; main differences : no yyyy parameter, no time loop , no mask and no missing
; values (++ to be checked ) in data input (condmag.nc), no scale factor,
; no OFFSET, no save of weight and addresses
;
; @version
; $Id$
;
;-
PRO condmag_on_orca, orcares, method, gridtype
;
  compile_opt idl2, strictarrsubs
;
;----
; check input parameters
;----
;
; check orcares definition
;
  CASE orcares OF
     'ORCA2': BEGIN
                 msg = 'iii : valid orcares parameter = ' +  orcares
                 PRINT, msg
                 filename_oce = 'meshmask_bab.nc'
              END
      ELSE  : BEGIN
                 msg = 'eee : invalid orcares parameter = ' + orcares
                 PRINT, msg
                 msg = 'eee : orcares must be ORCA2'
                 PRINT, msg
                 RETURN
              END
  ENDCASE
;
; check method definition
  CASE method OF
     'bilinear': BEGIN
                    msg = 'iii : valid method parameter = ' + method
                    PRINT, msg
                 END
     'imoms3'  : BEGIN
                    msg = 'iii : valid method parameter = ' + method
                    PRINT, msg
                 END
     ELSE      : BEGIN
                    msg = 'eee : invalid method parameter = ' + method 
                    PRINT, msg
                    msg = 'eee : method must be bilinear or imoms3'
                    PRINT, msg
                    RETURN
                 END
  ENDCASE
;
; check gridtype definition
  gridtype = strupcase(gridtype)
  CASE gridtype OF
     'T' : BEGIN
              msg = 'iii : valid gridtype parameter = ' + gridtype
              PRINT, msg
           END
     'U' : BEGIN
              msg = 'iii : valid gridtype parameter = ' + gridtype
              PRINT, msg
           END
     'V' : BEGIN
              msg = 'iii : valid gridtype parameter = ' + gridtype
              PRINT, msg
           END
     ELSE: BEGIN
              msg = 'eee : invalid gridtype parameter = ' + gridtype 
              PRINT, msg
              msg = 'eee : gridtype must be T, U or V'
              PRINT, msg
              RETURN
           END
  ENDCASE
;
; check for input files
;
; test if ${GEOMAG_ID} defined
  geomag_id_env=GETENV('GEOMAG_ID')
  CASE geomag_id_env OF
     ''  :  BEGIN
              msg = 'eee : ${GEOMAG_ID} is not defined'
              PRINT, msg
              RETURN
            END
     ELSE: BEGIN
             msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env
             PRINT, msg
           END 
  ENDCASE
;
  iodirin = isadirectory(geomag_id_env)
;
; existence and protection of ${GEOMAG_ID}
  IF (FILE_TEST(iodirin, /DIRECTORY,/EXECUTABLE , /READ) EQ 0) THEN BEGIN
     msg = 'eee : the directory' + iodirin  + ' is not accessible.'
     PRINT, msg
     RETURN
  ENDIF
;
; conductivity and magnetic field
  filename_condmag = 'condmag.nc'
;
; check if this file exists
  fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0,/MUST_EXIST)
  IF fullfilename_condmag[0] EQ '' THEN BEGIN
     msg = 'eee : the file ' + fullfilename_condmag + ' was not found.'
     PRINT, msg
     RETURN
  ENDIF
;
; protection
  IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN
     msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.'
     PRINT, msg
     RETURN
  ENDIF
;
; mesh mask
; check if this file exists
  fullfilename_oce = isafile(iodirin + filename_oce, NEW=0,/MUST_EXIST, RECURSIVE=0)
  IF fullfilename_oce[0] EQ '' THEN BEGIN
     msg = 'eee : the file ' + fullfilename_oce + ' was not found.'
     PRINT, msg
     RETURN
  ENDIF
;
; protection
  IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN
     msg = 'eee : the file ' + fullfileoce_condmag[0] + ' is not readable.'
     PRINT, msg
     RETURN
  ENDIF
;
; test if ${GEOMAG_OD} defined
  geomag_od_env=GETENV('GEOMAG_OD')
  CASE geomag_od_env OF
     '' : BEGIN
             msg = 'eee : ${GEOMAG_OD} is not defined'
             PRINT, msg
             RETURN
          END
     ELSE: BEGIN
             msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env
             PRINT, msg
           END 
  ENDCASE
;
; check if output data will be possible
  iodirout = isadirectory(geomag_od_env)
;
; existence and protection
  IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN
     msg = 'eee : the directory' + iodirout  + ' was not found.'
     PRINT, msg
     RETURN
  ENDIF
;
; build output filenames
  filename_cond_sed = 'cond_sed' + '_' + orcares +'.nc'
  fullfilename_cond_sed = iodirout + filename_cond_sed
;
; in order to avoid unexpected overwritten
  IF (FILE_TEST(fullfilename_cond_sed) EQ 1) THEN BEGIN
     msg = 'eee : the file ' + fullfilename_cond_sed  + ' already exists.'
     PRINT, msg
     RETURN
  ENDIF
;
  filename_br = 'Br' + '_' + orcares +'.nc'
  fullfilename_br = iodirout + filename_br
;
; in order to avoid unexpected overwritten
  IF (FILE_TEST(fullfilename_br) EQ 1) THEN BEGIN
     msg = 'eee : the file ' + fullfilename_br  + ' already exists.'
     PRINT, msg
     RETURN
  ENDIF
;
; d'après ncdump -h /usr/work/sur/fvi/OPA/geomag/condmag.nc
  condmaglonname = 'lo'
  condmaglatname = 'la'
  varname_cond_sed = 'cond_sed'
  varname_br = 'Br'
;
;----
; conductivity and magnetic field grid parameters
;----
;
  get_gridparams, fullfilename_condmag[0], $
     condmaglonname, condmaglatname, $
     condmaglon, condmaglat, jpia, jpja, 1,/DOUBLE
;
;----
; Oceanic grid parameters
;----
;
  olonname = 'glam' + STRLOWCASE(gridtype)
  olatname = 'gphi' + STRLOWCASE(gridtype)
  get_gridparams, fullfilename_oce[0], $
     olonname, olatname, $
     olon, olat, jpio, jpjo, 2,/DOUBLE
;
;---------------
; Compute weight and address
;---------------
;
  CASE method OF
     'bilinear': compute_fromreg_bilinear_weigaddr, $
                    condmaglon, condmaglat, olon, olat, weig, addr
     'imoms3'  : compute_fromreg_imoms3_weigaddr, $
                    condmaglon, condmaglat, olon, olat, weig, addr
  ENDCASE
;
; reading condmag file
  netcdf_id_condmag = NCDF_OPEN(fullfilename_condmag[0], /NOWRITE)
  msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read'
  PRINT, msg
  
  varinq_cond_sed = NCDF_VARINQ(netcdf_id_condmag, varname_cond_sed)
  varinq_br = NCDF_VARINQ(netcdf_id_condmag, varname_br)
;
;---------------------------
; Creation of the NetCdf file for cond_sed and Br
;---------------------------
;
  netcdf_id_cond_sed = NCDF_CREATE(fullfilename_cond_sed, /clobber)
  NCDF_CONTROL, netcdf_id_cond_sed, /NOFILL
  netcdf_id_br = NCDF_CREATE(fullfilename_br, /clobber)
  NCDF_CONTROL, netcdf_id_br, /NOFILL
;
; dimension
  dimidx = NCDF_DIMDEF(netcdf_id_cond_sed, 'x' ,  jpio)
  dimidy = NCDF_DIMDEF(netcdf_id_cond_sed, 'y' ,  jpjo)
  dimidt = NCDF_DIMDEF(netcdf_id_cond_sed, 'lo', /UNLIMITED)
  dimidx = NCDF_DIMDEF(netcdf_id_br, 'x' ,  jpio)
  dimidy = NCDF_DIMDEF(netcdf_id_br, 'y' ,  jpjo)
  dimidt = NCDF_DIMDEF(netcdf_id_br, 'lo', /UNLIMITED)
;
; global attributes
  NCDF_ATTPUT, netcdf_id_cond_sed, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformité
  NCDF_ATTPUT, netcdf_id_cond_sed, 'file_name'  , filename_cond_sed, /GLOBAL
  NCDF_ATTPUT, netcdf_id_cond_sed, 'Title'      , 'Conductance', /GLOBAL
  NCDF_ATTPUT, netcdf_id_br, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformité
  NCDF_ATTPUT, netcdf_id_br, 'file_name'  , filename_Br, /GLOBAL
  NCDF_ATTPUT, netcdf_id_br, 'Title'      , 'Magnetic field', /GLOBAL
;
; declaration of variables
; 4 common variables for the two files to produce
  varid = lonarr(3)
;
  varid[0] = NCDF_VARDEF(netcdf_id_cond_sed, 'nav_lon'  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[0], 'units'    , 'degrees_east'
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[0], 'valid_min', min(olon, max = omax),/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[0], 'valid_max', omax,/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[0], 'long_name', 'Longitude at t-point'
  varid[0] = NCDF_VARDEF(netcdf_id_br, 'nav_lon'  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id_br, varid[0], 'units'    , 'degrees_east'
  NCDF_ATTPUT, netcdf_id_br, varid[0], 'valid_min', min(olon, max = omax),/FLOAT
  NCDF_ATTPUT, netcdf_id_br, varid[0], 'valid_max', omax,/FLOAT
  NCDF_ATTPUT, netcdf_id_br, varid[0], 'long_name', 'Longitude at t-point'
;
  varid[1] = NCDF_VARDEF(netcdf_id_cond_sed, 'nav_lat'  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[1], 'units'    , 'degrees_north'
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[1], 'valid_min', min(olat, max = omax),/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[1], 'valid_max', omax,/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid[1], 'long_name', 'Latitude at t-point'
  varid[1] = NCDF_VARDEF(netcdf_id_br, 'nav_lat'  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id_br, varid[1], 'units'    , 'degrees_north'
  NCDF_ATTPUT, netcdf_id_br, varid[1], 'valid_min', min(olat, max = omax),/FLOAT
  NCDF_ATTPUT, netcdf_id_br, varid[1], 'valid_max', omax,/FLOAT
  NCDF_ATTPUT, netcdf_id_br, varid[1], 'long_name', 'Latitude at t-point'
;
  varid[2] = NCDF_VARDEF(netcdf_id_cond_sed, 'time'     , [dimidt], /FLOAT)
  varid[2] = NCDF_VARDEF(netcdf_id_br, 'time'     , [dimidt], /FLOAT)
;
; each of the two files to produce contains a specific variable
;
  varid_cond_sed = lonarr(1)
  varid_cond_sed[0] = NCDF_VARDEF(netcdf_id_cond_sed, varname_cond_sed, [dimidx, dimidy, dimidt], /FLOAT)
;
; variable 3
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'long_name', 'Conductance' ; ++ non cf 
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'units', 'siemens' ; ++ récupérer l'unite de cond_sed dans condmag.nc
; pour min  et max, il faut avoir lu la variable ... cf. plus bas ++
; donc pour l'instant on les met à valeur manquante
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_min', !VALUES.F_NAN ,/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_max', !VALUES.F_NAN ,/FLOAT
;
  varid_br = lonarr(1)
  varid_br[0] = NCDF_VARDEF(netcdf_id_br, varname_br, [dimidx, dimidy, dimidt], /FLOAT)
;
; variable 3
  NCDF_ATTPUT, netcdf_id_br, varid_br[0], 'long_name', 'magnetic field' ; ++ non cf 
  NCDF_ATTPUT, netcdf_id_br, varid_br[0], 'units', 'tesla' ; ++ récupérer l'unité de Br dans condmag.nc
; pour min  et max, il faut avoir lu la variable ... cf. plus bas ++
; donc pour l'instant on les met à valeur manquante
  NCDF_ATTPUT, netcdf_id_br, varid_br[0], 'valid_min', !VALUES.F_NAN ,/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_max', !VALUES.F_NAN ,/FLOAT
;
;---------------------------
; end of header definition, writing of the NetCdf files
;---------------------------
  NCDF_CONTROL, netcdf_id_cond_sed, /ENDEF
  NCDF_CONTROL, netcdf_id_br, /ENDEF
  NCDF_CLOSE, netcdf_id_condmag ; ++ si le close n'est pas ici ça merde mais je ne sais pas pourquoi
;
; grid
  NCDF_VARPUT, netcdf_id_cond_sed, 'nav_lon', olon
  NCDF_VARPUT, netcdf_id_cond_sed, 'nav_lat', olat
  NCDF_VARPUT, netcdf_id_cond_sed, varid[2], FLOAT(0.5) ; ++ c'est quoi cette valeur
  NCDF_VARPUT, netcdf_id_br, 'nav_lon', olon
  NCDF_VARPUT, netcdf_id_br, 'nav_lat', olat
  NCDF_VARPUT, netcdf_id_br, varid[2], FLOAT(0.5) ; ++ c'est quoi cette valeur

;---------------------------
; réouverture du fichier ... mais pourquoi on le relit !!++
  netcdf_id_condmag = NCDF_OPEN(fullfilename_condmag[0],/NOWRITE)
  msg = 'iii : ' + fullfilename_condmag[0] + ' reopened for read'
  PRINT, msg
;
  NCDF_VARGET, netcdf_id_condmag, varname_cond_sed, varin_cond_sed ; , COUNT = [jpia, jpja, 1], OFFSET = [0, 0, 0]
  NCDF_VARGET, netcdf_id_condmag, varname_br, varin_br ; , COUNT = [jpia, jpja, 1], OFFSET = [0, 0, 0]
;
; do the interpolation
      varin_cond_sed = TOTAL(weig*varin_cond_sed[addr], 1)
      varin_cond_sed = REFORM(varin_cond_sed, jpio, jpjo, /OVER)
      varin_br = TOTAL(weig*varin_br[addr], 1)
      varin_br = REFORM(varin_br, jpio, jpjo, /OVER)
;
      NCDF_CLOSE, netcdf_id_condmag
      varout_cond_sed = TEMPORARY(varin_cond_sed)
      varout_br = TEMPORARY(varin_br)
; compute min max
      minarr_cond_sed = min(varout_cond_sed, max = maxarr_cond_sed)
      minarr_br = min(varout_br, max = maxarr_br)
;
; put back the masked value
;++      IF bad[0] NE -1 THEN varout_cond_sed[TEMPORARY(bad)] = 32767
;++      IF bad[0] NE -1 THEN varout_br[TEMPORARY(bad)] = 32767
;
; write the data
      NCDF_VARPUT, netcdf_id_cond_sed, varname_cond_sed, varout_cond_sed, COUNT = [jpio, jpjo, 1], OFFSET = [0, 0, 0]
      NCDF_VARPUT, netcdf_id_br, varname_br, varout_br, COUNT = [jpio, jpjo, 1], OFFSET = [0, 0, 0]
;
; update min max attributes
  NCDF_CONTROL, netcdf_id_cond_sed, /REDEF
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_min', minarr_cond_sed,/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_max', maxarr_cond_sed,/FLOAT
  NCDF_CONTROL, netcdf_id_cond_sed, /ENDEF
  NCDF_CONTROL, netcdf_id_br, /REDEF
  NCDF_ATTPUT, netcdf_id_br, varid_br[0], 'valid_min', minarr_br,/FLOAT
  NCDF_ATTPUT, netcdf_id_br, varid_br[0], 'valid_max', maxarr_br,/FLOAT
  NCDF_CONTROL, netcdf_id_br, /ENDEF
;
;---------------------------
; close the netcdf files
  NCDF_CLOSE, netcdf_id_cond_sed
  NCDF_CLOSE, netcdf_id_br
;
  msg = 'iii : ' + fullfilename_cond_sed + ' created'
  PRINT, msg
  msg = 'iii : ' + fullfilename_br + ' created'
  PRINT, msg
;
END