source: trunk/SRC/Obsolete/meshlec.pro @ 134

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments 
; lecture du mask des sorties d'OPA. les sources se trouvent ds les
; repertoires sur maia du type:
; /nom_exp/RESTARTS
;
; @obsolete
;
; @examples
; IDL> meshmask[,' nomfich']
;
; @param nomfich {in}{required} string, c''est le nom du fichier a lire. Par defaut, c''est meshmask.
;
; @keyword GLAMBOUNDARY {in} un vecteur de 2 elements specifaint le min et le
; max qui doivent etre imposes en longitude (obligatoire si le
; tableau depasse 360 degres)
;
; @keyword pasblabla {in} pour suprimer les blablas
;
; @keyword DOUBLE {in} pour forcer a lire les tableaux en double
;        precision. ce Mot clef est maintenant active
;        automatiquement.
;
; @keyword GETDIMENSIONS
;
; @uses common.pro
;
; @restrictions
;  La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh
;  ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette
;  routine. pour attribuer automatiquement ces valeurs au maximum
;  possible les mettre toutes a -1 et meshlec les calculera.
;
; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
;     Marina Levy : lecture en double precision (cas calcul sur shine)
;
; @version $Id$
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
pro meshlec, nomfich, PASBLABLA = pasblabla,DOUBLE = double, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = GETDIMENSIONS
;
  compile_opt idl2, strictarrsubs, obsolete
;
@common
   tempsun = systime(1)         ; pour key_performance
;---------------------------------------------------------
;---------------------------------------------------------
   jpiglo      = 0L
   jpjglo      = 0L
   jpkglo      = 0L
   tab    = 'aaaaa'
;---------------------------------------------------------
;  definition du domaine de la grille surlequel sont
;  effectuees les sorties.les indices des tableaux commencant
;  a 1: cf le fichier wrivr2.F ds WKOPA sur le cray
;----------------------------------------------------------
; LECTURE DU MASK trouve ds les fichiers restart
;----------------------------------------------------------
;----------------------------------------------------------
;
;----------------------------------------------------------
;       constitution de l'adresse s_fichier et
;     ouverture du fichier a l'adresse s_fichier
;-------------------------------------------------------

   IF n_params() EQ 0 then nomfich = 'meshmask'
   s_fichier = isafile(file = nomfich, iodir = iodir)
   if not keyword_set(pasblabla) then print, ' '
   if not keyword_set(pasblabla) then print,'adresse du fichier: ',s_fichier

   openr, numlec, s_fichier, /get_lun, /f77_unformatted, /swap_if_little_endian
   filepamameters = fstat(numlec)
;---------------------------------------------------------
;     lecture
;---------------------------------------------------------
   readu, numlec, jpiglo, jpjglo, jpkglo
   if not keyword_set(pasblabla) then print, 'taille de la grille d''origine: ',jpiglo,', ',jpjglo,', ',jpkglo
;
   if keyword_set(getdimensions) then begin
      free_lun,numlec
      close, numlec
      return
   endif

;---------------------------------------------------------
; on determine si le fichier a ete ecrit en double precision on non...
   sizenumber = 8l
   sizefile8 = 4l+3l*4l+4l $
    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
    +4l*(4l+jpiglo*jpjglo*jpkglo*sizenumber+4l) $
    +1l*(4l+jpiglo*jpjglo*sizenumber+4l) $
    +4l*(4l+jpkglo*sizenumber+4l)
   if filepamameters.size GE sizefile8 THEN double = 1
;    sizenumber = 4l
;    sizefile4 = 4l+3l*4l+4l $
;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
;     +4l*(4l+jpiglo*jpjglo*jpkglo*sizenumber+4l) $
;     +1l*(4l+jpiglo*jpjglo*sizenumber+4l) $
;     +4l*(4l+jpkglo*sizenumber+4l)
; print, filepamameters.size,  sizefile4,  sizefile8
;    case filepamameters.size of
;       sizefile8:double = 1
;       sizefile4:double = 0
;       ELSE:BEGIN
;          nothing = report('The OPA Mesh file as not the good size!')
;          free_lun,numlec
;          close, numlec
;          return
;       END
;    endcase
;
;
   if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0
   if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1
   if ixminmesh EQ -1 THEN ixminmesh = 0
   IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1
   if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0
   IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1
   if iyminmesh EQ -1 THEN iyminmesh = 0
   IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1
   if n_elements(izminmesh) EQ 0 THEN izminmesh = 0
   IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1
   if izminmesh EQ -1 THEN izminmesh = 0
   IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1
;
;

   jpi    = long(ixmaxmesh-ixminmesh+1)
   jpj    = long(iymaxmesh-iyminmesh+1)
   jpk    = long(izmaxmesh-izminmesh+1)
;-------------------------------------------------------
; doit-on reellement lire la grille?
;-------------------------------------------------------
   meshparameters = {jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $
                     , jpi:jpi, jpj:jpj, jpk:jpk $
                     , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $
                     , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
                     , izminmesh:izminmesh, izmaxmesh:izmaxmesh $
                     , key_shift:key_shift}
;
;-------------------------------------------------------
   noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...')
;-------------------------------------------------------

   IF NOT keyword_set(double) THEN BEGIN
      z3d    = fltarr(jpiglo, jpjglo, jpkglo)
      z2d    = fltarr(jpiglo, jpjglo)
      z1d    = fltarr(jpkglo)
   ENDIF ELSE BEGIN
      z3d    = dblarr(jpiglo, jpjglo, jpkglo)
      z2d    = dblarr(jpiglo, jpjglo)
      z1d    = dblarr(jpkglo)
   ENDELSE

   if not keyword_set(pasblabla) then print, ' '
   readu, numlec, tab,z2d
   GLAMT=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMT[25,31]: ',GLAMT[25,31]
   readu, numlec, tab,z2d
   GLAMU=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMU[25,31]: ',GLAMU[25,31]
   readu, numlec, tab,z2d
   GLAMV=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMV[25,31]: ',GLAMV[25,31]
   readu, numlec, tab,z2d
   GLAMF=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMF[25,31]: ',z2d[25,31]

   if not keyword_set(pasblabla) then print, ' '
   readu, numlec, tab, z2d
   GPHIT=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIT[25,31]: ',GPHIT[25,31]
   readu, numlec, tab, z2d
   GPHIU=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIU[25,31]: ',GPHIU[25,31]
   readu, numlec, tab, z2d
   GPHIV=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIV[25,31]: ',GPHIV[25,31]
   readu, numlec, tab, z2d
   GPHIF=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIF[25,31]: ',z2d[25,31]

   if not keyword_set(pasblabla) then print, ' '
   readu, numlec, tab, z2d
   E1T=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1T[25,5]: ',z2d[25,5]
   readu, numlec, tab, z2d
   E1U=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1U[25,5]: ',z2d[25,5]
   readu, numlec, tab, z2d
   E1V=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1V[25,5]: ',z2d[25,5]
   readu, numlec, tab, z2d
   E1F=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1F[25,5]: ',z2d[25,5]

   if not keyword_set(pasblabla) then print, ' '
   readu, numlec, tab, z2d
   E2T=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2T[25,5]: ',z2d[25,5]
   readu, numlec, tab, z2d
   E2U=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2U[25,5]: ',z2d[25,5]
   readu, numlec, tab, z2d
   E2V=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2V[25,5]: ',z2d[25,5]
   readu, numlec, tab, z2d
   E2F=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2F[25,5]: ',z2d[25,5]

   if not keyword_set(pasblabla) then print, ' '
   readu, numlec, tab, z3d
   TMASK=byte(z3d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh,izminmesh:izmaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur TMASK[25,5,0]: ',TMASK[25,5,0]
   readu, numlec, tab, z3d
   UMASKred=byte(z3d[ixmaxmesh,iyminmesh:iymaxmesh,izminmesh:izmaxmesh])
   umaskred = reform(umaskred)
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur UMASK[25,5,0]: ',z3d[25,5,0]
   readu, numlec, tab, z3d
   VMASKred=byte(z3d[ixminmesh:ixmaxmesh,iymaxmesh,izminmesh:izmaxmesh])
   vmaskred = reform(vmaskred)
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur VMASK[25,5,0]: ',z3d[25,5,0]
   readu, numlec, tab, z3d
   fmaskredy=byte(z3d[ixmaxmesh,iyminmesh:iymaxmesh,izminmesh:izmaxmesh])
   coast = where(fmaskredy NE 0 and fmaskredy NE 1)
   IF coast[0] NE -1 THEN fmaskredy[coast] = 0b
   fmaskredx=byte(z3d[ixminmesh:ixmaxmesh,iymaxmesh,izminmesh:izmaxmesh])
   coast = where(fmaskredx NE 0 and fmaskredx NE 1)
   IF coast[0] NE -1 THEN fmaskredx[coast] = 0b
   fmaskredx = reform(fmaskredx)
   fmaskredy = reform(fmaskredy)
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur FMASK[25,5,0]: ',z3d[25,5,0]
;
   if not keyword_set(pasblabla) then print, ' '
   readu, numlec, tab, z2d
;FF=z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh]
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur FF[25,5]: ',z2d[25,5]
   readu, numlec, tab, z1d
   GDEPT=float(z1d[izminmesh:izmaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GDEPT[1]: ',GDEPT[1]
   readu, numlec, tab, z1d
   GDEPW=float(z1d[izminmesh:izmaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GDEPW[1]: ',GDEPW[1]
   readu, numlec, tab, z1d
   E3T=float(z1d[izminmesh:izmaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E3T[3]: ',E3T[3]
   readu, numlec, tab, z1d
   E3W=float(z1d[izminmesh:izmaxmesh])
   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E3W[3]: ',E3W[3]
;-------------------------------------------------------
   free_lun,numlec
   close, numlec
;-------------------------------------------------------
;-------------------------------------------------------
; bornes de glam qui ne doivent pas depasser 360 degres..
;-------------------------------------------------------
;    minglam = min(glamt, max = maxglam)
;    if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $
;     nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget))
;-------------------------------------------------------
   if keyword_set(glamboundary) then begin
      if glamboundary[0] NE glamboundary[1] then begin
         glamt = glamt MOD 360
         smaller = where(glamt LT glamboundary[0])
         if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360
         bigger = where(glamt GE glamboundary[1])
         if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360
         glamu = glamu MOD 360
         smaller = where(glamu LT glamboundary[0])
         if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360
         bigger = where(glamu GE glamboundary[1])
         if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360
         glamv = glamv MOD 360
         smaller = where(glamv LT glamboundary[0])
         if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360
         bigger = where(glamv GE glamboundary[1])
         if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360
         glamf = glamf MOD 360
         smaller = where(glamf LT glamboundary[0])
         if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360
         bigger = where(glamf GE glamboundary[1])
         if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360
      endif
   endif
;-------------------------------------------------------
; shift en x
;-------------------------------------------------------
   if keyword_set(key_shift) AND jpi NE 1 then begin
      glamt = shift(glamt,key_shift , 0)
      gphit = shift(gphit,key_shift , 0)
      e1t = shift(e1t, key_shift, 0)
      e2t = shift(e2t,key_shift , 0)
      glamu = shift(glamu, key_shift, 0)
      gphiu = shift(gphiu, key_shift, 0)
      e1u = shift(e1u,key_shift , 0)
      e2u = shift(e2u, key_shift, 0)
      glamv = shift(glamv, key_shift, 0)
      gphiv = shift(gphiv, key_shift, 0)
      e1v = shift(e1v,key_shift , 0)
      e2v = shift(e2v, key_shift, 0)
      glamf = shift(glamf, key_shift, 0)
      gphif = shift(gphif, key_shift, 0)
      e1f = shift(e1f, key_shift , 0)
      e2f = shift(e2f, key_shift, 0)
      if jpk EQ 1 then begin
         tmask = shift(tmask, key_shift, 0)
         vmaskred = shift(vmaskred, key_shift)
         fmaskredx = shift(fmaskredx, key_shift)
      ENDIF ELSE BEGIN
         tmask = shift(tmask, key_shift, 0,0)
         vmaskred = shift(vmaskred, key_shift, 0)
         fmaskredx = shift(fmaskredx, key_shift, 0)
      ENDELSE
   endif
;-------------------------------------------------------
   key_yreverse = 0
   key_zreverse = 0
   key_partialstep = 0
   key_stride = [1, 1, 1]
   key_gridtype = 'c'
;
   if not keyword_set(pasblabla) then print,'lecture '+nomfich+' finie'
   widget_control, noticebase, bad_id = toto, /destroy
   if keyword_set(key_performance) THEN print, 'temps meshlec', systime(1)-tempsun

   return
end