;+ ; ; @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 supprimer 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 ; ; @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 sur lequel 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