Changeset 81


Ignore:
Timestamp:
03/26/08 20:24:00 (16 years ago)
Author:
kolasinski
Message:

First temporary commit for the new reading procedures - Still some bugs

Location:
branches/procs-read
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • branches/procs-read/data_read.pro

    r60 r81  
    6666      ENDIF       
    6767   ENDIF  
    68    varexp = cmd.exp    
     68;   varexp = cmd.exp    
    6969 
    7070; define data base 
    … …  
    9191   ENDIF  
    9292 
    93 ; define data domain 
    94   
    95    IF hotyp NE '-' THEN BEGIN   ; time serie 
    96       ; get hovmoeller box + stride 
    97       box_plot = def_box(cmd.plt, dimplot, legbox, time_stride) 
    98       IF debug_w THEN print, '    box_plot for domdef in data_read = ', box_plot 
    99       CASE n_elements(box_plot) OF 
    100          4: domdef, box_plot[0:3], /MEMEINDICES 
    101          6: domdef, box_plot[0:5], /MEMEINDICES 
    102       ENDCASE  
    103       ; time grid for hovmoeller 
    104       time1 = 1+delta_t1 
    105       IF cmd.out EQ 'cdf' THEN nb_cycles = 1 
    106       timearr = compute_time(cmd.timave, cmd.date1, cmd.spec) 
    107       jpt = timearr.count 
    108       time2 = jpt+delta_t1 
    109       time = timearr.scale 
    110       niveau = 1 
    111 ;      print, '  Data_read, hotyp/jpt = ', hotyp, jpt 
    112 ;      print, '  Data_read, box = ', box_plot[0:3] 
    113    ENDIF ELSE BEGIN             ; single date 
    114       jpt = 1 
    115       CASE strmid(cmd.timave, strlen(cmd.timave)-1, 1) OF 
    116          'y': time1 = delta_t1+1 
    117          'm': BEGIN 
    118             CASE strmid(cmd.timave, strlen(cmd.timave)-2, 1) OF 
    119                'm': time1 = delta_t1+long(strmid(cmd.date1, 0, 2))  ; mean month 
    120                ELSE: time1 = delta_t1+1 
    121             ENDCASE  
    122             END  
    123          'd': time1 = delta_t1+1 
    124          ELSE: time1 = delta_t1+1 
    125       ENDCASE  
    126       time2 = time1 
    127       time = time1 
    128       print, '    Setting domdef in data_read, plt typ =', plttyp    
    129       CASE plttyp OF 
    130          'plt': BEGIN 
    131             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
    132             IF strmid(cmd.var, 0,2) EQ 'vo' THEN BEGIN ; get right zindex for vertical 
    133                CASE vert_type OF 
    134                   'level':domdef,[ box_plot[0:3],vert_mean],/zindex, /MEMEINDICES 
    135                   'z':domdef,[ box_plot[0:3],vert_mean], /MEMEINDICES 
    136                   ELSE: domdef, box_plot[0:3], /MEMEINDICES 
    137                ENDCASE  
    138             ENDIF ELSE BEGIN 
    139                domdef, box_plot[0:3], /MEMEINDICES 
    140             ENDELSE 
    141          END  
    142          'pltz': BEGIN 
    143             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
    144             domdef, box_plot[0:3], /MEMEINDICES 
    145          END  
    146          'plt1d': BEGIN 
    147             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
    148             domdef, box_plot[0:3], /MEMEINDICES 
    149             ; force all_data read 
    150             all_data = 1 
    151          END  
    152          ELSE: domdef, /MEMEINDICES 
    153       ENDCASE  
    154    ENDELSE  
    155  
    156 ;   print, jpt, time1, time2, hotyp 
     93; define horizontal domain 
     94; if needed, vertical domain will be defined in the call of read_ncdf 
     95   box_plot = def_box(cmd.plt, dimplot, legbox, time_stride) 
     96   ;domdef, [box_plot[0:3], vert_mean], /MEMEINDICES 
     97   ;domdef, box_plot[0:3], /MEMEINDICES 
     98   CASE n_elements(box_plot) OF 
     99      4 : IF vert_switch GE 1 THEN box_plot = [box_plot, vert_mean] 
     100      6 : IF vert_switch GE 1 THEN box_plot = [box_plot[0:3], vert_mean] 
     101   ENDCASE 
     102 
     103; define timetsteps time1 and time2 
     104    IF hotyp NE '-' THEN BEGIN 
     105       time1 = delta_t1 
     106       timearr = compute_time(cmd.timave, cmd.date1, cmd.spec) 
     107       time2 =  timearr.count+delta_t1 
     108    ENDIF ELSE BEGIN 
     109       IF strpos(cmd.timave, 'mm') NE -1 THEN time1 = delta_t1+long(strmid(cmd.date1, 0, 2)) $ 
     110       ELSE  time1 = delta_t1 
     111       time2 = time1 
     112    ENDELSE 
     113 
     114;    IF hotyp NE '-' THEN BEGIN   ; time serie 
     115;       ; get hovmoeller box + stride 
     116;       box_plot = def_box(cmd.plt, dimplot, legbox, time_stride) 
     117;       IF debug_w THEN print, '    box_plot for domdef in data_read = ', box_plot 
     118;       CASE n_elements(box_plot) OF 
     119;          4: domdef, box_plot[0:3], /MEMEINDICES 
     120;          6: domdef, box_plot[0:5], /MEMEINDICES 
     121;       ENDCASE  
     122;       ; time grid for hovmoeller 
     123;       time1 = 1+delta_t1 
     124;       IF cmd.out EQ 'cdf' THEN nb_cycles = 1 
     125;       timearr = compute_time(cmd.timave, cmd.date1, cmd.spec) 
     126;       jpt = timearr.count 
     127;       time2 = jpt+delta_t1 
     128;       time = timearr.scale 
     129;       niveau = 1 
     130; ;      print, '  Data_read, hotyp/jpt = ', hotyp, jpt 
     131; ;      print, '  Data_read, box = ', box_plot[0:3] 
     132;    ENDIF ELSE BEGIN             ; single date 
     133;       jpt = 1 
     134;       CASE strmid(cmd.timave, strlen(cmd.timave)-1, 1) OF 
     135;          'y': time1 = delta_t1+1 
     136;         'm': BEGIN 
     137;             CASE strmid(cmd.timave, strlen(cmd.timave)-2, 1) OF 
     138;                'm': time1 = delta_t1+long(strmid(cmd.date1, 0, 2))  ; mean month 
     139;                ELSE: time1 = delta_t1+1 
     140;             ENDCASE  
     141;             END  
     142;          'd': time1 = delta_t1+1 
     143;          ELSE: time1 = delta_t1+1 
     144;       ENDCASE  
     145;       time2 = time1 
     146;       time = time1 
     147;       print, '    Setting domdef in data_read, plt typ =', plttyp    
     148;       CASE plttyp OF 
     149;          'plt': BEGIN 
     150;             box_plot = def_box(cmd.plt, 2, legbox, time_stride) 
     151;             IF strmid(cmd.var, 0,2) EQ 'vo' THEN BEGIN ; get right zindex for vertical 
     152;                CASE vert_type OF 
     153;                   'level':domdef,[ box_plot[0:3],vert_mean],/zindex, /MEMEINDICES 
     154;                   'z':domdef,[ box_plot[0:3],vert_mean], /MEMEINDICES 
     155;                   ELSE: domdef, box_plot[0:3], /MEMEINDICES 
     156;                ENDCASE  
     157;             ENDIF ELSE BEGIN 
     158;                domdef, box_plot[0:3], /MEMEINDICES 
     159;             ENDELSE 
     160;          END  
     161;          'pltz': BEGIN 
     162;             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
     163;             domdef, box_plot[0:3], /MEMEINDICES 
     164;          END  
     165;          'plt1d': BEGIN 
     166;             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
     167;             domdef, box_plot[0:3], /MEMEINDICES 
     168;             ; force all_data read 
     169;             all_data = 1 
     170;          END  
     171;          ELSE: domdef, /MEMEINDICES 
     172;      ENDCASE  
     173;    ENDELSE  
     174;    ;box_plot = [ box_plot[0:3],vert_mean] 
     175;    print, jpt, time1, time2, hotyp 
    157176 
    158177; save time in common fld_att 
    … …  
    166185      '@@': BEGIN 
    167186         IF debug_w THEN print,  'keyword_set(all_data) : ', keyword_set(all_data) 
    168          fldr = macro_read(file_name, cmd.var, ncdf_db, TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
     187         fldr = macro_read(file_name, cmd.var, ncdf_db, BOXZOOM = box_plot, TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
    169188         IF stddev_diff EQ 1 THEN fldr.origin =  'diff' 
    170189      END  
    … …  
    294313 
    295314               ENDIF ELSE BEGIN  ; general case 
    296                   fldr = nc_read(file_name, cmd.var, ncdf_db, $ 
    297                                  TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
    298  
     315                  fldr = nc_read(file_name, cmd.var, ncdf_db, BOXZOOM = box_plot, $ 
     316                                  TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
     317                  ;fldr = nc_read(file_name, cmd.var, ncdf_db, $ 
     318                  ;               TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)               
    299319                  ; perform mean if hovmoeller diagram is required 
    300320                  ; It is absolutely necessary for the division case with  
    … …  
    315335      END   
    316336   ENDCASE   
     337 
     338; Pb with varexp which is always updated in read_ncdf 
     339; For Seb, varexp is the name of the file 
     340; For Eric, varexp is the name of the experiment 
     341   varexp = cmd.exp 
    317342 
    318343; if tkeavt, take log 

  • branches/procs-read/macro_read.pro

    r48 r81  
    22; Read Macro field Data 
    33;---------------------------   
    4 FUNCTION macro_read, file_name, var_name, ncdf_db, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA =  all_data, _extra = ex 
     4FUNCTION macro_read, file_name, var_name, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA =  all_data, _extra = ex 
    55@common 
    66@com_eg 
    … …  
    3434            ELSE:  
    3535         ENDCASE  
    36          strcall = 'field_lec = '+macro_def[2]+'(file_name, ncdf_db, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, _extra = ex)' 
     36         strcall = 'field_lec = '+macro_def[2]+'(file_name, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, _extra = ex)' 
    3737         res = execute(strcall) 
    3838         fieldo = {name: '', data: field_lec.data, legend: '', units: '', origin: '', dim: 0, direc:''} 

  • branches/procs-read/macros/make_eos.pro

    r2 r81  
    4949;------------------------------------------------------------ 
    5050;------------------------------------------------------------ 
    51 FUNCTION make_eos, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ZMTYP = zmtyp 
     51FUNCTION make_eos, file_name, ncdf_db,  BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, ZMTYP = zmtyp 
    5252; 
    5353@common 
    … …  
    5757; Read T and S 
    5858; 
    59    tn = nc_read(file_name,'votemper', ncdf_db, TIME_1 = time_1,  TIME_2 =  time_2, no_mean = 1) 
    60    sn = nc_read(file_name,'vosaline', ncdf_db, TIME_1 = time_1,  TIME_2 =  time_2, no_mean = 1) 
     59   tn = nc_read(file_name,'votemper', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, no_mean = 1) 
     60   sn = nc_read(file_name,'vosaline', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, no_mean = 1) 
    6161 
    6262; declarations 

  • branches/procs-read/macros/make_msf.pro

    r2 r81  
    22; make MSF 
    33; 
    4 FUNCTION make_msf, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ZMTYP = zmtyp 
     4FUNCTION make_msf, file_name, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, ZMTYP = zmtyp 
    55 
    66@common 
    … …  
    1414; 
    1515   file_nam = strmid(file_name, 0, strlen(file_name)-4) 
    16    v = nc_read(file_nam+'V.nc','vomecrty', ncdf_db, TIME_1 = time_1,  TIME_2 =  time_2) 
     16   v = nc_read(file_nam+'V.nc','vomecrty', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2) 
    1717 
    1818   idx = where(v.data EQ valmask) 

  • branches/procs-read/macros/make_stddev.pro

    r17 r81  
    22; make std dev from 2D macro_base_fld monthly time serie (defined in macro_read) 
    33; 
    4 FUNCTION make_stddev, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, ZMTYP = zmtyp 
     4FUNCTION make_stddev, file_name, ncdf_db,  BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, ZMTYP = zmtyp 
    55 
    66@common 
    … …  
    1212   IF debug_w THEN print, 'keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA)  
    1313 
    14    mfld = nc_read(file_name, macro_base_fld, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data) 
    15  
     14   mfld = nc_read(file_name, macro_base_fld, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data) 
    1615   mfld.data = trends(mfld.data, '412', 'xyt') 
    1716 

  • branches/procs-read/meshes/mesh_from_file.pro

    r69 r81  
    1515   print,'    Model inits for ', model, ' from file: ', s_file 
    1616 
    17    sm_file = hom_idl+'grids/grids_'+model+'.nc'  
    18    res = find(sm_file) 
     17;   sm_file = hom_idl+'grids/grids_'+model+'.nc'  
     18;   res = find(sm_file) 
    1919   varexpp = varexp 
    20    IF res NE 'NOT FOUND' THEN BEGIN  
    21  
    22       initncdf, s_file, GLAMBOUNDARY = glamboundary_box 
    23  
    24       print, '    Found mask from ',sm_file  
    25  
    26       tmask = byte(read_ncdf('sftlf', 0, 0, /timestep, file = sm_file, /nostruct)) 
    27        
    28       idx = where(tmask EQ valmask) 
    29  
    30       IF idx(0) NE -1 THEN tmask(idx) = 0. 
    31       idx =  where(tmask LE 50.) 
    32       tmask(idx) = 0. 
    33       tmask =  tmask < 1 
    34       tmask =  1-tmask 
    35       triangles=triangule() 
    36  
    37    ENDIF ELSE BEGIN  
     20;   IF res NE 'NOT FOUND' THEN BEGIN  
     21; 
     22;      initncdf, s_file, GLAMBOUNDARY = glamboundary_box 
     23; 
     24;      print, '    Found mask from ',sm_file  
     25; 
     26;      tmask = byte(read_ncdf('sftlf', 0, 0, /timestep, file = sm_file, /nostruct)) 
     27;       
     28;      idx = where(tmask EQ valmask) 
     29; 
     30;      IF idx(0) NE -1 THEN tmask(idx) = 0. 
     31;      idx =  where(tmask LE 50.) 
     32;      tmask(idx) = 0. 
     33;      tmask =  tmask < 1 
     34;      tmask =  1-tmask 
     35;      triangles=triangule() 
     36; 
     37;   ENDIF ELSE BEGIN  
    3838      CASE var_name OF 
    3939         '@@voenergy': var_local = 'so' 
    … …  
    6363      ; build grid 
    6464      CASE mesh_type of 
    65          'atm': initncdf, s_file, GLAMBOUNDARY = glamboundary_box, ZAXISNAME = name_level 
     65         'atm': BEGIN 
     66            ;;initncdf, s_file, USEASMASK = var_local, missing_value = valmask, GLAMBOUNDARY = glamboundary_box, ZAXISNAME = name_level 
     67            initncdf, s_file, GLAMBOUNDARY = glamboundary_box, ZAXISNAME = name_level 
     68         END 
    6669         ELSE: initncdf, s_file, USEASMASK = var_local, missing_value = valmask, GLAMBOUNDARY = glamboundary_box, ZAXISNAME = 'depth' 
    6770      ENDCASE  
    … …  
    7174      IF debug_w THEN print,  '     gdept :', size(gdept) 
    7275       
    73    ENDELSE  
     76;   ENDELSE  
    7477 
    7578   varexp = varexpp 
    … …  
    8487   IF debug_w THEN print, '  ' 
    8588 
     89   stop 
     90 
    8691  return 
    8792END  

  • branches/procs-read/meshes/mesh_orca.pro

    r26 r81  
    5454 
    5555 
    56    IF keyword_set(NO_SHIFT) THEN key_shift = 0 
    57  
     56;   IF keyword_set(NO_SHIFT) THEN key_shift = 0 
     57   IF keyword_set(NO_SHIFT) THEN shift = 0 ;; key_shift = 0 
    5858   no_lon_shift = 0 
    5959 
    60    key_periodique = 1 
     60;   key_periodique = 1 
     61   periodic = 1 
    6162 
    6263; reduce grid in latitude ? 
    … …  
    120121      'ORCA05':  BEGIN 
    121122         mesmsk = 'micromeshmaskORCA05.nc' 
    122          cmd_grid =  'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box' 
     123         ;;cmd_grid =  'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box' 
     124         cmd_grid =  'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box, PERIODIC = periodic' 
    123125      END 
    124126      'ORCA_R2': BEGIN 
    … …  
    130132                  'V3': mesmsk = 'meshmask_ORCA_R2_V2.nc' 
    131133               ENDCASE 
    132                IF whole_arrays EQ 0 THEN BEGIN  
    133                    cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box, onearth = onearth' 
     134               IF whole_arrays EQ 0 THEN BEGIN 
     135                   
     136                  cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box, PERIODIC=periodic' 
     137                  ;;cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = 
     138                  ;;glamboundary_box, onearth = onearth' 
    134139               ENDIF ELSE BEGIN  
    135                    cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box' 
     140                  ;;cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box 
     141                  cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box, SHIFT=shift, PERIODIC=periodic' 
    136142               ENDELSE  
    137143            END  
    138144            'L300': BEGIN 
    139145               mesmsk = 'meshmask.orca.2d.L300.nc' 
    140                cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box' 
     146               ;;cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box' 
     147               cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box, PERIODIC=periodic' 
    141148            END  
    142149         ENDCASE  
    … …  
    144151      'ORCA_R4': BEGIN 
    145152         mesmsk = 'meshmask_orca4.nc' 
    146          cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box'  
     153         ;;cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box' 
     154         cmd_grid = 'ncdf_meshread, sm_file, GLAMBOUNDARY = glamboundary_box, PERIODIC=periodic'  
    147155      END  
    148156   ENDCASE 
    … …  
    153161   sm_file = hom_idl+'grids/'+mesmsk 
    154162   res = find(sm_file) 
    155  
     163    
    156164   IF res NE 'NOT FOUND' THEN BEGIN 
    157        res_grid =  execute(cmd_grid) 
     165      res_grid =  execute(cmd_grid) 
    158166   ENDIF ELSE BEGIN 
    159        stop, 'No meshmask file found for ORCA mesh config' 
     167      stop, 'No meshmask file found for ORCA mesh config' 
    160168   ENDELSE 
    161169 
    … …  
    243251   print, '    End of ORCA mesh config ' 
    244252 
     253 
    245254END 

  • branches/procs-read/nc_read.pro

    r41 r81  
    33; EG 24/2/99 
    44;---------------------------   
    5 FUNCTION nc_read, file_name, var_name, ncdf_db, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, NO_MEAN = no_mean 
     5FUNCTION nc_read, file_name, var_name, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, NO_MEAN = no_mean 
    66 
    77; arguments = file_name, varname 
    … …  
    1111@common 
    1212@com_eg 
    13  
    14 ;;   print,  'key_yreverse : ', key_yreverse   
    15  
     13   
    1614; inits 
    17    IF debug_w THEN print, ' ' 
    18    IF debug_w THEN print, '  ENTER nc_read' 
     15   IF debug_w THEN print, '   ' 
     16   IF debug_w THEN print, '   ENTER nc_read...' 
    1917 
    2018   IF NOT keyword_set(TIME_1) THEN time_1 =  1 
    2119   IF NOT keyword_set(TIME_2) THEN time_2 =  time_1 
    22 ; 
     20 
     21 
    2322; decide which subdomain of data to read 
    24 ; 
    25  
    26    IF debug_w THEN print, '     keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA)  
    27  
     23   IF debug_w THEN print, '     keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA) 
    2824   IF keyword_set(ALL_DATA) THEN BEGIN 
    29       print, '    Reading whole domain' 
    30       firstx = 0 
    31       firsty = 0 
    32       firstz = 0 
    33       nx = jpi 
    34       ny = jpj 
    35       nz = jpk 
    36       lastx = jpi-1 
    37       lasty = jpj-1 
    38       lastz = jpk-1 
     25      tout = 1 
     26   ENDIF ELSE tout = 0 
     27   CASE vert_type OF 
     28      'z' : zindex = 0 
     29      'level' : zindex = 1 
     30      ELSE: zindex = 0 
     31   ENDCASE 
     32 
     33;   IF keyword_set(ALL_DATA) THEN BEGIN 
     34;      print, '    Reading whole domain' 
     35;      firstx = 0 
     36;      firsty = 0 
     37;      firstz = 0 
     38;      nx = jpi 
     39;      ny = jpj 
     40;      nz = jpk 
     41;      lastx = jpi-1 
     42;      lasty = jpj-1 
     43;      lastz = jpk-1 
    3944;Rajout MK 
    4045;Trouver une meilleure place 
    41       ixminmesh = 0 
    42       iyminmesh = 0 
    43       ixmaxmesh = jpi-1 
    44       iymaxmesh = jpj-1 
     46;      ixminmesh = 0 
     47;      iyminmesh = 0 
     48;      ixmaxmesh = jpi-1 
     49;      iymaxmesh = jpj-1 
    4550;Fin rajout 
    46    ENDIF ELSE BEGIN 
    47       grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
    48       mask = 1 
    49    ENDELSE 
    50  
    51    IF debug_w THEN print, '     nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz' 
    52    IF debug_w THEN print, '     ', nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
    53  
    54    IF debug_w THEN print,  '     key_offset =', key_offset 
    55    IF debug_w THEN print,  '     jpi, jpj, jpk =',  jpi, jpj, jpk 
     51;   ENDIF ELSE BEGIN 
     52;      grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
     53;      mask = 1 
     54;   ENDELSE 
     55;   IF debug_w THEN print, '     nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz' 
     56;   IF debug_w THEN print, '     ', nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
     57 
     58;   IF debug_w THEN print,  '     key_offset =', key_offset 
     59;   IF debug_w THEN print,  '     jpi, jpj, jpk, jpt =',  jpi, jpj, jpk, jpt 
    5660 
    5761 
    5862; define reading boundaries 
    59  
    60    x_count = nx 
    61    y_count = ny 
    62    z_count = nz 
    63  
    64    IF debug_w THEN print,  '     nx,ny,nz =', nx,ny,nz  
     63;   x_count = nx 
     64;   y_count = ny 
     65;   z_count = nz 
     66 
     67;   IF debug_w THEN print,  '     nx,ny,nz =', nx,ny,nz  
    6568 
    6669; force izmaxmesh to lastz 
    67    jpk = lastz+1 
     70;   jpk = lastz+1 
    6871 
    6972; dealing with longitude periodicity 
    70    IF x_count NE jpi THEN BEGIN 
    71       key_shift_read = 0 
    72    ENDIF ELSE key_shift_read = key_shift 
    73    IF debug_w THEN print,  '     key_shift_read=', key_shift_read 
    74  
    75    x_offset = key_offset[0]+ firstx+(key_shift_read-key_shift) 
    76    y_offset = key_offset[1]+ firsty 
    77    z_offset = key_offset[2]+ firstz 
    78    IF debug_w THEN print,  '     offset =', x_offset, y_offset, z_offset 
    79  
    80 ; pour trouver un fichier, tu as 
    81 ; findfile (tres pratique) et dialog_pickfile (interactif) 
    82  
    83 ; pour verifier si il y a une variable je fais 
    84  
    85 ; inq=ncdf_inquire(cdfid) 
    86 ; for varid=0,inq.nvars-1 do BEGIN 
    87 ;    varinq=ncdf_varinq(cdfid,varid)  
    88 ;    if varinq.name eq nom then goto, variabletrouvee 
    89 ; endfor 
    90 ; return, -1 
    91 ; variabletrouvee: 
     73;   IF x_count NE jpi THEN BEGIN 
     74;      key_shift_read = 0 
     75;   ENDIF ELSE key_shift_read = key_shift 
     76;   IF debug_w THEN print,  '     key_shift_read=', key_shift_read 
     77 
     78;   x_offset = key_offset[0]+ firstx+(key_shift_read-key_shift) 
     79;   y_offset = key_offset[1]+ firsty 
     80;   z_offset = key_offset[2]+ firstz 
     81;   IF debug_w THEN print,  '     offset =', x_offset, y_offset, z_offset 
    9282 
    9383; local directory 
    9484 
    95    IF strpos(ncdf_db, ':') GE 1 THEN directory = (str_sep(ncdf_db, ':'))[1] $ 
     85   IF strpos(ncdf_db, ':') GE 1 THEN directory = (strsplit(ncdf_db, ':', /EXTRACT))[1] $ 
    9686    ELSE directory = ncdf_db 
    9787 
    9888; existence of file 
    9989 
    100    list = findfile(directory+file_name, count = nb_file) 
    101     
    102    IF nb_file EQ 0 THEN BEGIN 
    103       print, '' 
    104       print, '   ***** file ', file_name, ' not found in ', $ 
    105        directory 
    106       print, '' 
    107       field = {data: -1.0} 
    108       return,  field 
    109    ENDIF  
     90;   list = findfile(directory+file_name, count = nb_file) 
     91    
     92;   IF nb_file EQ 0 THEN BEGIN 
     93;      print, '' 
     94;      print, '   ***** file ', file_name, ' not found in ', $ 
     95;       directory 
     96;      print, '' 
     97;      field = {data: -1.0} 
     98;      return,  field 
     99;   ENDIF  
    110100 
    111101; ouverture fichier netCDF + contenu 
    112    cdfid=ncdf_open(directory+file_name)  
    113    inq=ncdf_inquire(cdfid) 
     102;   cdfid=ncdf_open(directory+file_name)  
     103;   inq=ncdf_inquire(cdfid) 
    114104; que contient la variable 
    115    varid = ncdf_varid(cdfid, var_name) 
     105;   varid = ncdf_varid(cdfid, var_name) 
    116106 
    117107   name_suff = '' 
    118  
    119    IF varid EQ -1 THEN BEGIN 
    120       print, '' 
    121       print, '   ***** field ', var_name, ' not found in file ', $ 
    122        file_name 
    123       print, '' 
    124       field = {data: -1.0} 
    125       return,  field 
    126    ENDIF  
    127    varinq=ncdf_varinq(cdfid, var_name) 
     108   data_direc = '' 
     109 
     110;   IF varid EQ -1 THEN BEGIN 
     111;      print, '' 
     112;      print, '   ***** field ', var_name, ' not found in file ', $ 
     113;       file_name 
     114;      print, '' 
     115;      field = {data: -1.0} 
     116;      return,  field 
     117;   ENDIF  
     118;   varinq=ncdf_varinq(cdfid, var_name) 
    128119 
    129120; test sur la dimension 
    130    err_mess = '' 
    131    field_dim = n_elements(varinq.dim) 
     121;   err_mess = '' 
     122;   field_dim = n_elements(varinq.dim) 
    132123 
    133124; get unlimited record variable 
    134    IF inq.recdim NE -1 THEN BEGIN 
    135       ncdf_diminq, cdfid,  inq.recdim,  name_time, nb_time 
    136       ;;ncdf_varget, cdfid, inq.recdim, time_array 
    137       ;;nb_time = (size(time_array))(1) 
    138    ENDIF ELSE BEGIN 
    139       ; Look for a potential record dimension  
    140       IF debug_w THEN print, '    Warning : no unlimited record in netCDF file' 
    141       ; Look for the names of all dimensions and all variables 
    142       list_dims = ncdf_listdims(cdfid) 
    143       list_vars = ncdf_listvars(cdfid) 
    144       ; Propose all the names and make the user choose the right one 
    145       ; for the record dimension 
    146       print, 'In the file ',  file_name,  ', the dimensions are named :' 
    147       print, list_dims 
    148       print, 'In the file ',  file_name,  ', the variables are named :' 
    149       print, list_vars 
    150       name_time = xquestion('Among the preceding names, which one could refer to the record dimension ?',  /chkwid) 
    151       IF name_time NE '-' THEN BEGIN 
    152          dimidt = ncdf_dimid(cdfid, name_time) 
    153          ncdf_diminq, cdfid,  dimidt,  name_time, nb_time 
    154          print, 'You chose ', name_time,  ' as a record dimension and its size is ',  nb_time 
    155          inq.recdim = dimidt 
    156       ENDIF ELSE BEGIN 
    157          print, 'No record dimension considered in the file' 
    158          nb_time =  0 
    159       ENDELSE 
    160  
    161       IF varinq.dim[varinq.ndims-1] NE dimidt THEN STOP,  $ 
    162        'Post_it cannot handle variables whose record dimension is not the last one' 
    163  
    164 ;      ncdf_diminq,cdfid,(n_elements(varinq.dim)-1), name_time, nb_time      
    165 ;      dimidl = ncdf_dimid(cdfid, name_time) 
    166 ;      ncdf_diminq,cdfid,dimidl, name_time, nb_time 
    167 ;      varidl = ncdf_varid(cdfid, name_time) 
    168 ;      ncdf_varget, cdfid, varidl, time_array 
    169  
    170    ENDELSE  
    171  
    172    IF jpt GT nb_time THEN BEGIN  
    173       jpt =  nb_time 
    174       print, 'Warning : the specified jpt does not correspond to the real time dimension in the file !' 
    175       print, 'New jpt =',  jpt 
    176    ENDIF 
    177    ; defs for @read_ncdf_varget 
    178  
    179 ;   key_stride = time_stride 
    180    key_stride = 1 
    181    if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] 
    182    key_stride = 1l > long(key_stride) 
    183    jpitotal = long(ixmaxmesh-ixminmesh+1) 
    184    key_shift = long(testvar(var = key_shift)) 
    185     
    186    key = long(key_shift MOD jpitotal) 
    187    if key LT 0 then key = key+jpitotal 
    188    ixmin = ixminmesh-ixmindta 
    189    iymin = iyminmesh-iymindta 
    190    firsttps = time_1-1 
    191    lasttps = time_2-1 
    192  
    193    name = varid 
    194  
     125;   IF inq.recdim NE -1 THEN BEGIN 
     126;      ncdf_diminq, cdfid,  inq.recdim,  name_time, nb_time 
     127;      ;;ncdf_varget, cdfid, inq.recdim, time_array 
     128;      ;;nb_time = (size(time_array))(1) 
     129;   ENDIF ELSE BEGIN 
     130;      ; Look for a potential record dimension  
     131;      IF debug_w THEN print, '    Warning : no unlimited record in netCDF file' 
     132;      ; Look for the names of all dimensions and all variables 
     133;      list_dims = ncdf_listdims(cdfid) 
     134;      list_vars = ncdf_listvars(cdfid) 
     135;      ; Propose all the names and make the user choose the right one 
     136;      ; for the record dimension 
     137;      print, 'In the file ',  file_name,  ', the dimensions are named :' 
     138;      print, list_dims 
     139;      print, 'In the file ',  file_name,  ', the variables are named :' 
     140;      print, list_vars 
     141;      name_time = xquestion('Among the preceding names, which one could refer to the record dimension ?',  /chkwid) 
     142;      IF name_time NE '-' THEN BEGIN 
     143;         dimidt = ncdf_dimid(cdfid, name_time) 
     144;         ncdf_diminq, cdfid,  dimidt,  name_time, nb_time 
     145;         print, 'You chose ', name_time,  ' as a record dimension and its size is ',  nb_time 
     146;         inq.recdim = dimidt 
     147;      ENDIF ELSE BEGIN 
     148;         print, 'No record dimension considered in the file' 
     149;         nb_time =  0 
     150;      ENDELSE 
     151 
     152;      IF varinq.dim[varinq.ndims-1] NE dimidt THEN STOP,  $ 
     153;       'Post_it cannot handle variables whose record dimension is not the last one' 
     154 
     155;;      ncdf_diminq,cdfid,(n_elements(varinq.dim)-1), name_time, nb_time      
     156;;      dimidl = ncdf_dimid(cdfid, name_time) 
     157;;      ncdf_diminq,cdfid,dimidl, name_time, nb_time 
     158;;      varidl = ncdf_varid(cdfid, name_time) 
     159;;      ncdf_varget, cdfid, varidl, time_array 
     160 
     161;   ENDELSE  
     162 
     163;   IF jpt GT nb_time THEN BEGIN  
     164;      jpt =  nb_time 
     165;      print, 'Warning : the specified jpt does not correspond to the real time dimension in the file !' 
     166;      print, 'New jpt =',  jpt 
     167;   ENDIF 
     168;   ; defs for @read_ncdf_varget 
     169 
     170;;   key_stride = time_stride 
     171;   key_stride = 1 
     172;   if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] 
     173;   key_stride = 1l > long(key_stride) 
     174;   jpitotal = long(ixmaxmesh-ixminmesh+1) 
     175;   key_shift = long(testvar(var = key_shift)) 
     176    
     177;   key = long(key_shift MOD jpitotal) 
     178;   if key LT 0 then key = key+jpitotal 
     179;   ixmin = ixminmesh-ixmindta 
     180;   iymin = iyminmesh-iymindta 
     181;   firsttps = time_1-1 
     182;   lasttps = time_2-1 
     183 
     184;   name = varid 
     185 
     186;   IF debug_w THEN print, '     ' 
     187;   IF debug_w THEN print, '     key=', key  
     188;   IF debug_w THEN print, '     jpitotal=', jpitotal 
     189;   IF debug_w THEN print, '     ixmindta,iymindta,izmindta  =', ixmindta,iymindta,izmindta 
     190;   IF debug_w THEN print, '     ixminmesh,iyminmesh=', ixminmesh,iyminmesh 
     191;   IF debug_w THEN print, '     ixmaxmesh,iymaxmesh=', ixmaxmesh,iymaxmesh  
     192;   IF debug_w THEN print, '     izminmesh,izmaxmesh=', izminmesh,izmaxmesh  
     193;   IF debug_w THEN print, '     ixmin,iymin=', ixmin,iymin  
     194;   IF debug_w THEN print, '     firsttps,lasttps=', firsttps, lasttps 
     195;   IF debug_w THEN print, '     key_shift in nc_read=', key_shift 
     196;   IF debug_w THEN print, '     key_yreverse, firsty, lasty',key_yreverse, firsty, lasty 
     197;   IF debug_w THEN print, '     ' 
     198 
     199   ncdf_getatt, directory+file_name, var_name, add_offset = add_offset, scale_factor = scale_factor, $ 
     200    missing_value = missing_value, units = units, calendar = calendar, long_name = long_name 
     201 
     202; By default units for the Z axis are meters 
     203   IF n_elements(gdept) GT 1 AND name_level NE '-' THEN BEGIN 
     204     ncdf_getatt, directory+file_name, name_level, units = units_depth 
     205   ENDIF ELSE units_depth = 'm' 
     206 
     207; Read the data with a call to read_ncdf 
    195208   IF debug_w THEN print, '     ' 
    196    IF debug_w THEN print, '     key=', key  
    197    IF debug_w THEN print, '     jpitotal=', jpitotal 
     209   IF debug_w THEN print, '     Reading raw data from ',  file_name 
     210 
     211   lec_data = read_ncdf(var_name, time_1-1, time_2-1, BOXZOOM = boxzoom, FILENAME = directory+file_name, $ 
     212                        /TIMESTEP, /ADDSCL_BEFORE, TOUT = tout, /NOSTRUCT, /CONT_NOFILL, GRID = vargrid, $ 
     213                        ZINVAR = zinvar, ZINDEX = zindex) 
     214 
     215   print,  'zinvar', zinvar 
     216    
     217   IF debug_w THEN print, '     ' 
     218;   IF debug_w THEN print, '     key=', key  
     219;   IF debug_w THEN print, '     jpitotal=', jpitotal 
    198220   IF debug_w THEN print, '     ixmindta,iymindta,izmindta  =', ixmindta,iymindta,izmindta 
    199221   IF debug_w THEN print, '     ixminmesh,iyminmesh=', ixminmesh,iyminmesh 
    200222   IF debug_w THEN print, '     ixmaxmesh,iymaxmesh=', ixmaxmesh,iymaxmesh  
    201223   IF debug_w THEN print, '     izminmesh,izmaxmesh=', izminmesh,izmaxmesh  
    202    IF debug_w THEN print, '     ixmin,iymin=', ixmin,iymin  
    203    IF debug_w THEN print, '     firsttps,lasttps=', firsttps, lasttps 
     224;   IF debug_w THEN print, '     ixmin,iymin=', ixmin,iymin  
     225;   IF debug_w THEN print, '     firsttps,lasttps=', firsttps, lasttps 
     226   IF debug_w THEN print, '     jpt=', jpt 
    204227   IF debug_w THEN print, '     key_shift in nc_read=', key_shift 
    205    IF debug_w THEN print, '     key_yreverse, firsty, lasty',key_yreverse, firsty, lasty 
     228;   IF debug_w THEN print, '     key_yreverse, firsty, lasty',key_yreverse, firsty, lasty 
     229   IF debug_w THEN print, '     key_yreverse',key_yreverse  
    206230   IF debug_w THEN print, '     ' 
    207231 
    208    CASE n_elements(varinq.dim) OF 
    209  
    210       ;; fichier 2d : surface 
    211       2: BEGIN 
    212          print, '    Reading ', var_name, ' (2D data)   from ', file_name 
    213          @read_ncdf_varget 
    214          lec_data = res 
    215          data_direc = 'xy' 
    216          niveau = 1 
    217          END  
    218  
    219       ;; fichier 3d : 2 cas = 1/ 2d espace + temps 2/ 3d espace 
    220       3: BEGIN  
    221          ;; si varinq.dim contient la dim infinie (no 3) -> temps 
    222          dim_3 = '3d' 
    223          IF nb_time GE 1 THEN dim_3 = '2d' 
    224          IF dim_3 EQ '2d' THEN BEGIN  
    225           
    226             IF time_2 EQ time_1 THEN BEGIN  
    227                print, '    Reading ', var_name, ' (2D data - index ', strcompress(string(time_1)),') from ', file_name 
    228                @read_ncdf_varget 
    229                lec_data = res 
    230                data_direc = 'xy' 
    231             ENDIF ELSE BEGIN  
    232                print, '    Reading ', var_name, ' (2D data time serie)', strcompress(string(time_1)),'-', strtrim(string(time_2), 2), ' [every ',strtrim(string(time_stride), 2), ']   from ', file_name 
    233                IF debug_w THEN print,  '     x_offset et y_offset et time_1 :',  x_offset,  y_offset,  time_1 
    234                @read_ncdf_varget 
    235                lec_data = res 
    236                data_direc = 'xyt'  
    237             ENDELSE  
    238  
    239          ENDIF ELSE BEGIN  
    240             print, '    Reading ', var_name, ' (3D data)', '   from ', file_name 
    241             @read_ncdf_varget 
    242             lec_data = res 
    243             data_direc = 'xyz' 
    244          ENDELSE 
    245       END  
    246       ;; fichier 4d : volume + temps 
    247       4: BEGIN 
    248          IF time_2 EQ time_1 THEN BEGIN  
    249             print, '    Reading ', var_name, ' (3D data - index ', strcompress(string(time_1)),') from ', file_name 
    250             ; read vertical levels 
    251             IF debug_w THEN print, '   mesh_type= ',mesh_type  
    252             IF mesh_type EQ 'atm' THEN BEGIN 
    253                ncdf_diminq,cdfid,2, name_level, nb_level 
    254                IF debug_w THEN print, 'name_level, nb_level = ', name_level, nb_level 
    255  
    256                dimidl = ncdf_dimid(cdfid, name_level) 
    257                ncdf_diminq,cdfid,dimidl, name_level, nb_level 
    258  
    259                varidl = ncdf_varid(cdfid, name_level) 
    260                ; make sure name_level is in hPa 
    261                ncdf_attget, cdfid, varidl, 'units', val_unit 
    262                IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar'  AND string(val_unit) NE 'mbar' THEN BEGIN  
    263                   print, '  vertical levels coordinate not obvious = ', name_level 
    264                   print, '  ... using <levels> ...' 
    265                   varidl = ncdf_varid(cdfid, 'levels') 
    266                ENDIF  
     232   ;lec_data = struct_data.tab 
     233   ;data_direc = 'xy' 
     234   ;field_dim = (size(struct.tab))[0] 
     235   field_dim = (size(lec_data))[0] 
     236   IF field_dim LE 0 THEN stop,  '  Something wrong happened in read_ncdf ' 
     237 
     238   ; build 3D atmospheric mask if needed and mask values 
     239;   IF mesh_type EQ 'atm' THEN BEGIN 
     240;      ;tmask = make_3d_atm_msk(directory+file_name, field_dim) 
     241;      CASE atmos_msk OF  
     242;         0: print, '          [atmos grid : take all points] ' ; take all points 
     243;         1: BEGIN 
     244;                                ; take ocean points only 
     245;            idx = where(tmask EQ 0) 
     246;            stop 
     247;            IF idx NE -1 THEN lec_data(idx) = 1.e20 
     248;            print, '          [atmos grid : take ocean points only] ' 
     249;         END 
     250;         2: BEGIN 
     251;                                ; take land points only 
     252;            idx = where(tmask GT 0) 
     253;            IF idx NE -1 THEN lec_data(idx) = 1.e20 
     254;            print, '          [atmos grid : take land points only] ' 
     255;         END 
     256;      ENDCASE  
     257;   ENDIF 
     258    
     259; Average data along vertical if needed and update some features 
     260; needed for plt (data_direc, name_suff)    
     261   update_data, TAB = lec_data, VNAME = var_name, BOXZOOM = boxzoom, ZUNITS = units_depth, $ 
     262    ZINVAR = zinvar, SUFFIX = name_suff, D_DIREC = data_direc, $ 
     263    TIME_1 = time_1, TIME_2 =  time_2, NO_MEAN = no_mean 
     264   field_dim = (size(lec_data))[0] 
     265   ;stop 
     266 
     267;   CASE n_elements(varinq.dim) OF 
     268;   CASE no_read_ncdf_varget OF 
     269;   CASE field_dim OF 
     270;      1: BEGIN 
     271;         print,  ' !!!! 1D case not implemented !!!!' 
     272;         print,  ' !!!! read_ncdf cannot handle 1D data !!!!' 
     273;         END 
     274; 
     275;      ;; fichier 2d : surface 
     276;      2: BEGIN 
     277;         print, '    Found ', var_name, ' (2D data)   from ', file_name 
     278;        ;@read_ncdf_varget 
     279;         ;lec_data = res         
     280;         data_direc = 'xy' 
     281;         niveau = 1 
     282;         END  
     283 
     284;     ;; fichier 3d : 2 cas = 1/ 2d espace + temps 2/ 3d espace 
     285;     3: BEGIN  
     286;         ;; si varinq.dim contient la dim infinie (no 3) -> temps 
     287;         dim_3 = '3d' 
     288;         IF nb_time GE 1 THEN dim_3 = '2d' 
     289;         IF jpt GE 1 THEN dim_3 = '2d' 
     290;         IF dim_3 EQ '2d' THEN BEGIN  
     291;         IF jpt EQ 1 THEN BEGIN 
     292; 
     293;            print, '    Found ', var_name, ' (3D data)', '   from ', file_name 
     294;            ;@read_ncdf_varget 
     295;            ;lec_data = res 
     296;            data_direc = 'xyz' 
     297; 
     298;            IF time_2 EQ time_1 THEN BEGIN  
     299;               print, '    Found ', var_name, ' (2D data - index ', strcompress(string(time_1)),') from ', file_name 
     300;               ;@read_ncdf_varget 
     301;               ;lec_data = res 
     302;               data_direc = 'xy' 
     303;            ENDIF ELSE BEGIN  
     304;               print, '    Found ', var_name, ' (2D data time serie)', strcompress(string(time_1)),'-', strtrim(string(time_2), 2), ' [every ',strtrim(string(time_stride), 2), ']   from ', file_name 
     305;               ;IF debug_w THEN print,  '     x_offset et y_offset et time_1 :',  x_offset,  y_offset,  time_1 
     306;               ;@read_ncdf_varget 
     307;               ;lec_data = res 
     308;               data_direc = 'xyt'  
     309;            ENDELSE  
     310; 
     311;         ENDIF ELSE IF jpt GT 1 THEN BEGIN 
     312;         ENDELSE 
     313;      END  
     314; 
     315;      ;; fichier 4d : volume + temps 
     316;      4: BEGIN 
     317; 
     318;         IF time_2 EQ time_1 THEN BEGIN  
     319;            print, '    Found ', var_name, ' (3D data - index ', strcompress(string(time_1)),') from ', file_name 
     320;            ; read vertical levels 
     321;            IF debug_w THEN print, '   mesh_type= ',mesh_type  
     322;            IF mesh_type EQ 'atm' THEN BEGIN 
     323;               ncdf_diminq,cdfid,2, name_level, nb_level 
     324;               IF debug_w THEN print, 'name_level, nb_level = ', name_level, nb_level 
     325; 
     326;               dimidl = ncdf_dimid(cdfid, name_level) 
     327;               ncdf_diminq,cdfid,dimidl, name_level, nb_level 
     328; 
     329;               varidl = ncdf_varid(cdfid, name_level) 
     330;               ; make sure name_level is in hPa 
     331;               ncdf_attget, cdfid, varidl, 'units', val_unit 
     332;               IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar'  AND string(val_unit) NE 'mbar' THEN BEGIN  
     333;                  print, '  vertical levels coordinate not obvious = ', name_level 
     334;                  print, '  ... using <levels> ...' 
     335;                  varidl = ncdf_varid(cdfid, 'levels') 
     336;               ENDIF  
    267337  
    268                ncdf_varget, cdfid, varidl, gdept  
    269                gdept = gdept 
    270                gdepw = gdept 
    271                e3t = shift(gdept, 1)-gdept 
    272                e3t[0] = e3t[1] 
    273                e3w = e3t 
    274                jpk = nb_level 
    275                tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) 
    276             ENDIF  
    277             IF debug_w THEN print,  ' going into @read_ncdf_varget '  
    278             @read_ncdf_varget 
    279             lec_data = res 
     338;               ncdf_varget, cdfid, varidl, gdept  
     339;               gdept = gdept 
     340;               gdepw = gdept 
     341;               e3t = shift(gdept, 1)-gdept 
     342;               e3t[0] = e3t[1] 
     343;               e3w = e3t 
     344;               jpk = nb_level 
     345;               tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) 
     346;            ENDIF  
     347;            IF debug_w THEN print,  ' going into @read_ncdf_varget '  
     348;            @read_ncdf_varget 
     349;            lec_data = res 
    280350;            lec_data =  reform(lec_data,x_count, y_count, z_count)  
    281             data_direc = 'xyz' 
    282             IF debug_w THEN help, lec_data 
     351;            data_direc = 'xyz' 
     352;            IF debug_w THEN help, lec_data 
    283353 
    284354;           vertical average ? 
    285             IF vert_switch ge 1 AND NOT keyword_set(no_mean) THEN BEGIN  
    286                old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    287                print, '      Average in vertical domain ', vert_type, vert_mean 
    288                IF mesh_type EQ 'atm' THEN BEGIN 
    289                   CASE atmos_msk OF  
    290                      0: print, '          [atmos grid : take all points] ' ; take all points 
    291                      1: BEGIN 
    292                         ; take ocean points only 
    293                         idx = where(tmask EQ 0) 
    294                         lec_data(idx) = 1.e20 
    295                         print, '          [atmos grid : take ocean points only] ' 
    296                      END 
    297                      2: BEGIN 
    298                         ; take land points only 
    299                         idx = where(tmask GT 0) 
    300                         lec_data(idx) = 1.e20 
    301                         print, '          [atmos grid : take land points only] ' 
    302                      END 
    303                   ENDCASE  
    304                   CASE vert_type OF 
    305                      'z': BEGIN ;  depth range 
    306                         IF debug_w THEN print, ' performing average' 
    307                         zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
    308                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    309                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
    310                         ENDIF ELSE BEGIN 
    311                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    312                         ENDELSE  
    313                      END  
    314                      ELSE: BEGIN ; levels range 
    315                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    316                            zmean = lec_data(*, *, vert_mean[0]) 
    317                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
    318                         ENDIF ELSE BEGIN 
    319                            zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) 
    320                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    321                         ENDELSE  
    322                      END  
    323                   ENDCASE  
    324 ;                  tmask = reform(tmask(*, *, 0), jpi, jpj) 
    325                ENDIF ELSE BEGIN 
    326                   ; ocean = always mask 
    327 ;                  idx = where(tmask EQ 0) 
    328 ;                  lec_data(idx) = 1.e20 
    329                   CASE vert_type OF 
    330                      'z': BEGIN 
    331                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    332                            zmean = lec_data ; right depth already selected 
    333                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    334                         ENDIF ELSE BEGIN 
    335                            zmean = moyenne(lec_data, 'z', boite = vert_mean) 
    336                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    337                         ENDELSE  
    338                      END  
    339                      ELSE: BEGIN ; case level (zindex) 
    340                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    341                            zmean = lec_data ; right depth already selected 
    342                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    343                         ENDIF ELSE BEGIN 
    344                            zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) 
    345                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    346                         ENDELSE  
    347                      END  
    348                   ENDCASE  
    349                ENDELSE  
    350                lec_data = zmean 
    351                domdef, old_boite 
    352                field_dim = field_dim - 1 
    353                data_direc = 'xy' 
    354                nzt = 1 
    355                firstz = 0 
    356                lastz = 0 
    357             ENDIF  
    358  
    359          ENDIF ELSE BEGIN  
    360             print, '    Reading ', var_name, ' (3D data time serie)',  strcompress(string(time_1)),'-', strtrim(string(time_2), 2), '   from ', file_name 
    361             ; read vertical levels 
    362             IF mesh_type EQ 'atm' THEN BEGIN 
    363                ncdf_diminq,cdfid,2, name_level, nb_level 
     355;            IF vert_switch ge 1 AND NOT keyword_set(no_mean) THEN BEGIN  
     356;               old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
     357;               print, '      Average in vertical domain ', vert_type, vert_mean 
     358;               IF mesh_type EQ 'atm' THEN BEGIN 
     359;                  CASE atmos_msk OF  
     360;                     0: print, '          [atmos grid : take all points] ' ; take all points 
     361;                     1: BEGIN 
     362;                        ; take ocean points only 
     363;                        idx = where(tmask EQ 0) 
     364;                        lec_data(idx) = 1.e20 
     365;                        print, '          [atmos grid : take ocean points only] ' 
     366;                     END 
     367;                     2: BEGIN 
     368;                        ; take land points only 
     369;                        idx = where(tmask GT 0) 
     370;                        lec_data(idx) = 1.e20 
     371;                        print, '          [atmos grid : take land points only] ' 
     372;                     END 
     373;                  ENDCASE  
     374;                  CASE vert_type OF 
     375;                     'z': BEGIN ;  depth range 
     376;                        IF debug_w THEN print, ' performing average' 
     377;                        zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
     378;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     379;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
     380;                        ENDIF ELSE BEGIN 
     381;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     382;                        ENDELSE  
     383;                     END  
     384;                     ELSE: BEGIN ; levels range 
     385;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     386;                           zmean = lec_data(*, *, vert_mean[0]) 
     387;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
     388;                        ENDIF ELSE BEGIN 
     389;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) 
     390;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     391;                        ENDELSE  
     392;                     END  
     393;                  ENDCASE  
     394;;                  tmask = reform(tmask(*, *, 0), jpi, jpj) 
     395;               ENDIF ELSE BEGIN 
     396;;                  ; ocean = always mask 
     397;;                  idx = where(tmask EQ 0) 
     398;;                  lec_data(idx) = 1.e20 
     399;                  CASE vert_type OF 
     400;                     'z': BEGIN 
     401;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     402;                           zmean = lec_data ; right depth already selected 
     403;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
     404;                        ENDIF ELSE BEGIN 
     405;                           zmean = moyenne(lec_data, 'z', boite = vert_mean) 
     406;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     407;                        ENDELSE  
     408;                     END  
     409;                     ELSE: BEGIN ; case level (zindex) 
     410;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     411;                           zmean = lec_data ; right depth already selected 
     412;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
     413;                        ENDIF ELSE BEGIN 
     414;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) 
     415;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     416;                        ENDELSE  
     417;                     END  
     418;                  ENDCASE  
     419;               ENDELSE  
     420;               lec_data = zmean 
     421;               domdef, old_boite 
     422;               field_dim = field_dim - 1 
     423;               data_direc = 'xy' 
     424;               nzt = 1 
     425;               firstz = 0 
     426;               lastz = 0 
     427;            ENDIF  
     428; 
     429;         ENDIF ELSE BEGIN  
     430;            print, '    Reading ', var_name, ' (3D data time serie)',  strcompress(string(time_1)),'-', strtrim(string(time_2), 2), '   from ', file_name 
     431;            ; read vertical levels 
     432;            IF mesh_type EQ 'atm' THEN BEGIN 
     433;               ncdf_diminq,cdfid,2, name_level, nb_level 
    364434 
    365435               ; get name of vertical level from metadata 
    366                file_grid_config = hom_def+'grid_config.def' 
    367                spawn, 'grep -i "\ '+meshlec_type+' " '+file_grid_config, line 
    368                line = strcompress(strtrim(line[0], 2))  
    369                length = strlen(line) 
    370     
    371                IF length EQ 0 THEN BEGIN 
    372                   print, '    *** nc_read : define name_level from grid ', meshlec_type, $ 
    373                    ' in file ', file_grid_config 
    374                   name_level = '' 
    375                   return, -1 
    376                ENDIF ELSE BEGIN  
    377                   argvar = str_sep(line, ' ')  
    378                   name_level = argvar[1] 
    379                ENDELSE 
    380                dimidl = ncdf_dimid(cdfid, name_level) 
    381                ncdf_diminq,cdfid,dimidl, name_level, nb_level 
    382  
    383                varidl = ncdf_varid(cdfid, name_level) 
    384                ; make sure name_level is in hPa 
    385                ncdf_attget, cdfid, varidl, 'units', val_unit 
    386                IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar' AND string(val_unit) NE 'mbar' THEN BEGIN  
    387                   print, '  vertical levels coordinate not obvious' 
    388                   print, '  ... using <levels> ...' 
    389                   varidl = ncdf_varid(cdfid, 'levels') 
    390                ENDIF  
    391                ncdf_varget, cdfid, varidl, gdept  
    392                gdepw = gdept 
    393                e3t = shift(gdept, 1)-gdept 
    394                e3t[0] = e3t[1] 
    395                e3w = e3t 
    396                jpk = nb_level 
    397                tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) 
    398             ENDIF  
    399             @read_ncdf_varget 
    400             lec_data = res 
    401             data_direc = 'xyzt' 
    402  
    403 ;           vertical average ? 
    404             IF vert_switch ge 1  AND NOT keyword_set(no_mean) THEN BEGIN  
    405                old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    406                print, '      Average in vertical domain ', vert_type, vert_mean 
    407                IF mesh_type EQ 'atm' THEN BEGIN 
    408                   CASE atmos_msk OF  
    409                      0: print, '          [take all points] ' ; take all points 
    410                      1: BEGIN 
    411                         ; take ocean points only 
    412                         idx = where(tmask EQ 0) 
    413                         lec_data(idx) = 1.e20 
    414                         print, '          [take ocean points only] ' 
    415                      END 
    416                      2: BEGIN 
    417                         ; take land points only 
    418                         idx = where(tmask GT 0) 
    419                         lec_data(idx) = 1.e20 
    420                         print, '          [take land points only] ' 
    421                      END 
    422                   ENDCASE  
    423                   CASE vert_type OF 
    424                      'z': BEGIN ; levels 
    425                          
    426                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    427                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
    428                         ENDIF ELSE BEGIN 
    429                            zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
    430                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean(0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']' 
    431                         ENDELSE  
    432                      END  
    433                      ELSE: BEGIN ; levels 
    434                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    435                            zmean = lec_data(*, *, vert_mean[0]) 
    436                            name_suff = ' at '+vert_type+' '+strtrim(string(long(vert_mean(0))), 2)+' ['+strtrim(string(gdept(vert_mean(0))), 2)+' hPa]' 
    437                         ENDIF ELSE BEGIN 
    438                            zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) 
    439                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    440                         ENDELSE  
    441                      END  
    442                   ENDCASE  
    443 ;                  tmask = reform(tmask(*, *, 0), jpi, jpj) 
    444                ENDIF ELSE BEGIN 
    445                   ; ocean = always mask 
    446                   ; idx = where(tmask EQ 0) 
    447                   ; lec_data(idx) = 1.e20 
    448                   CASE vert_type OF 
    449                      'z': BEGIN 
    450                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    451                            zmean = lec_data 
    452                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    453                         ENDIF ELSE BEGIN 
    454                            zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
    455                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    456                         ENDELSE  
    457                      END  
    458                      ELSE: BEGIN 
    459                         IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    460                            zmean = lec_data 
    461                            name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    462                         ENDIF ELSE BEGIN 
    463                            zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) 
    464                            name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    465                         ENDELSE  
    466                      END  
    467                   ENDCASE  
    468                ENDELSE  
    469                lec_data = zmean 
    470                domdef, old_boite 
    471                field_dim = field_dim - 1 
    472                data_direc = 'xyt' 
    473                nzt = 1 
    474                firstz = 0 
    475                lastz = 0 
    476             ENDIF  
    477          ENDELSE  
    478       END  
    479      ;; erreur si dim > 4 
    480       ELSE: BEGIN 
    481          err_mess =  '  *** nc_read : ERROR dimension > 4' 
    482          lec_data = -1.0 
    483       END 
    484    ENDCASE  
    485  
    486 ; scaling ? 
    487    FOR i = 0, varinq.natts-1 DO BEGIN 
    488       att_txt = ncdf_attname(cdfid, varid, i) 
    489       IF att_txt EQ 'scale_factor' THEN BEGIN 
    490          ncdf_attget, cdfid, varid, att_txt, valscale 
    491          lec_data = lec_data*valscale 
    492       ENDIF  
    493    ENDFOR  
     436;               file_grid_config = hom_def+'grid_config.def' 
     437;               spawn, 'grep -i "\ '+meshlec_type+' " '+file_grid_config, line 
     438;               line = strcompress(strtrim(line[0], 2))  
     439;               length = strlen(line) 
     440    
     441;               IF length EQ 0 THEN BEGIN 
     442;                  print, '    *** nc_read : define name_level from grid ', meshlec_type, $ 
     443;                   ' in file ', file_grid_config 
     444;                  name_level = '' 
     445;                  return, -1 
     446;               ENDIF ELSE BEGIN  
     447;                  argvar = strsplit(line, ' ', /EXTRACT)  
     448;                  name_level = argvar[1] 
     449;               ENDELSE 
     450;               dimidl = ncdf_dimid(cdfid, name_level) 
     451;               ncdf_diminq,cdfid,dimidl, name_level, nb_level 
     452; 
     453;               varidl = ncdf_varid(cdfid, name_level) 
     454;               ; make sure name_level is in hPa 
     455;               ncdf_attget, cdfid, varidl, 'units', val_unit 
     456;               IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar' AND string(val_unit) NE 'mbar' THEN BEGIN  
     457;                  print, '  vertical levels coordinate not obvious' 
     458;                  print, '  ... using <levels> ...' 
     459;                  varidl = ncdf_varid(cdfid, 'levels') 
     460;               ENDIF  
     461;               ncdf_varget, cdfid, varidl, gdept  
     462;               gdepw = gdept 
     463;               e3t = shift(gdept, 1)-gdept 
     464;               e3t[0] = e3t[1] 
     465;               e3w = e3t 
     466;               jpk = nb_level 
     467;               tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) 
     468;            ENDIF  
     469;            @read_ncdf_varget 
     470;            lec_data = res 
     471;            data_direc = 'xyzt' 
     472; 
     473;;           vertical average ? 
     474;            IF vert_switch ge 1  AND NOT keyword_set(no_mean) THEN BEGIN  
     475;               old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
     476;               print, '      Average in vertical domain ', vert_type, vert_mean 
     477;               IF mesh_type EQ 'atm' THEN BEGIN 
     478;                  CASE atmos_msk OF  
     479;                     0: print, '          [take all points] ' ; take all points 
     480;                     1: BEGIN 
     481;                        ; take ocean points only 
     482;                        idx = where(tmask EQ 0) 
     483;                        lec_data(idx) = 1.e20 
     484;                        print, '          [take ocean points only] ' 
     485;                     END 
     486;                     2: BEGIN 
     487;                        ; take land points only 
     488;                        idx = where(tmask GT 0) 
     489;                        lec_data(idx) = 1.e20 
     490;                        print, '          [take land points only] ' 
     491;                     END 
     492;                  ENDCASE  
     493;                  CASE vert_type OF 
     494;                     'z': BEGIN ; levels 
     495;                         
     496;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     497;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
     498;                        ENDIF ELSE BEGIN 
     499;                           zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
     500;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean(0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']' 
     501;                        ENDELSE  
     502;                     END  
     503;                     ELSE: BEGIN ; levels 
     504;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     505;                           zmean = lec_data(*, *, vert_mean[0]) 
     506;                           name_suff = ' at '+vert_type+' '+strtrim(string(long(vert_mean(0))), 2)+' ['+strtrim(string(gdept(vert_mean(0))), 2)+' hPa]' 
     507;                        ENDIF ELSE BEGIN 
     508;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) 
     509;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     510;                        ENDELSE  
     511;                     END  
     512;                  ENDCASE  
     513;;                  tmask = reform(tmask(*, *, 0), jpi, jpj) 
     514;               ENDIF ELSE BEGIN 
     515;                  ; ocean = always mask 
     516;                  ; idx = where(tmask EQ 0) 
     517;                  ; lec_data(idx) = 1.e20 
     518;                  CASE vert_type OF 
     519;                     'z': BEGIN 
     520;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     521;                           zmean = lec_data 
     522;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
     523;                        ENDIF ELSE BEGIN 
     524;                           zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
     525;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     526;                        ENDELSE  
     527;                     END  
     528;                     ELSE: BEGIN 
     529;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
     530;                           zmean = lec_data 
     531;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
     532;                        ENDIF ELSE BEGIN 
     533;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) 
     534;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
     535;                        ENDELSE  
     536;                     END  
     537;                  ENDCASE  
     538;               ENDELSE  
     539;               lec_data = zmean 
     540;               domdef, old_boite 
     541;               field_dim = field_dim - 1 
     542;               data_direc = 'xyt' 
     543;               nzt = 1 
     544;               firstz = 0 
     545;               lastz = 0 
     546;            ENDIF  
     547;         ENDELSE  
     548;      END  
     549;     ;; erreur si dim > 4 
     550;      ELSE: BEGIN 
     551;         err_mess =  '  *** nc_read : ERROR dimension > 4' 
     552;         lec_data = -1.0 
     553;      END 
     554;   ENDCASE  
    494555 
    495556; Attribut du champ 
    … …  
    498559   field.name = var_name 
    499560   field.origin = directory+file_name 
     561    
    500562; 
    501563; get long name  
    502564;   result = '???' 
    503    FOR i = 0, varinq.natts-1 DO BEGIN 
    504       att_txt = ncdf_attname(cdfid, varid, i) 
    505       IF att_txt EQ 'long_name' OR att_txt EQ 'title' THEN ncdf_attget, cdfid, varid, att_txt, result 
    506    ENDFOR  
    507    field.legend = strtrim(string(result),1)+name_suff 
     565;   FOR i = 0, varinq.natts-1 DO BEGIN 
     566;      att_txt = ncdf_attname(cdfid, varid, i) 
     567;      IF att_txt EQ 'long_name' OR att_txt EQ 'title' THEN ncdf_attget, cdfid, varid, att_txt, result 
     568;   ENDFOR  
     569    
     570;   field.legend = strtrim(string(result),1)+name_suff 
     571 
     572;;;;   ncdf_getatt, directory+file_name, var_name, add_offset = add_offset, scale_factor = scale_factor, missing_value = missing_value, units = units, calendar = calendar, long_name = long_name 
     573   field.legend = long_name+name_suff 
     574 
     575; scaling ? 
     576;   FOR i = 0, varinq.natts-1 DO BEGIN 
     577;      att_txt = ncdf_attname(cdfid, varid, i) 
     578;      IF att_txt EQ 'scale_factor' THEN BEGIN 
     579;         ncdf_attget, cdfid, varid, att_txt, valscale 
     580;         lec_data = lec_data*valscale 
     581;      ENDIF  
     582;   ENDFOR  
     583 
    508584 
    509585; get unit 
    510586; if it exists 
    511   isunits = ncdf_attinq(cdfid,  varid,  'units') 
    512   IF isunits.datatype NE 'UNKNOWN' THEN BEGIN 
    513      ncdf_attget, cdfid, varid, 'units', result 
    514      field.units = strtrim(string(result),1) 
    515   ENDIF ELSE BEGIN 
    516      print,  'No units for the variable ',  field.name 
    517      print,  ' Or units were forgotten in the file !' 
    518   ENDELSE 
    519   field.dim = field_dim  
    520  
     587;  isunits = ncdf_attinq(cdfid,  varid,  'units') 
     588;  IF isunits.datatype NE 'UNKNOWN' THEN BEGIN 
     589;     ncdf_attget, cdfid, varid, 'units', result 
     590;     field.units = strtrim(string(result),1) 
     591;  ENDIF ELSE BEGIN 
     592;     print,  'No units for the variable ',  field.name 
     593;     print,  ' Or units were forgotten in the file !' 
     594;  ENDELSE 
     595   field.units = units 
     596   field.dim = field_dim  
    521597 
    522598;  get valmask (might need valmask = float(string(valmask)) 
    523599 
    524    valmask = 1.e20 
    525    FOR i = 0, varinq.natts-1 DO BEGIN 
    526       att_txt = ncdf_attname(cdfid, varid, i) 
    527       IF att_txt EQ 'missing_value' OR att_txt EQ 'mask value' OR att_txt EQ '_FillValue' THEN ncdf_attget, cdfid, varid, att_txt, valmask 
    528    ENDFOR  
    529  
     600;   valmask = 1.e20 
     601;   FOR i = 0, varinq.natts-1 DO BEGIN 
     602;      att_txt = ncdf_attname(cdfid, varid, i) 
     603;      IF att_txt EQ 'missing_value' OR att_txt EQ 'mask value' OR att_txt EQ '_FillValue' THEN ncdf_attget, cdfid, varid, att_txt, valmask 
     604;   ENDFOR  
     605   valmask = missing_value 
    530606; set valmask to 1.e20 
    531607 
    … …  
    558634   print, '      = ',field.legend, '    [min/max = ',minf , maxf,'  ', field.units,' - masked values = ',valmask, chardim, ']' 
    559635 
    560  
    561    ncdf_close, cdfid 
    562  
    563 ;key_offset = key_offset_orig 
    564 ;jpi = jpi_orig 
    565 ;jpj = jpj_orig 
    566 ;jpk = jpk_orig 
     636;   ncdf_close, cdfid 
    567637 
    568638   IF keyword_set(key_yreverse) THEN print, '    key_yreverse active : ',  key_yreverse 
    … …  
    570640   IF debug_w THEN print,  '  ...EXIT nc_read' 
    571641   IF debug_w THEN print,  '  ' 
    572  
     642    
    573643   return, field 
    574644 
    575  
    576645END  

  • branches/procs-read/trends.pro

    r73 r81  
    360360         print,  '     -> writing 1D data to ',  asciidir+filename 
    361361         print,  ' ' 
    362          printf, nuldat, fld*scale, format = '(f8.3)' 
     362         printf, nuldat, fld*scale, format = '(f8.4)' 
    363363         free_lun, nuldat 
    364364         close, nuldat 
Note: See TracChangeset for help on using the changeset viewer.