source: trunk/condmag_on_orca.pro @ 14

;+
;
; @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' or 'ORCA025'
;
; @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 or V
;
; @keyword DRAKKAR_EXP {type=string}
; code for Drakkar experiment
; only used when orcares = ORC025
; must be G42 ++ G70
;
; @restrictions
;  - Requires SAXO tools
;  - must have condmag.nc in ${GEOMAG_ID}/
;  - must have meshmask for ORCA2 or ORCA025 in ${GEOMAG_ID}/
;  - must not have cond_sed_*.nc in ${GEOMAG_OD}/
;  - must not have Br_*.nc in ${GEOMAG_OD}/
;
; @todo
; provide tools to plot output files
; 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 choose in the directory defined in
; ${GEOMAG_ID}/
; for ORCA2 filename is meshmask_bab.nc
; ++ au pif entre mesh_hgr, mesh_zgr et mask
; for ORCA0252 filename is ORCA025-G42_mesh_hgr.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
; reee522 2006-12-13T13:50:32Z rhodes (IRIX64)
; add optional keyword drakkar_exp
; reee522 2006-12-13T12:09:05Z rhodes (IRIX64)
; externalisation of netcdf writing
; reee522 2006-11-23T13:34:31Z rhodes (IRIX64)
; add ORCA025 (beginning)
; reee522 2006-11-23T13:28:07Z rhodes (IRIX64)
; information about the opening of meshmask file was missing ... because
; the open is hidden in the call of get_gridparam
; 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, DRAKKAR_EXP = drakkar_exp
;
  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'
                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
                    msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp
                    PRINT, msg
                 END
              END
     'ORCA025': BEGIN
                 msg = 'iii : valid orcares parameter = ' + orcares
                 PRINT, msg
                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
                    msg = 'iii : DRAKKAR_EXP keyword set'
                    PRINT, msg
                    msg = 'iii : DRAKKAR_EXP = ' + drakkar_exp
                    PRINT, msg
                    CASE drakkar_exp OF
                       'G42' : BEGIN
                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
                                  PRINT, msg
                               END
                       'G70' : BEGIN
                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
                                  PRINT, msg
                               END
                       ELSE : BEGIN
                                  msg = 'eee : invalid DRAKKAR_EXP keyword = ' + drakkar_exp
                                  PRINT, msg
                                  RETURN
                              END
                    ENDCASE
                    filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc'
                 END
              END
      ELSE  : BEGIN
                 msg = 'eee : invalid orcares parameter = ' + orcares
                 PRINT, msg
                 msg = 'eee : orcares must be ORCA2 or ORCA025'
                 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
;
; d'après ncdump -h /usr/work/sur/fvi/OPA/geomag/condmag.nc
  condmaglonname = 'lo'
  condmaglatname = 'la'
;
;----
; 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
  msg = 'iii : ' + fullfilename_oce[0] + ' opened for read'
  PRINT, msg
;
;---------------
; 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
;
  varname_cond_sed = 'cond_sed'
  varname_br = 'Br'
;
  varinq_cond_sed = NCDF_VARINQ(netcdf_id_condmag, varname_cond_sed)
  NCDF_VARGET, netcdf_id_condmag, varname_cond_sed, varin_cond_sed
;
  varinq_br = NCDF_VARINQ(netcdf_id_condmag, varname_br)
  NCDF_VARGET, netcdf_id_condmag, varname_br, varin_br
;
  NCDF_CLOSE, netcdf_id_condmag
;
;---------------------------
; 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)
;
      varout_cond_sed = TEMPORARY(varin_cond_sed)
      varout_br = TEMPORARY(varin_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
;
;
; produce outputs
  condmag_output, orcares, $
     varinq_cond_sed, 'Conductance', 'Conductance', 'siemens', $
     jpio, jpjo, olon, olat, $
     varout_cond_sed
  condmag_output, orcares, $
     varinq_br, 'Magnetic field', 'magnetic field', 'tesla', $
     jpio, jpjo, olon, olat, $
     varout_br
;
END