Index: /trunk/condmag_cond_sed_on_orca.pro
===================================================================
--- /trunk/condmag_cond_sed_on_orca.pro	(revision 4)
+++ /trunk/condmag_cond_sed_on_orca.pro	(revision 4)
@@ -0,0 +1,346 @@
+;+
+;
+; @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
