source: trunk/condmag_cond_sed_on_orca.pro @ 5

;+
;
; @file_comments
; interpolate condmag.nc file (sediment and magnetic fields) on ORCA grid
; and produce one file cond_sed_<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}/
;
; @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}/
; see geomag_env.sh
;
; @examples
; IDL> .run condmag_cond_sed_on_orca
; IDL> condmag_cond_sed_on_orca, 'ORCA2', 'bilinear','T'
;
; @history
; 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_cond_sed_on_orca, orcares, method, gridtype
;
  compile_opt idl2, strictarrsubs
;
;----
; check input parameters
;----
;
; check orcares definition
;
  CASE orcares OF
     'ORCA2': BEGIN
                 filename_oce = 'meshmask_bab.nc'
                 PRINT, 'valid orcares parameter'
              END
      ELSE  : BEGIN
                 msg = 'orcares must be ORCA2'
                 PRINT, msg
                 RETURN
              END
  ENDCASE
;
; check method definition
  CASE method OF
     'bilinear': BEGIN
                    PRINT, 'valid method parameter'
                 END
     'imoms3'  : BEGIN
                    PRINT, 'valid method parameter'
                 END
     ELSE      : BEGIN
                    msg = 'method must be bilinear or imoms3'
                    PRINT, msg
                    RETURN
                 END
  ENDCASE
;
; check gridtype definition
  gridtype = strupcase(gridtype)
  CASE gridtype OF
     'T' : BEGIN
              PRINT, 'valid gridtype parameter'
           END
     'U' : BEGIN
              PRINT, 'valid gridtype parameter'
           END
     'V' : BEGIN
              PRINT, 'valid gridtype parameter'
           END
     ELSE: BEGIN
              msg = '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 = '${GEOMAG_ID} is not defined'
              PRINT, msg
              RETURN
            END
     ELSE: BEGIN
             msg = '${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 = '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)
  IF fullfilename_condmag[0] EQ '' THEN BEGIN
     msg = 'the file ' + fullfilename_condmag + ' was not found.'
     PRINT, msg
     RETURN
  ENDIF
;
; protection
  IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN
     msg = '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)
  IF fullfilename_oce[0] EQ '' THEN BEGIN
     msg = 'the file ' + fullfilename_oce + ' was not found.'
     PRINT, msg
     RETURN
  ENDIF
;
; protection
  IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN
     msg = '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 = '${GEOMAG_OD} is not defined'
             PRINT, msg
             RETURN
          END
     ELSE: BEGIN
             msg = '${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 = 'the directory' + iodirout  + ' was not found.'
     PRINT, msg
     RETURN
  ENDIF
;
; build output filename
  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 = 'the file ' + fullfilename_cond_sed  + ' 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'
;
;----
; 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])
  varinq_cond_sed = NCDF_VARINQ(netcdf_id_condmag, varname_cond_sed)
;
;---------------------------
; Creation of the NetCdf file for cond_sed
;---------------------------
;
  netcdf_id_cond_sed = NCDF_CREATE(fullfilename_cond_sed, /clobber)
  NCDF_CONTROL, netcdf_id_cond_sed, /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)
;
; global attributes
  NCDF_ATTPUT, netcdf_id_cond_sed, 'Conventions', 'GDT 1.2', /GLOBAL
  NCDF_ATTPUT, netcdf_id_cond_sed, 'file_name'  , filename_cond_sed, /GLOBAL
  NCDF_ATTPUT, netcdf_id_cond_sed, 'Title'      , 'Monthly Levitus Sea Salinity corrected', /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[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[2] = NCDF_VARDEF(netcdf_id_cond_sed, 'time'     , [dimidt], /FLOAT)
;
; each of the 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', 'Salinity' ; ++ faux et non cf
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'units', 'siemens'
; 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
;---------------------------
; end of header definition, writing of the NetCdf file
;---------------------------
  NCDF_CONTROL, netcdf_id_cond_sed, /ENDEF
  NCDF_CLOSE, netcdf_id_condmag ; ++ si le close n'est pa ici ça merde mais je ne sais pas pourquoi
; grid
;---------------------------
  NCDF_VARPUT, netcdf_id_cond_sed, 'nav_lon', TEMPORARY(olon)
  NCDF_VARPUT, netcdf_id_cond_sed, 'nav_lat', TEMPORARY(olat)
  NCDF_VARPUT, netcdf_id_cond_sed, varid[2], FLOAT(0.5) ; ++ c'est quoi cette valeur
;---------------------------
; data
;---------------------------
; réouverture du fichier ... mais pourquoi on le relit !!++
      netcdf_id_condmag = NCDF_OPEN(fullfilename_condmag[0])
;
      NCDF_VARGET, netcdf_id_condmag, varname_cond_sed, varin_cond_sed ; , 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)
;
      NCDF_CLOSE, netcdf_id_condmag
      varout_cond_sed = TEMPORARY(varin_cond_sed)
; compute min max
      minarr = min(varout_cond_sed, max = maxarr)
;
; put back the masked value
;++      IF bad[0] NE -1 THEN varout_cond_sed[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]
;
; update min max attributes
  NCDF_CONTROL, netcdf_id_cond_sed, /REDEF
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_min', minarr,/FLOAT
  NCDF_ATTPUT, netcdf_id_cond_sed, varid_cond_sed[0], 'valid_max', maxarr,/FLOAT
  NCDF_CONTROL, netcdf_id_cond_sed, /ENDEF
;
;---------------------------
; close the netcdf file
  NCDF_CLOSE, netcdf_id_cond_sed
  msg =' done'
  PRINT, msg
;
END