Changeset 86


Ignore:
Timestamp:
04/07/08 15:50:01 (16 years ago)
Author:
kolasinski
Message:

Merged procs-read branch changes r80:85 into teh trunk

Location:
trunk
Files:
22 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/config/domain_boxes.def

    r78 r86  
    167167  160_140sW                200/220/-10/10/20/28                  1 
    168168  west_pac                120/180/-2/2                  1 
     169  west_pac_atm                120/180/-2/2/0/100000                  1 
    169170  west2_pac                140/180/-30/30                  1 
    170171  westw_pac                120/165/-2/2                  1 

  • trunk/config/fld_glo_mmx.def

    r78 r86  
    175175 #  vo fields 
    176176 # 
    177  votemper     -2         30       24.5       28      -2      30      1 
     177 temp    @votemper  
     178 votemper     -2         30       15       28      -2      30      1 
    178179 vosaline     30       40       31       34         33     34   2  
    179180 vodenpot     20       30    20       30    20       30      4 

  • trunk/procs/com_eg.pro

    r41 r86  
    44; attributs horizontal+vertical means 
    55; 
    6 COMMON zm_att, box_h, depth_z, zoom_z, diaznl_idx, box_plot, legbox, hpa_max, hpa_min, max_spec, spec_win, vert_type, vert_mean,  vert_switch, glamboundary_box, msf_mean, name_level 
     6COMMON zm_att, box_h, depth_z, zoom_z, diaznl_idx, box_plot, legbox, pres_max, pres_min, max_spec, spec_win, vert_type, vert_mean,  vert_switch, glamboundary_box, msf_mean, name_level 
    77; 
    88; attributs hoevmoeller 

  • trunk/procs/data_read.pro

    r60 r86  
    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 
     93; define horizontal domain and vertical domain if needed (used in 
     94; read_ncdf and update_data) 
     95 
     96   box_plot = def_box(cmd.plt, dimplot, legbox, time_stride) 
     97   CASE n_elements(box_plot) OF 
     98      4 : IF vert_switch GE 1 THEN box_plot = [box_plot, vert_mean] 
     99      6 : IF vert_switch GE 1 THEN box_plot = [box_plot[0:3],  vert_mean] ELSE box_plot = box_plot[0:3] 
     100   ENDCASE 
     101; Exceptions : need to read all the vertical data (and not a subset) 
     102; for the pltz case. Other cases ?    
     103   IF plttyp EQ 'pltz' THEN box_plot = box_plot[0:3] 
     104 
     105; define timetsteps time1 and time2 
     106 
     107   IF hotyp NE '-' THEN BEGIN 
     108      time1 = delta_t1+1 
    106109      timearr = compute_time(cmd.timave, cmd.date1, cmd.spec) 
    107       jpt = timearr.count 
    108       time2 = jpt+delta_t1 
     110      time2 =  timearr.count+delta_t1 
    109111      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  
     112   ENDIF ELSE BEGIN 
     113      IF strpos(cmd.timave, 'mm') NE -1 THEN time1 = delta_t1+long(strmid(cmd.date1, 0, 2)) $ 
     114      ELSE  time1 = delta_t1+1 
    126115      time2 = time1 
    127116      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 
     117   ENDELSE 
    157118 
    158119; save time in common fld_att 
    … …  
    166127      '@@': BEGIN 
    167128         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) 
     129         fldr = macro_read(file_name, cmd.var, ncdf_db, BOXZOOM = box_plot, TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
    169130         IF stddev_diff EQ 1 THEN fldr.origin =  'diff' 
    170131      END  
    … …  
    174135               idx = strpos(cmd.var, '=f(') 
    175136               var1 = strmid(cmd.var, 0, idx) 
    176                fldr1 = nc_read(file_name, var1, ncdf_db, $ 
     137               fldr1 = nc_read(file_name, var1, ncdf_db, BOXZOOM = box_plot, $ 
    177138                               TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
    178139               file_name1 = file_name 
    … …  
    185146               IF strpos(cmd.var, '=f(next)') EQ -1 THEN BEGIN  ; type 1 y=f(x) on 1 line 
    186147                  var2 = strmid(cmd.var, idx+3, strlen(cmd.var)-idx-4) 
    187                   fldr2 = nc_read(file_name, var2, ncdf_db, $ 
     148                  fldr2 = nc_read(file_name, var2, ncdf_db, BOXZOOM = box_plot, $ 
    188149                                  TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
    189150                  datyp2 = datyp 
    … …  
    294255 
    295256               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  
     257                  fldr = nc_read(file_name, cmd.var, ncdf_db, BOXZOOM = box_plot, $ 
     258                                  TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex) 
     259                  ;fldr = nc_read(file_name, cmd.var, ncdf_db, $ 
     260                  ;               TIME_1 = time1,  TIME_2 =  time2, ALL_DATA = all_data, _extra = ex)               
    299261                  ; perform mean if hovmoeller diagram is required 
    300262                  ; It is absolutely necessary for the division case with  
    … …  
    316278   ENDCASE   
    317279 
    318 ; if tkeavt, take log 
    319  
    320 ;    CASE cmd.var OF  
    321 ;       'votkeavt': BEGIN 
    322 ;          idx = where(fldr.data LE valmask/10) 
    323 ;          idxm = where(fldr.data EQ valmask) 
    324 ;          fldr.data(idx) = alog10(fldr.data(idx)) 
    325 ;          ; on T grid 
    326 ;          vargrid = 'T' 
    327 ;          fldr.data = (fldr.data+shift(fldr.data, 0, 0, 1))*.5 
    328 ;          fldr.data(idxm) = valmask 
    329 ;       END  
    330 ;       ELSE:  
    331 ;    ENDCASE   
     280; Redefinition of time because it is updated in read_ncdf 
     281; Useful for pltt routine 
     282 
     283   IF hotyp NE '-' THEN BEGIN 
     284      time = timearr.scale 
     285   ENDIF ELSE BEGIN 
     286      time = time1 
     287   ENDELSE 
    332288 
    333289 
    334290; density projection 
     291 
    335292   IF splot EQ 1 THEN BEGIN 
    336293      IF cmd.timave EQ '1y' AND 1 THEN BEGIN  
    … …  
    354311            print,  '   Date1 = ',  cmd.date1 
    355312            sfild = fldr 
    356             bin_sigma, cmd, sfild    ;;;;;;;;;;;;;; prob 
     313            bin_sigma, cmd, sfild, BOXZOOM = box_plot, ALL_DATA = all_data ;;;;;;;;;;;;;; prob 
    357314            cumulative = cumulative + sfild.data 
    358315         ENDFOR 
    … …  
    369326               print, 'data_read: Performing monthly bining. Indices: first, last, current = 0,',  timedim-1,  i 
    370327               sfild2 = {name: sfild.name, data: sfild.data(*, *, *, i), legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim-1} 
    371                bin_sigma, cmd, sfild2 
     328               bin_sigma, cmd, sfild2, BOXZOOM = box_plot, ALL_DATA = all_data  
    372329               fldrd(*, *, *, i)= sfild2.data 
    373330            ENDFOR  
    … …  
    375332         ENDIF ELSE BEGIN 
    376333            sfild = fldr 
    377             bin_sigma, cmd, sfild 
     334            bin_sigma, cmd, sfild, BOXZOOM = box_plot, ALL_DATA = all_data 
     335            ;;bin_sigma, cmd, sfild, ALL_DATA = all_data 
    378336            fldr = sfild 
    379337         ENDELSE  
    … …  
    413371 
    414372      ; perform difference 
    415  
    416373      diff = fldr.data-field2.data 
    417374      IF (where(fldr.data EQ valmask))[0] NE -1 THEN $ 
    … …  
    430387       
    431388      ; read field2 
    432  
    433389      field2 = data_read(cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex) 
    434390      IF n_elements(field2.data) EQ 1 THEN return, field2 
    … …  
    471427          
    472428         ; decode field2 
    473  
    474429         argvar = strsplit(cmdl, '/', /EXTRACT) 
    475430 
    … …  
    482437          
    483438         ; read field2 
    484  
    485439         field2 = data_read(cmd2, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex) 
    486440         IF n_elements(field2.data) EQ 1 THEN return, field2 
    487441 
    488442         ; perform difference 
    489  
    490443         diff = fldr.data-field2.data 
    491444         IF (where(fldr.data EQ valmask))[0] NE -1 THEN $ 
    … …  
    508461ENDELSE 
    509462 
     463; Pb with varexp which is always updated in read_ncdf 
     464; For Seb, varexp is the name of the file. For Eric, varexp is the name of the experiment 
     465varexp = cmd.exp 
     466 
    510467IF debug_w THEN print, '       varexp before exit in data_read = ', varexp 
    511468IF debug_w THEN print, '   ...EXIT data_read' 

  • trunk/procs/macro_read.pro

    r48 r86  
    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:''} 

  • trunk/procs/macros/make_depth.pro

    r68 r86  
    11;------------------------------------------------------------ 
    22;------------------------------------------------------------ 
    3 FUNCTION make_depth, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ZMTYP = zmtyp 
     3FUNCTION make_depth, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ZMTYP = zmtyp, _EXTRA = ex 
    44; 
    55@common 

  • trunk/procs/macros/make_eos.pro

    r2 r86  
    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, ALL_DATA = all_data, no_mean = 1) 
     60   sn = nc_read(file_name,'vosaline', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, no_mean = 1) 
    6161 
    6262; declarations 

  • trunk/procs/macros/make_msf.pro

    r2 r86  
    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, ALL_DATA = all_data, _extra = ex) 
    1717 
    1818   idx = where(v.data EQ valmask) 

  • trunk/procs/macros/make_stddev.pro

    r17 r86  
    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 

  • trunk/procs/mask_z.pro

    r2 r86  
    169169               CASE strmid(cmd.plt, 1, 1) OF  
    170170                  't': boite_pltz = box_h 
    171                   ELSE : boite_pltz = [box_h, hpa_min, hpa_max ] 
     171                  ELSE :  boite_pltz = [box_h, pres_min, pres_max ] 
    172172               ENDCASE  
    173173            ENDIF  
    … …  
    191191               CASE strmid(cmd.plt, 1, 1) OF  
    192192                  't': boite_pltz = zbox 
    193                   ELSE : IF n_elements(boite_pltz) EQ 4 THEN boite_pltz = [zbox, hpa_min, hpa_max] 
     193                  ELSE : IF n_elements(boite_pltz) EQ 4 THEN boite_pltz = [zbox, pres_min, pres_max] 
    194194               ENDCASE  
    195195            ENDIF  

  • trunk/procs/meshes/mesh_from_file.pro

    r69 r86  
    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 

  • trunk/procs/meshes/mesh_orca.pro

    r26 r86  
    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 

  • trunk/procs/msf_simple.pro

    r2 r86  
    4545;  integration zonale du flux ! 
    4646; 
    47 fm = total(z, 1)  
     47fm = total(z, 1, /NAN)  
    4848; 
    4949;  calcul de la msf en integrant depuis le fond 

  • trunk/procs/nc_read.pro

    r41 r86  
    11;---------------------------   
    2 ; Reading ORCA netcdf files 
    3 ; EG 24/2/99 
     2; Reading any netcdf file 
    43;---------------------------   
    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 
     4FUNCTION 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 
    65 
    76; arguments = file_name, varname 
    … …  
    1110@common 
    1211@com_eg 
    13  
    14 ;;   print,  'key_yreverse : ', key_yreverse   
    15  
     12   
    1613; inits 
    17    IF debug_w THEN print, ' ' 
    18    IF debug_w THEN print, '  ENTER nc_read' 
     14   IF debug_w THEN print, '   ' 
     15   IF debug_w THEN print, '   ENTER nc_read...' 
    1916 
    2017   IF NOT keyword_set(TIME_1) THEN time_1 =  1 
    2118   IF NOT keyword_set(TIME_2) THEN time_2 =  time_1 
    22 ; 
    23 ; decide which subdomain of data to read 
    24 ; 
    25  
    26    IF debug_w THEN print, '     keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA)  
    27  
    28    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 
    39 ;Rajout MK 
    40 ;Trouver une meilleure place 
    41       ixminmesh = 0 
    42       iyminmesh = 0 
    43       ixmaxmesh = jpi-1 
    44       iymaxmesh = jpj-1 
    45 ;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 
    5619 
    5720 
    58 ; 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  
    65  
    66 ; force izmaxmesh to lastz 
    67    jpk = lastz+1 
    68  
    69 ; 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: 
     21; decide which subdomain of data to read 
     22   IF debug_w THEN print, '     keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA) 
     23   IF keyword_set(ALL_DATA) THEN BEGIN 
     24      tout = 1 
     25   ENDIF ELSE tout = 0 
     26   CASE vert_type OF 
     27      'z' : zindex = 0 
     28      'level' : zindex = 1 
     29      ELSE: zindex = 0 
     30   ENDCASE 
    9231 
    9332; local directory 
    94  
    95    IF strpos(ncdf_db, ':') GE 1 THEN directory = (str_sep(ncdf_db, ':'))[1] $ 
     33   IF strpos(ncdf_db, ':') GE 1 THEN directory = (strsplit(ncdf_db, ':', /EXTRACT))[1] $ 
    9634    ELSE directory = ncdf_db 
    9735 
    98 ; existence of file 
     36; Find the variable's attribute 
     37   ncdf_getatt, directory+file_name, var_name, add_offset = add_offset, scale_factor = scale_factor, $ 
     38    missing_value = missing_value, units = units, calendar = calendar, long_name = long_name 
    9939 
    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  
     40; By default units for the Z axis are meters 
     41   IF n_elements(gdept) GT 1 THEN BEGIN 
     42      IF name_level NE '-' THEN $ 
     43       ncdf_getatt, directory+file_name, name_level, units = units_depth ELSE units_depth = 'm' 
     44   ENDIF ELSE units_depth = '' 
    11045 
    111 ; ouverture fichier netCDF + contenu 
    112    cdfid=ncdf_open(directory+file_name)  
    113    inq=ncdf_inquire(cdfid) 
    114 ; que contient la variable 
    115    varid = ncdf_varid(cdfid, var_name) 
     46; Check consistency between time_2 computed from param in cmdline and the max number 
     47; of time steps found in the file. 
     48   jpt_max = find_jptmax(file_name, IODIRECTORY = directory) 
     49   IF time_2-time_1+1 GT jpt_max AND jpt_max NE -1 THEN time_2 = time_1+jpt_max-1 
    11650 
    117    name_suff = '' 
     51; Read the data with a call to read_ncdf 
     52   IF debug_w THEN print, '     ' 
     53   IF debug_w THEN print, '     Reading raw data from ',  file_name 
    11854 
    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) 
    128  
    129 ; test sur la dimension 
    130    err_mess = '' 
    131    field_dim = n_elements(varinq.dim) 
    132  
    133 ; 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 
     55   lec_data = read_ncdf(var_name, time_1-1, time_2-1, BOXZOOM = boxzoom, FILENAME = directory+file_name, $ 
     56                        /TIMESTEP, /ADDSCL_BEFORE, TOUT = tout, /NOSTRUCT, /CONT_NOFILL, GRID = vargrid, $ 
     57                        ZINVAR = zinvar, ZINDEX = zindex) 
     58   field_dim = (size(lec_data))[0] 
     59   IF field_dim LE 0 THEN stop,  '  Something wrong happened in read_ncdf ' 
    19460 
    19561   IF debug_w THEN print, '     ' 
    196    IF debug_w THEN print, '     key=', key  
    197    IF debug_w THEN print, '     jpitotal=', jpitotal 
     62   IF debug_w THEN print,  '    boxzoom=', boxzoom 
     63   IF debug_w THEN print,  '    firstxt,firstyt,firstzt=', firstxt, firstyt, firstzt 
     64   IF debug_w THEN print,  '    lastxt,lastyt,lastzt=', firstxt, firstyt, firstzt    
     65   IF debug_w THEN print,  '    zinvar=', zinvar 
    19866   IF debug_w THEN print, '     ixmindta,iymindta,izmindta  =', ixmindta,iymindta,izmindta 
    19967   IF debug_w THEN print, '     ixminmesh,iyminmesh=', ixminmesh,iyminmesh 
    20068   IF debug_w THEN print, '     ixmaxmesh,iymaxmesh=', ixmaxmesh,iymaxmesh  
    20169   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 
    204    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 
     70   IF debug_w THEN print, '     jpt=', jpt 
     71   IF debug_w THEN print, '     key_shift=', key_shift 
     72   IF debug_w THEN print, '     key_yreverse=',key_yreverse  
    20673   IF debug_w THEN print, '     ' 
     74    
     75; Average data along vertical if needed and update some features 
     76; needed for plt (data_direc, name_suff) 
     77   name_suff = '' 
     78   data_direc = '' 
     79   update_data, TAB = lec_data, VNAME = var_name, BOXZOOM = boxzoom, ZUNITS = units_depth, $ 
     80    ZINVAR = zinvar, SUFFIX = name_suff, D_DIREC = data_direc, $ 
     81    TIME_1 = time_1, TIME_2 =  time_2, NO_MEAN = no_mean 
     82   field_dim = (size(lec_data))[0] 
    20783 
    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  
    267   
    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 
    280 ;            lec_data =  reform(lec_data,x_count, y_count, z_count)  
    281             data_direc = 'xyz' 
    282             IF debug_w THEN help, lec_data 
    283  
    284 ;           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 
    364  
    365                ; 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  
    494  
    495 ; Attribut du champ 
    496  
     84; Field attributes 
    49785   field = {name: '', data: lec_data, legend: '', units: '', origin: '', dim: 0, direc: data_direc} 
    49886   field.name = var_name 
    49987   field.origin = directory+file_name 
    500 ; 
    501 ; get long name  
    502 ;   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 
     88   field.legend = long_name+name_suff 
     89   field.units = units 
     90   field.dim = field_dim  
    50891 
    509 ; get unit 
    510 ; 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  
     92; get valmask (might need valmask = float(string(valmask)) 
    52093 
    521  
    522 ;  get valmask (might need valmask = float(string(valmask)) 
    523  
    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  
    530 ; set valmask to 1.e20 
    531  
     94; valmask = 1.e20 
     95   IF size(missing_value,  /TYPE) EQ 4 OR size(missing_value,  /TYPE) EQ 5 THEN BEGIN 
     96      valmask = missing_value 
    53297; ensure valmask is positive 
    533  
    534    IF valmask LT 0 THEN BEGIN 
    535       print, '      *** Warning valmask is negative - changing sign: ', valmask 
    536       idx_t = where (field.data EQ valmask) 
    537       IF idx_t(0) NE -1 THEN field.data(idx_t) = -valmask 
    538       valmask = -valmask 
    539       idx_t=0 ; free memory 
     98      IF valmask LT 0 THEN BEGIN 
     99         print, '      *** Warning valmask is negative - changing sign: ', valmask 
     100         idx_t = where (field.data EQ valmask) 
     101         IF idx_t(0) NE -1 THEN field.data(idx_t) = -valmask 
     102         valmask = -valmask 
     103         idx_t=0                ; free memory 
     104      ENDIF 
    540105   ENDIF 
    541106 
    … …  
    558123   print, '      = ',field.legend, '    [min/max = ',minf , maxf,'  ', field.units,' - masked values = ',valmask, chardim, ']' 
    559124 
    560  
    561    ncdf_close, cdfid 
    562  
    563 ;key_offset = key_offset_orig 
    564 ;jpi = jpi_orig 
    565 ;jpj = jpj_orig 
    566 ;jpk = jpk_orig 
    567  
    568    IF keyword_set(key_yreverse) THEN print, '    key_yreverse active : ',  key_yreverse 
    569  
    570125   IF debug_w THEN print,  '  ...EXIT nc_read' 
    571126   IF debug_w THEN print,  '  ' 
    572  
     127    
    573128   return, field 
    574129 
    575  
    576130END  

  • trunk/procs/plt_map.pro

    r70 r86  
    598598         ENDELSE  
    599599         END  
    600       'atm': BEGIN & type_yz = ',type_yz=''hPa''' & zoom_txt = ', zoom='+string(hpa_max) & END  
     600      'atm': BEGIN & type_yz = ',type_yz=''Pa''' & zoom_txt = ', zoom='+string(pres_max) & END  
    601601      ELSE:type_yz = '' 
    602602   ENDCASE  

  • trunk/procs/trends.pro

    r73 r86  
    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 
    … …  
    375375      ENDELSE 
    376376   ENDIF  
    377  
    378    return, fld 
    379377    
    380378   IF debug_w THEN print, '   Vargrid : ',  vargrid 
    381379   IF debug_w THEN print, '   Leaving trends.pro' 
    382380 
     381   return, fld 
     382 
    383383END   
    384384       

  • trunk/procs/yfx.pro

    r71 r86  
    3535               IF debug_w THEN print, '   cmd1_back = ', cmd1_back 
    3636               IF debug_w THEN print, '   cmd2_back = ', cmd2_back 
    37                IF debug_w THEN print, '   grid1_full = ', grid1_full 
    38                IF debug_w THEN print, '   grid2_full = ', grid2_full 
     37               IF debug_w AND sw_diffg THEN print, '   grid1_full = ', grid1_full 
     38               IF debug_w AND sw_diffg THEN print, '   grid2_full = ', grid2_full 
    3939               IF sw_diffg EQ 1 THEN BEGIN 
    4040                  jptb = jpt 

  • trunk/tools/density_binning/density_bin_IDL_gm/bin_sigma.pro

    r2 r86  
    1 PRO bin_sigma, cmd, pfildi 
     1PRO bin_sigma, cmd, pfildi, BOXZOOM = boxzoom, ALL_DATA = all_data 
    22 
    33@common 
    … …  
    2828;;------------------------------------- 
    2929; 
    30    grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz 
    31  
     30    
    3231; a- read density and bowl [if sig_bowl=1] 
    3332 
    … …  
    3534   file_nam = strmid(file_name, 0, strlen(file_name)-4)+'T.nc' 
    3635   vargrid = 'T' 
    37    sig = make_eos(file_nam, '', time_1 = time1_r, time_2 = time2_r) 
     36   sig = make_eos(file_nam, '', time_1 = time1_r, time_2 = time2_r, BOXZOOM = boxzoom, ALL_DATA = all_data) 
    3837   IF sig_bowl EQ 1 THEN BEGIN  
    3938      sobwlmax = make_sobwlmax(file_nam, '', time_1 = time1_r, time_2 = time2_r) 
    … …  
    4140      sobwlmax =  {data:0} 
    4241   ENDELSE  
     42 
     43   mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 
    4344 
    4445; b- define density grid 

  • trunk/usr/.hist_post_it

    r78 r86  
    44 ------------- 
    55  
    6 field = data_read(cmd,'-','pltz', 2, 1, ZMTYP = 'yz_merid_180W') 
     6field = data_read(cmd,'t','pltt', 1, 1, ZMTYP = 't_nino_3') 
    77  
    8  boite_pltz =       178.000      182.000     -30.0000      30.0000      0.00000 
    9       400.000 
    10 pltz,{a:fld, d:'March 1949', n:'Meridional_Velocity (m/s) merid_180W', u:'m/s', g:'V'}, -999999., 0.500000,int = 0.01000000,/yz, boite=boite_pltz, zoom= 400,petit = [ 1, 1, 1],/rempli, LANDSCAPE = 1, NOCOLORBAR = 0, NOERASE = 0,xthick=2,ythick=2,zthick=2,sepdate=' ',label=3,c_annotation=c_anot_str, format = '(f6.2)',discret=coul_gmt[idx_cb1],div=n_elements(idx_cb0)-1,cb_label=levels_gmt[idx_cb0],ndayspm=calendar_type,cont_thick=2,type_yz='m',marge=[0,0,2,2], xcharsize= 1.00000, ycharsize= 1.00000,cell_fill=2 
     8boite_pltspec=      0.00000      0.00000     0.169815      51.6981 
     9plttg,power,period, 4.19495e-06, 25.1951,boite=boite_pltspec,reverse_y=1,nocontour=1,petit = [ 1, 1, 1],/rempli, LANDSCAPE = 1, NOCOLORBAR = 1, NOERASE = 0,xthick=2,ythick=2,zthick=2,sepdate=' ',label=3,c_annotation=c_anot_str, format = '(f5.0)',discret=coul_gmt[idx_cb1],div=n_elements(idx_cb0)-1,cb_label=levels_gmt[idx_cb0],ndayspm=calendar_type,cont_thick=2,type_yz='Pa',marge=[0,0,2,2], xcharsize= 1.00000, ycharsize= 1.00000,cell_fill=2 

  • trunk/usr/init.pro

    r78 r86  
    1515; path definition 
    1616; 
    17 !path = expand_path('+' + '/.autofs/home/mklod/SVN/POST_IT') $ ; procs local home 
     17!path = expand_path('+' + '/.autofs/home/mklod/SVN/TRUNK') $ ; procs local home 
    1818      + ':' + expand_path('+' + '/.autofs/home/mklod/SAXO_new_svn') $ ; saxo home 
    1919      + ':' + expand_path('+' + !dir) 
    … …  
    2929; define default directories 
    3030; 
    31 homedir = isadirectory('/.autofs/home/mklod/SVN/POST_IT/', title = 'Select the default HOME directory') 
    32 iodir = isadirectory('/.autofs/home/mklod/SVN/POST_IT/', title = 'Select the default IO directory') 
    33 psdir = isadirectory('/.autofs/home/mklod/SVN/POST_IT/out/post_out/', title = 'Select the default postscripts directory') 
    34 imagedir = isadirectory('/.autofs/home/mklod/SVN/POST_IT/out/post_out/', title = 'Select the default images directory') 
    35 animdir = isadirectory('/.autofs/home/mklod/SVN/POST_IT/out/post_out/', title = 'Select the default animations directory') 
     31homedir = isadirectory('/.autofs/home/mklod/SVN/TRUNK/', title = 'Select the default HOME directory') 
     32iodir = isadirectory('/.autofs/home/mklod/SVN/TRUNK/', title = 'Select the default IO directory') 
     33psdir = isadirectory('/.autofs/home/mklod/SVN/TRUNK/out/post_out/', title = 'Select the default postscripts directory') 
     34imagedir = isadirectory('/.autofs/home/mklod/SVN/TRUNK/out/post_out/', title = 'Select the default images directory') 
     35animdir = isadirectory('/.autofs/home/mklod/SVN/TRUNK/out/post_out/', title = 'Select the default animations directory') 
    3636; 
    3737; define printer parameters 

  • trunk/usr/plt_def.pro

    r78 r86  
    2828; vertical domain 
    2929 
    30    depth_z = 400 
    31    zoom_z = 400 
    32  
    33    hpa_min = 10 
    34    hpa_max = 500 
     30   depth_z = 400 ; m 
     31   zoom_z = 400 ; m 
     32 
     33   pres_min = 1000  ; Pa 
     34   pres_max = 100000  ; Pa 
    3535    
    3636   msf_mean =  0 
     37 
    3738; vertical average for 3D fields  
    3839 
    39    vert_type = 'z'          ;  'z' for depth/altitude or 'level' or '0' for nothing 
    40    vert_mean = [0, 10]   ; [depth1,depth2] or [level1,level2] in C notation 0-jpk-1  
    41    ;vert_mean = [80000, 100000] 
     40   vert_type = '0'          ;  'z' for depth/altitude or 'level' or '0' for nothing 
     41   vert_mean = [0, 200]   ; [depth1,depth2] or [level1,level2] in C notation 0-jpk-1  
     42;   vert_mean = [90000, 100000] 
    4243 
    4344; density domain (sigma) + delta sigma 
    4445 
    4546   sig_min = 20. 
    46    sig_max = 27. 
     47   sig_max = 28. 
    4748   sig_del = 0.2 
    4849 
    … …  
    231232   data_domain = 'zonal' 
    232233   data_domain = 'pacific' 
    233 ;   data_domain = 'global' 
     234   data_domain = 'global' 
    234235 
    235236; grids list (IPCC atmos regular) 
    … …  
    237238 
    238239; machine type ('x' or 'WIN') 
    239         dev_type='x' 
     240   dev_type='x' 
    240241 
    241242; debug mode 
    242      debug_w = 1 
     243   debug_w = 1 
    243244 
    244245END 

  • trunk/usr/post_it.pro

    r78 r86  
    2020;;'ERA40 = local:/home2/mkdlod/database/CORREL/',  $ 
    2121'CDH4 = local:/home2/mkdlod/database/CDH4/', $ 
    22 '2L24DEV = local:/home2/mkdlod/database/', $ 
     22'2L24 = local:/home2/mkdlod/database/', $ 
    2323'NCEPDEV = local:/home2/mkdlod/database/', $ 
    2424'ERA40TDEV = local:/home2/mkdlod/database/', $  
    … …  
    3131;'ncdf_db = local:/.autofs/home/mklod/IDL_WORK/', $ 
    3232;;'ncdf_db = local:/home/ericg/backups/lemnos/home/database/',  $ 
    33 ;'ncdf_db = local:/home2/mkdlod/database/TESTS/', $ 
    34 'ncdf_db = local:/home2/mkdlod/database/TiKE/', $ 
     33'ncdf_db = local:/home2/mkdlod/database/TESTS/', $ 
     34;'ncdf_db = local:/home2/mkdlod/database/TiKE/', $ 
    3535;'ncdf_db = local:/home2/mkdlod/database/', $ 
    3636; Input databases 
    … …  
    6666other_file = 'correl' 
    6767other_file = 'post_it2' 
    68 ;other_file = 'post_it_tests' 
     68other_file = 'post_it_tests' 
    6969;other_file = '-' 
    7070 
Note: See TracChangeset for help on using the changeset viewer.