source: trunk/initorca2_bab.pro @ 20

; +
;
; @file_comments
; initialisation on ORCA grid
;
; @param orcares {in}{required}{type=string}
; code of ORCA grid to use for output results
; must be 'ORCA2' or 'ORCA025'
;
; @keyword DRAKKAR_EXP {type=string}
; code for Drakkar experiment
; only used when orcares = ORCA025
; must be G42 ++ G70
;
; @keyword EXTRA
; must exist in order to be able to use this procedure as an argument
;
; @examples
; appel classique
; IDL> initocemesh,'ORCA2'
;
; to plot ${GEOMAG_OD}/cond_sed_ORCA2.nc :
; IDL> myfile=getenv('GEOMAG_OD') + 'cond_sed_ORCA2.nc'
; IDL> xxx, myfile, 'initocemesh', '''ORCA2''',format='(i5)',
;
; @history
; F. Vivier 2003
; F. Pinsard adaptation to new SAXO
; according to SAXO ./SRC/Tests/tst_initorca2.pro
; fplod 2007-05-15T13:43:21Z aedon.locean-ipsl.upmc.fr (Darwin)
; transformation in procedure with parameter orcares
; fplod 2007-06-05T12:22:38Z aedon.locean-ipsl.upmc.fr (Darwin)
; (following Sebastien Masson advices)
; . add EXTRA keyword in order to be able to use this procedure as an argument
; of <pro>xxx</pro>.
; . add  EXTRA to the call of ncdf_meshlec
; replace initorca2_bab by initocemesh in filename and procedure name
;
; @todo
;
; zminmesh and izmaxmesh may be set to 0,0 because we only work with
; surface problem (no z) whatever the resolution of ORCA. this will safe lot
; of memory
;
; use <pro>ncdf_getaxis</pro> to recover meshmask dimension
; on x, retrieve -1 and -1 (bande de recouvrement haut et bas)
; on y, for ORCA2, retrieve -2 (bande de recouvrement ++ et ++)
;
; reminder : ORCA2 = 182 x 149 and ORCA025 1442 x 1021
;
; for ORCA025, use /keepcont in <pro>triangule</pro> call.
;
; @version
; $Id$
;
;-
;
PRO initocemesh, orcares, DRAKKAR_EXP = drakkar_exp, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4cal
@cm_4data
;
; check input parameters
CASE orcares OF
   'ORCA2': BEGIN
               msg = 'iii : valid orcares parameter = ' + orcares
               ras = report(msg)
               filename_oce = 'meshmask_bab.nc'
               IF keyword_set(DRAKKAR_EXP) THEN BEGIN
                  msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp
                  ras = report(msg)
               END
            END
     'ORCA025': BEGIN
                 msg = 'iii : valid orcares parameter = ' + orcares
                 ras = report(msg)
                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
                    msg = 'iii : DRAKKAR_EXP keyword set'
                    ras = report(msg)
                    msg = 'iii : DRAKKAR_EXP = ' + drakkar_exp
                    ras = report(msg)
                    CASE drakkar_exp OF
                       'G42' : BEGIN
                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
                                  ras = report(msg)
                               END
                       'G70' : BEGIN
                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
                                  ras = report(msg)
                               END
                       ELSE : BEGIN
                                  msg = 'eee : invalid DRAKKAR_EXP keyword = ' + drakkar_exp
                                  ras = report(msg)
                                  RETURN
                              END
                    ENDCASE
                    filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc'
                 END
              END
    ELSE  : BEGIN
               msg = 'eee : invalid orcares parameter = ' + orcares
               ras = report(msg)
               msg = 'eee : orcares must be ORCA2 or ORCA025'
               ras = report(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'
              ras = report(msg)
              RETURN
            END
     ELSE: BEGIN
             msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env
             ras = report(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.'
     ras = report(msg)
     RETURN
  ENDIF
;
; 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.'
     ras = report(msg)
     RETURN
  ENDIF
;
; protection
  IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN
     msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.'
     ras = report(msg)
     RETURN
  ENDIF
;
;--------------------------------------------------------
; boundaries of the read grid regarding to the original grid
;----------------------------------------------------------
;
CASE orcares OF
   'ORCA2': BEGIN
ixminmesh  =1
ixmaxmesh  =180
;
iyminmesh  =0
iymaxmesh  =147
;
izminmesh  =-1
izmaxmesh  =-1
;
              END
   'ORCA025': BEGIN
ixminmesh  =1
ixmaxmesh  =1440
;
iyminmesh  =0
iymaxmesh  =1019
;
izminmesh  =-1
izmaxmesh  =-1
;
              END
      ELSE  : BEGIN
                 msg = 'eee : invalid orcares parameter = ' + orcares
                 ras = report(msg)
                 msg = 'eee : orcares must be ORCA2 or ORCA025'
                 ras = report(msg)
                 RETURN
              END
  ENDCASE
jpt = 1
time = 0
;------------------------------------------------------
; read the grid
;------------------------------------------------------
;
key_stride = [1, 1, 1]
ncdf_meshlec, fullfilename_oce[0], glamboundary = [20, 380], _extra = ex
;-------------------------------------------------------------
domdef
;
CASE orcares OF
  'ORCA2'  : BEGIN
                triangles_list = triangule()
             END
  'ORCA025': BEGIN
                triangles_list = triangule(/keepcont)
             END
      ELSE : BEGIN
                msg = 'eee : invalid orcares parameter = ' + orcares
                ras = report(msg)
                msg = 'eee : orcares must be ORCA2 or ORCA025'
                ras = report(msg)
                RETURN
             END
ENDCASE
;
; boundaries of the data regarding to the original grid
jpidta = jpiglo
jpjdta = jpjglo
jpkdta = jpkglo
ixmindta = 0
ixmaxdta = jpidta-1
iymindta = 0
iymaxdta = jpjdta-1
izmindta = 0
izmaxdta = jpkdta-1
;
@updateold
;
END