Ignore:
Timestamp:
08/30/07 14:44:23 (17 years ago)
Author:
smasson
Message:

bugfix for interpolation from ORCA2 without mask

Location:
trunk/SRC/ToBeReviewed/LECTURE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro

    r240 r271  
    3131; Useless, defined for compatibility 
    3232; 
     33; @keyword ADDSCL_BEFORE {default=0}{type=scalar: 0 or 1} 
     34; put 1 to apply add_offset ad scale factor on data before looking for 
     35; missing values 
     36; 
    3337; @keyword BOXZOOM 
    3438; Contain the boxzoom on which we have to do the reading 
    … …  
    4145; 
    4246; @keyword INIT {default=0}{type=scalar: 0 or 1} 
    43 ; To call automatically initncdf, filename and thus 
     47; To call automatically initncdf with filename as input argument and thus 
    4448; redefine all the grid parameters 
    4549; 
    … …  
    6266; but only the array referring to the field. 
    6367; 
    64 ; @keyword TIMEVAR {type=string} 
    65 ; It define the name of the variable that 
    66 ; contains the time axis. This keyword can be useful if there 
    67 ; is no unlimited dimension or if the time axis selected by default 
    68 ; (the first 1D array with unlimited dimension) is not the good one. 
    69 ; 
    7068; @keyword ZETAFILENAME {default=FILENAME}{type=string} 
    7169; For ROMS outputs. The filename of the file where zeta variable should be read 
    … …  
    7573; 
    7674; @keyword _EXTRA 
    77 ; Used to pass keywords 
     75; Used to pass keywords to <pro>isafile</pro>, <pro>initncdf</pro>, 
     76; <pro>ncdf_gettime</pro> and <pro>domdef</pro> 
    7877; 
    7978; @returns 
    … …  
    9594; 
    9695FUNCTION read_ncdf, name, beginning, ending, compatibility, BOXZOOM = boxzoom, FILENAME = filename $ 
    97                     , PARENTIN = parentin, TIMESTEP = timestep, TIMEVAR = timevar $ 
     96                    , PARENTIN = parentin, TIMESTEP = timestep, ADDSCL_BEFORE = addscl_before $ 
    9897                    , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $ 
    9998                    , GRID = grid, CALLITSELF = callitself $ 
    … …  
    115114;  print,filename 
    116115; is parent a valid widget ? 
    117   if keyword_set(parentin) then BEGIN 
     116  IF keyword_set(parentin) THEN BEGIN 
    118117    parent = long(parentin) 
    119118    parent = parent*widget_info(parent, /managed) 
    … …  
    123122; Opening of the name file 
    124123;------------------------------------------------------------ 
    125   if size(filename, /type) NE 7 then $ 
     124  IF size(filename, /type) NE 7 THEN $ 
    126125    return, report('read_ncdf cancelled') 
    127126  IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null' 
    128127  cdfid = ncdf_open(filename) 
    129   contient = ncdf_inquire(cdfid) 
     128  inq = ncdf_inquire(cdfid) 
    130129;------------------------------------------------------------ 
    131130; we check if the variable name exists in the file. 
    132131;------------------------------------------------------------ 
    133   if ncdf_varid(cdfid, name) EQ -1 then BEGIN 
     132  IF ncdf_varid(cdfid, name) EQ -1 THEN BEGIN 
    134133    ncdf_close, cdfid 
    135134    return, report('variable '+name+' !C not found in the file '+filename) 
    136135  ENDIF 
    137   varcontient = ncdf_varinq(cdfid, name) 
    138   IF varcontient.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data') 
     136  varinq = ncdf_varinq(cdfid, name) 
     137  IF varinq.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data') 
    139138; look for the dimension names 
    140   dimnames = strarr(varcontient.ndims) 
    141   FOR i = 0, varcontient.ndims-1 DO BEGIN 
    142     ncdf_diminq, cdfid, varcontient.dim[i], tmp, dimsize 
     139  dimnames = strarr(varinq.ndims) 
     140  FOR i = 0, varinq.ndims-1 DO BEGIN 
     141    ncdf_diminq, cdfid, varinq.dim[i], tmp, dimsize 
    143142    dimnames[i] = tmp 
    144143  ENDFOR 
    … …  
    146145; shall we redefine the grid parameters 
    147146;------------------------------------------------------------ 
    148   if keyword_set(init) THEN initncdf, filename, _extra = ex 
     147  IF keyword_set(init) THEN initncdf, filename, _extra = ex 
    149148;------------------------------------------------------------ 
    150149; check the time axis and the debut and ending dates 
    151150;------------------------------------------------------------ 
    152   if n_elements(beginning) EQ 0 then begin 
     151  IF n_elements(beginning) EQ 0 THEN BEGIN 
    153152    beginning = 0 
    154153    timestep = 1 
    155   endif 
    156   if keyword_set(timestep) then begin 
    157     firsttps = beginning[0] 
    158     if n_elements(ending) NE 0 then lasttps = ending[0] ELSE lasttps = firsttps 
    159     jpt = lasttps-firsttps+1 
    160     IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt) 
    161   ENDIF ELSE BEGIN 
    162     if keyword_set(parent) then BEGIN 
     154  ENDIF 
     155; define time and jpt 
     156  CASE 1 OF 
     157    keyword_set(timestep):BEGIN  
     158      firsttps = beginning[0] 
     159      IF n_elements(ending) NE 0 THEN lasttps = ending[0] ELSE lasttps = firsttps 
     160      jpt = lasttps-firsttps+1 
     161      IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt) 
     162    END 
     163    keyword_set(parent):BEGIN  
    163164      widget_control, parent, get_uvalue = top_uvalue 
    164165      filelist = extractatt(top_uvalue, 'filelist') 
    … …  
    168169      date1 = date2jul(beginning[0]) 
    169170      if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 
    170       firsttps = where(time EQ date1) & firsttps = firsttps[0] 
    171       lasttps = where(time EQ date2) & lasttps = lasttps[0] 
    172     ENDIF ELSE BEGIN 
    173       IF keyword_set(timevar) THEN BEGIN 
    174         timeid = ncdf_varid(cdfid, timevar) 
    175         IF timeid EQ -1 THEN BEGIN 
    176           ncdf_close, cdfid 
    177           return, report('the file '+filename+' as no variable ' + timevar $ 
    178                          + '. !C Use the TIMESTEP keyword') 
    179         endif 
    180         timecontient = ncdf_varinq(cdfid, timeid) 
    181         contient.recdim = timecontient.dim[0] 
    182       ENDIF ELSE BEGIN 
    183 ; we find the infinite dimension 
    184         timedim = contient.recdim 
    185         if timedim EQ -1 then BEGIN 
    186           ncdf_close, cdfid 
    187           return, report('the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword') 
    188         endif 
    189 ; we find the FIRST time axis 
    190         timeid = 0 
    191         repeat BEGIN       ; As long as we have not find a variable having only one dimension: the infinite one 
    192           timecontient = ncdf_varinq(cdfid, timeid) ; that the variable contain. 
    193           timeid = timeid+1 
    194         endrep until (n_elements(timecontient.dim) EQ 1 $ 
    195                       AND timecontient.dim[0] EQ contient.recdim) $ 
    196           OR timeid EQ contient.nvars+1 
    197 ; 
    198         if timeid EQ contient.nvars+1 then BEGIN 
    199           ncdf_close, cdfid 
    200           return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword') 
    201         endif 
    202         timeid = timeid-1 
    203       ENDELSE 
    204 ; we must found the time origin of the julian calendar used in the 
    205 ; time axis. 
    206 ; does the attribut units an dcalendar exist for the variable time axis? 
    207       if timecontient.natts EQ 0 then BEGIN 
     171      firsttps = (where(abs(time - date1) LT 0.9d/86400.d))[0]   ; beware of rounding errors... 
     172      lasttps = (where(abs(time - date2) LT 0.9d/86400.d))[0] 
     173      jpt = lasttps-firsttps+1 
     174    END 
     175    ELSE:BEGIN  
     176      time = ncdf_gettime(filename, cdfid, caller = 'read_ncdf', err = err, _extra = ex) 
     177      IF n_elements(err) NE 0 THEN BEGIN 
     178        dummy = report(err) 
    208179        ncdf_close, cdfid 
    209         return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable') 
    210       endif 
    211       attnames = strarr(timecontient.natts) 
    212       for attiq = 0, timecontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, timeid, attiq) 
    213       if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN 
    214         ncdf_close, cdfid 
    215         return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword') 
     180        return, -1 
    216181      ENDIF 
    217 ; 
    218 ; now we try to find the attribut called calendar... 
    219 ; the attribute "calendar" exists? 
    220 ; If no, we suppose that the calendar is gregorian calendar 
    221 ; 
    222       if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN 
    223         ncdf_attget, cdfid, timeid, 'calendar', value 
    224         value = string(value) 
    225         CASE value OF 
    226           'noleap':key_caltype = 'noleap' 
    227           '360d':key_caltype = '360d' 
    228           'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' 
    229           ELSE:BEGIN 
    230 ;            notused = report('Unknown calendar: '+value+', we use greg calendar.') 
    231             key_caltype = 'greg' 
    232           END 
    233         ENDCASE 
    234       ENDIF ELSE BEGIN 
    235 ;        notused = report('Unknown calendar, we use '+key_caltype+' calendar.') 
    236         IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' 
    237       ENDELSE 
    238 ; 
    239 ; now we take acre of units attribut 
    240       ncdf_attget, cdfid, timeid, 'units', value 
    241 ; 
    242 ; time_counter:units = "seconds since 0001-01-01 00:00:00" ; 
    243 ; time_counter:units = "hours since 0001-01-01 00:00:00" ; 
    244 ; time_counter:units = "days since 1979-01-01 00:00:00" ; 
    245 ; time_counter:units = "months since 1979-01-01 00:00:00" ; 
    246 ; time_counter:units = "years since 1979-01-01 00:00:00" ; 
    247 ; 
    248 ; we decript the "units" attribut to find the time origin 
    249       value = strtrim(strcompress(string(value)), 2) 
    250       mots = str_sep(value, ' ') 
    251       unite = mots[0] 
    252       unite = strlowcase(unite) 
    253       IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1) 
    254       IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7) 
    255       IF unite NE 'second' AND unite NE 'hour' AND unite NE 'day' $ 
    256          AND unite NE 'month' AND unite NE 'year' THEN BEGIN 
    257         ncdf_close, cdfid 
    258         return, report('time units does not start with seconds/hours/days/months/years') 
    259       ENDIF 
    260       IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN 
    261         ncdf_close, cdfid 
    262         return, report('attribut units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*') 
    263       ENDIF 
    264       depart = str_sep(mots[2], '-') 
    265       ncdf_varget, cdfid, timeid, time 
    266       time = double(time) 
    267       case unite of 
    268         'second':time = julday(depart[1], depart[2], depart[0])+time/86400.d 
    269         'hour':time = julday(depart[1], depart[2], depart[0])+time/24.d 
    270         'day':time = julday(depart[1], depart[2], depart[0])+time 
    271         'month':BEGIN 
    272           if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m 
    273             time = julday(depart[1], depart[2], depart[0])+round(time*30) $ 
    274             ELSE for t = 0, n_elements(time)-1 DO $ 
    275             time[t] = julday(depart[1]+time[t], depart[2], depart[0]) 
    276         END 
    277         'year':BEGIN 
    278           if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y 
    279             time = julday(depart[1], depart[2], depart[0])+round(time*365) $ 
    280             ELSE for t = 0, n_elements(time)-1 do $ 
    281             time[t] = julday(depart[1], depart[2], depart[0]+time[t]) 
    282         END 
    283         ELSE:BEGIN 
    284           ncdf_close, cdfid 
    285           return, report('The "units" attribute of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."') 
    286         end 
    287       ENDCASE 
     182; date1 
    288183      date1 = date2jul(beginning[0]) 
     184; date2 
    289185      if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 
    290       time = double(time) 
    291       firsttps = where(time GE date1) & firsttps = firsttps[0] 
     186; firsttps 
     187      firsttps = where(time GE (date1 - 0.9d/86400.d)) & firsttps = firsttps[0] 
    292188      if firsttps EQ -1 THEN BEGIN 
    293189        ncdf_close, cdfid 
    294190        return, report('date 1: '+strtrim(jul2date(date1), 1)+' is not found in the time axis.') 
    295191      ENDIF 
    296       lasttps = where(time LE date2) 
    297       if lasttps[0] EQ -1 THEN BEGIN 
     192; lasttps 
     193      lasttps = where(time LE (date2 + 0.9d/86400.d)) & lasttps = lasttps[n_elements(lasttps)-1] 
     194      if lasttps EQ -1 THEN BEGIN 
    298195        ncdf_close, cdfid 
    299196        return, report('the time axis has no date before date 2: '+strtrim(jul2date(date2), 1)) 
    300197      endif 
    301       lasttps = lasttps[n_elements(lasttps)-1] 
    302198      if lasttps LT firsttps then BEGIN 
    303199        ncdf_close, cdfid 
    304200        return, report('the time axis has no dates between date1 and  date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1)) 
    305201      endif 
    306     ENDELSE 
    307     time = time[firsttps:lasttps] 
    308     jpt = lasttps-firsttps+1 
    309   ENDELSE 
     202      time = time[firsttps:lasttps] 
     203      jpt = lasttps-firsttps+1 
     204    END 
     205  ENDCASE 
    310206;------------------------------------------------------------ 
    311207; Name of the grid on which the field refer to. 
    … …  
    316212; are we in one of the case corresponding to ROMS conventions? 
    317213      CASE 1 OF 
    318         dimnames[2 <(varcontient.ndims-1)] EQ 's_w':vargrid = 'W' 
     214        dimnames[2 <(varinq.ndims-1)] EQ 's_w':vargrid = 'W' 
    319215        dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_rho':vargrid = 'T' 
    320216        dimnames[0] EQ 'xi_u'   AND dimnames[1] EQ 'eta_u'  :vargrid = 'U' 
    … …  
    342238  ENDELSE 
    343239;--------------------------------------------------------------- 
    344 ; call the init function ? 
    345 ;--------------------------------------------------------------- 
    346  
    347 ;--------------------------------------------------------------- 
    348240; redefinition of the  domain 
    349241;--------------------------------------------------------------- 
    350   if keyword_set(tout) then begin 
     242  IF keyword_set(tout) THEN BEGIN 
    351243    nx = jpi 
    352244    ny = jpj 
    … …  
    413305;--------------------------------------------------------------------- 
    414306;--------------------------------------------------------------------- 
    415 ; We define global variable joined with the variable. 
     307; We define common (cm_4data) variables associated with the variable. 
    416308;--------------------------------------------------------------------- 
    417309; varname 
    418310  IF NOT keyword_set(callitself) THEN varname = name 
    419 ; varunit 
    420   if varcontient.natts NE 0 then begin 
    421     attnames = strarr(varcontient.natts) 
    422     for attiq = 0, varcontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, name, attiq) 
    423     lowattnames = strlowcase(attnames) 
    424 ; 
    425     found = (where(lowattnames EQ 'units'))[0] 
    426     IF found NE -1 then ncdf_attget, cdfid, name, attnames[found], value ELSE value = '' 
    427     IF NOT keyword_set(callitself) THEN varunit = strtrim(string(value), 2) 
    428 ; 
    429     found = (where(lowattnames EQ 'add_offset'))[0] 
    430     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], add_offset ELSE add_offset = 0. 
    431 ; 
    432     found = (where(lowattnames EQ 'scale_factor'))[0] 
    433     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], scale_factor ELSE scale_factor = 1. 
    434 ; 
    435     missing_value = 'no' 
    436     found = (where(lowattnames EQ '_fillvalue'))[0] 
    437     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value 
    438     found = (where(lowattnames EQ 'missing_value'))[0] 
    439     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value 
    440 ; 
    441   ENDIF ELSE BEGIN 
    442     IF NOT keyword_set(callitself) THEN varunit = '' 
    443     add_offset = 0. 
    444     scale_factor = 1. 
    445     missing_value = 'no' 
    446   ENDELSE 
     311; varunit and others... 
     312  ncdf_getatt, cdfid, name, add_offset = add_offset, scale_factor = scale_factor, missing_value = missing_value, units = units 
     313  IF NOT keyword_set(callitself) THEN varunit = units 
    447314; vardate 
    448315; We make a legible date in function of the specified language. 
    … …  
    455322; We apply the value valmask on land points. 
    456323  if NOT keyword_set(cont_nofill) then begin 
    457     valmask = 1e20 
     324    valmask = 1.e20 
    458325    case 1 of 
    459       varcontient.ndims eq 2:BEGIN ;xy array 
    460         mask = mask[*, *, 0] 
     326      varinq.ndims eq 2:BEGIN ;xy array 
     327        earth = where(mask[*, *, 0] EQ 0) 
     328      END 
     329      varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] EQ -1:BEGIN ;xyz array 
    461330        earth = where(mask EQ 0) 
    462331      END 
    463       varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array 
    464         earth = where(mask EQ 0) 
    465       END 
    466       varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array 
    467         mask = mask[*, *, 0] 
    468         earth = where(mask EQ 0) 
     332      varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1:BEGIN ;xyt array     
     333        earth = where(mask[*, *, 0] EQ 0) 
    469334        if earth[0] NE -1 then BEGIN 
    470335          earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt) 
    471336        END 
    472337      END 
    473       varcontient.ndims eq 4:BEGIN ;xyzt array 
     338      varinq.ndims eq 4:BEGIN ;xyzt array 
    474339        earth = where(mask EQ 0) 
    475340        if earth[0] NE -1 then BEGIN 
    … …  
    477342        END 
    478343      END 
    479     endcase 
     344    ENDCASE 
    480345  ENDIF ELSE earth = -1 
    481346; we look for  missing_value 
     347; we apply add_offset, scale_factor and missing_value 
     348  IF keyword_set(addscl_before) THEN BEGIN 
     349    if scale_factor NE 1 then res = temporary(res)*scale_factor 
     350    if add_offset NE 0 then res = temporary(res)+add_offset 
     351  ENDIF 
    482352  IF size(missing_value, /type) NE 7 then BEGIN 
    483     IF size(missing_value, /type) EQ 1 THEN BEGIN 
    484       missing_value = strlowcase(string(missing_value)) 
    485       IF strmid(missing_value, 0, 1, /reverse_offset) EQ 'f' THEN $ 
    486         missing_value = strmid(missing_value, 0, strlen(missing_value)-1) 
    487       IF isnumber(string(missing_value), tmp) EQ 1 THEN missing_value = tmp ELSE BEGIN 
    488         ras = report('Warning: missing value is not a number: ' + missing_value) 
    489         missing_value = - 1 
    490       ENDELSE 
    491     ENDIF 
    492 ;    if missing_value NE valmask then begin 
    493     if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $ 
    494     ELSE missing = where(abs(res) gt abs(missing_value)/10.) 
    495 ;    ENDIF ELSE missing = -1 
     353    IF finite(missing_value) EQ 1 THEN BEGIN 
     354        CASE 1 OF 
     355          missing_value GT 1.e6:missing = where(res GT missing_value-10.) 
     356          missing_value LT -1.e6:missing = where(res LT missing_value+10.) 
     357          abs(missing_value) LT 1.e-6:missing = where(res LT 1.e-6) 
     358          ELSE:missing = where(res EQ missing_value) 
     359        ENDCASE 
     360    ENDIF ELSE missing = where(finite(res) EQ 0) 
    496361  ENDIF ELSE missing = -1 
    497 ; we apply add_offset, scale_factor and missing_value 
    498   if scale_factor NE 1 then res = temporary(res)*scale_factor 
    499   if add_offset NE 0 then res = temporary(res)+add_offset 
     362  IF NOT keyword_set(addscl_before) THEN BEGIN 
     363    if scale_factor NE 1 then res = temporary(res)*scale_factor 
     364    if add_offset NE 0 then res = temporary(res)+add_offset 
     365  ENDIF 
    500366  if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan 
    501367  if earth[0] NE -1 then res[temporary(earth)] = 1.e20 
    … …  
    508374      ncdf_attget, cdfid, 'hc', hc, /global 
    509375; look for all variables names 
    510       allvarnames = strarr(contient.nvars) 
    511       FOR i = 0, contient.nvars-1 DO BEGIN 
     376      allvarnames = strarr(inq.nvars) 
     377      FOR i = 0, inq.nvars-1 DO BEGIN 
    512378        tmp = ncdf_varinq( cdfid, i) 
    513379        allvarnames[i] = tmp.name 

  • trunk/SRC/ToBeReviewed/LECTURE/read_ncdf_varget.pro

    r231 r271  
    9797;...................................................................... 
    9898;...................................................................... 
    99    varcontient.ndims eq 2:BEGIN ;xy array 
     99   varinq.ndims eq 2:BEGIN ;xy array 
    100100;...................................................................... 
    101101;...................................................................... 
    … …  
    185185;...................................................................... 
    186186;...................................................................... 
    187    varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array 
     187   varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] EQ -1:BEGIN ;xyz array 
    188188;...................................................................... 
    189189;...................................................................... 
    … …  
    276276;...................................................................... 
    277277;...................................................................... 
    278    varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array 
     278   varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1:BEGIN ;xyt array 
    279279;...................................................................... 
    280280;...................................................................... 
    … …  
    376376;...................................................................... 
    377377;...................................................................... 
    378    varcontient.ndims eq 4:BEGIN ;xyzt array 
     378   varinq.ndims eq 4:BEGIN ;xyzt array 
    379379;...................................................................... 
    380380;...................................................................... 
    … …  
    481481; we apply reverse 
    482482  IF keyword_set(key_yreverse) AND ny NE 1 THEN BEGIN 
    483     IF varcontient.ndims - ((where(varcontient.dim EQ contient.recdim))[0] NE -1) EQ 2 THEN $ 
     483    IF varinq.ndims - ((where(varinq.dim EQ inq.recdim))[0] NE -1) EQ 2 THEN $ 
    484484       res = reverse(reform(res, nx, ny, jpt, /overwrite),  2) $ 
    485485    ELSE res = reverse(reform(res, nx, ny, nz, jpt, /overwrite),  2) 
    486486  ENDIF 
    487487  if keyword_set(key_zreverse) AND nz NE 1 $ 
    488      AND varcontient.ndims - ((where(varcontient.dim EQ contient.recdim))[0] NE -1) EQ 3 THEN $ 
     488     AND varinq.ndims - ((where(varinq.dim EQ inq.recdim))[0] NE -1) EQ 3 THEN $ 
    489489        res = reverse(reform(res, nx, ny, nz, jpt, /overwrite),  3) 
    490490; 
Note: See TracChangeset for help on using the changeset viewer.