Changeset 84


Ignore:
Timestamp:
04/04/08 16:47:32 (16 years ago)
Author:
kolasinski
Message:

Clean up data_read, nc_read and update_data procedures

Location:
branches/procs-read
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/procs-read/data_read.pro

    r83 r84  
    9393; define horizontal domain and vertical domain if needed (used in 
    9494; read_ncdf and update_data) 
     95 
    9596   box_plot = def_box(cmd.plt, dimplot, legbox, time_stride) 
    9697   CASE n_elements(box_plot) OF 
    … …  
    103104 
    104105; define timetsteps time1 and time2 
     106 
    105107   IF hotyp NE '-' THEN BEGIN 
    106108      time1 = delta_t1+1 
    … …  
    114116      time = time1 
    115117   ENDELSE 
    116  
    117 ;    IF hotyp NE '-' THEN BEGIN   ; time serie 
    118 ;       ; get hovmoeller box + stride 
    119 ;       box_plot = def_box(cmd.plt, dimplot, legbox, time_stride) 
    120 ;       IF debug_w THEN print, '    box_plot for domdef in data_read = ', box_plot 
    121 ;       CASE n_elements(box_plot) OF 
    122 ;          4: domdef, box_plot[0:3], /MEMEINDICES 
    123 ;          6: domdef, box_plot[0:5], /MEMEINDICES 
    124 ;       ENDCASE  
    125 ;       ; time grid for hovmoeller 
    126 ;       time1 = 1+delta_t1 
    127 ;       IF cmd.out EQ 'cdf' THEN nb_cycles = 1 
    128 ;       timearr = compute_time(cmd.timave, cmd.date1, cmd.spec) 
    129 ;       jpt = timearr.count 
    130 ;       time2 = jpt+delta_t1 
    131 ;       time = timearr.scale 
    132 ;       niveau = 1 
    133 ; ;      print, '  Data_read, hotyp/jpt = ', hotyp, jpt 
    134 ; ;      print, '  Data_read, box = ', box_plot[0:3] 
    135 ;    ENDIF ELSE BEGIN             ; single date 
    136 ;       jpt = 1 
    137 ;       CASE strmid(cmd.timave, strlen(cmd.timave)-1, 1) OF 
    138 ;          'y': time1 = delta_t1+1 
    139 ;         'm': BEGIN 
    140 ;             CASE strmid(cmd.timave, strlen(cmd.timave)-2, 1) OF 
    141 ;                'm': time1 = delta_t1+long(strmid(cmd.date1, 0, 2))  ; mean month 
    142 ;                ELSE: time1 = delta_t1+1 
    143 ;             ENDCASE  
    144 ;             END  
    145 ;          'd': time1 = delta_t1+1 
    146 ;          ELSE: time1 = delta_t1+1 
    147 ;       ENDCASE  
    148 ;       time2 = time1 
    149 ;       time = time1 
    150 ;       print, '    Setting domdef in data_read, plt typ =', plttyp    
    151 ;       CASE plttyp OF 
    152 ;          'plt': BEGIN 
    153 ;             box_plot = def_box(cmd.plt, 2, legbox, time_stride) 
    154 ;             IF strmid(cmd.var, 0,2) EQ 'vo' THEN BEGIN ; get right zindex for vertical 
    155 ;                CASE vert_type OF 
    156 ;                   'level':domdef,[ box_plot[0:3],vert_mean],/zindex, /MEMEINDICES 
    157 ;                   'z':domdef,[ box_plot[0:3],vert_mean], /MEMEINDICES 
    158 ;                   ELSE: domdef, box_plot[0:3], /MEMEINDICES 
    159 ;                ENDCASE  
    160 ;             ENDIF ELSE BEGIN 
    161 ;                domdef, box_plot[0:3], /MEMEINDICES 
    162 ;             ENDELSE 
    163 ;          END  
    164 ;          'pltz': BEGIN 
    165 ;             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
    166 ;             domdef, box_plot[0:3], /MEMEINDICES 
    167 ;          END  
    168 ;          'plt1d': BEGIN 
    169 ;             box_plot = def_box(cmd.plt, 2, legbox, time_stride)  
    170 ;             domdef, box_plot[0:3], /MEMEINDICES 
    171 ;             ; force all_data read 
    172 ;             all_data = 1 
    173 ;          END  
    174 ;          ELSE: domdef, /MEMEINDICES 
    175 ;      ENDCASE  
    176 ;    ENDELSE  
    177 ;    ;box_plot = [ box_plot[0:3],vert_mean] 
    178 ;    print, jpt, time1, time2, hotyp 
    179118 
    180119; save time in common fld_att 
    … …  
    341280; Redefinition of time because it is updated in read_ncdf 
    342281; Useful for pltt routine 
     282 
    343283   IF hotyp NE '-' THEN BEGIN 
    344284      time = timearr.scale 
    … …  
    347287   ENDELSE 
    348288 
    349 ; if tkeavt, take log 
    350  
    351 ;    CASE cmd.var OF  
    352 ;       'votkeavt': BEGIN 
    353 ;          idx = where(fldr.data LE valmask/10) 
    354 ;          idxm = where(fldr.data EQ valmask) 
    355 ;          fldr.data(idx) = alog10(fldr.data(idx)) 
    356 ;          ; on T grid 
    357 ;          vargrid = 'T' 
    358 ;          fldr.data = (fldr.data+shift(fldr.data, 0, 0, 1))*.5 
    359 ;          fldr.data(idxm) = valmask 
    360 ;       END  
    361 ;       ELSE:  
    362 ;    ENDCASE   
    363  
    364289 
    365290; density projection 
     291 
    366292   IF splot EQ 1 THEN BEGIN 
    367293      IF cmd.timave EQ '1y' AND 1 THEN BEGIN  

  • branches/procs-read/nc_read.pro

    r83 r84  
    11;---------------------------   
    2 ; Reading ORCA netcdf files 
    3 ; EG 24/2/99 
     2; Reading any netcdf file 
    43;---------------------------   
    54FUNCTION 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 
    … …  
    3130   ENDCASE 
    3231 
    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 
    44 ;Rajout MK 
    45 ;Trouver une meilleure place 
    46 ;      ixminmesh = 0 
    47 ;      iyminmesh = 0 
    48 ;      ixmaxmesh = jpi-1 
    49 ;      iymaxmesh = jpj-1 
    50 ;Fin rajout 
    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 
    60  
    61  
    62 ; define reading boundaries 
    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  
    68  
    69 ; force izmaxmesh to lastz 
    70 ;   jpk = lastz+1 
    71  
    72 ; dealing with longitude periodicity 
    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 
    82  
    8332; local directory 
    84  
    8533   IF strpos(ncdf_db, ':') GE 1 THEN directory = (strsplit(ncdf_db, ':', /EXTRACT))[1] $ 
    8634    ELSE directory = ncdf_db 
    87  
    88 ; existence of file 
    89  
    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  
    100  
    101 ; ouverture fichier netCDF + contenu 
    102 ;   cdfid=ncdf_open(directory+file_name)  
    103 ;   inq=ncdf_inquire(cdfid) 
    104 ; que contient la variable 
    105 ;   varid = ncdf_varid(cdfid, var_name) 
    106  
    107    name_suff = '' 
    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) 
    119  
    120 ; test sur la dimension 
    121 ;   err_mess = '' 
    122 ;   field_dim = n_elements(varinq.dim) 
    123  
    124 ; get unlimited record variable 
    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, '     ' 
    19835 
    19936; Find the variable's attribute 
    … …  
    21956                        /TIMESTEP, /ADDSCL_BEFORE, TOUT = tout, /NOSTRUCT, /CONT_NOFILL, GRID = vargrid, $ 
    22057                        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 ' 
    22160 
     61   IF debug_w THEN print, '     ' 
     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    
    22265   IF debug_w THEN print,  '    zinvar=', zinvar 
    223    IF debug_w THEN print, '     ' 
    224 ;   IF debug_w THEN print, '     key=', key  
    225 ;   IF debug_w THEN print, '     jpitotal=', jpitotal 
    22666   IF debug_w THEN print, '     ixmindta,iymindta,izmindta  =', ixmindta,iymindta,izmindta 
    22767   IF debug_w THEN print, '     ixminmesh,iyminmesh=', ixminmesh,iyminmesh 
    22868   IF debug_w THEN print, '     ixmaxmesh,iymaxmesh=', ixmaxmesh,iymaxmesh  
    22969   IF debug_w THEN print, '     izminmesh,izmaxmesh=', izminmesh,izmaxmesh  
    230 ;   IF debug_w THEN print, '     ixmin,iymin=', ixmin,iymin  
    231 ;   IF debug_w THEN print, '     firsttps,lasttps=', firsttps, lasttps 
    23270   IF debug_w THEN print, '     jpt=', jpt 
    233    IF debug_w THEN print, '     key_shift in nc_read=', key_shift 
    234 ;   IF debug_w THEN print, '     key_yreverse, firsty, lasty',key_yreverse, firsty, lasty 
     71   IF debug_w THEN print, '     key_shift=', key_shift 
    23572   IF debug_w THEN print, '     key_yreverse=',key_yreverse  
    23673   IF debug_w THEN print, '     ' 
    237  
    238    ;lec_data = struct_data.tab 
    239    ;data_direc = 'xy' 
    240    ;field_dim = (size(struct.tab))[0] 
    241    field_dim = (size(lec_data))[0] 
    242    IF field_dim LE 0 THEN stop,  '  Something wrong happened in read_ncdf ' 
    243  
    244    ; build 3D atmospheric mask if needed and mask values 
    245 ;   IF mesh_type EQ 'atm' THEN BEGIN 
    246 ;      ;tmask = make_3d_atm_msk(directory+file_name, field_dim) 
    247 ;      CASE atmos_msk OF  
    248 ;         0: print, '          [atmos grid : take all points] ' ; take all points 
    249 ;         1: BEGIN 
    250 ;                                ; take ocean points only 
    251 ;            idx = where(tmask EQ 0) 
    252 ;            stop 
    253 ;            IF idx NE -1 THEN lec_data(idx) = 1.e20 
    254 ;            print, '          [atmos grid : take ocean points only] ' 
    255 ;         END 
    256 ;         2: BEGIN 
    257 ;                                ; take land points only 
    258 ;            idx = where(tmask GT 0) 
    259 ;            IF idx NE -1 THEN lec_data(idx) = 1.e20 
    260 ;            print, '          [atmos grid : take land points only] ' 
    261 ;         END 
    262 ;      ENDCASE  
    263 ;   ENDIF 
    26474    
    26575; Average data along vertical if needed and update some features 
    266 ; needed for plt (data_direc, name_suff)    
     76; needed for plt (data_direc, name_suff) 
     77   name_suff = '' 
     78   data_direc = '' 
    26779   update_data, TAB = lec_data, VNAME = var_name, BOXZOOM = boxzoom, ZUNITS = units_depth, $ 
    26880    ZINVAR = zinvar, SUFFIX = name_suff, D_DIREC = data_direc, $ 
    26981    TIME_1 = time_1, TIME_2 =  time_2, NO_MEAN = no_mean 
    27082   field_dim = (size(lec_data))[0] 
    271    ;stop 
    27283 
    273 ;   CASE n_elements(varinq.dim) OF 
    274 ;   CASE no_read_ncdf_varget OF 
    275 ;   CASE field_dim OF 
    276 ;      1: BEGIN 
    277 ;         print,  ' !!!! 1D case not implemented !!!!' 
    278 ;         print,  ' !!!! read_ncdf cannot handle 1D data !!!!' 
    279 ;         END 
    280 ; 
    281 ;      ;; fichier 2d : surface 
    282 ;      2: BEGIN 
    283 ;         print, '    Found ', var_name, ' (2D data)   from ', file_name 
    284 ;        ;@read_ncdf_varget 
    285 ;         ;lec_data = res         
    286 ;         data_direc = 'xy' 
    287 ;         niveau = 1 
    288 ;         END  
    289  
    290 ;     ;; fichier 3d : 2 cas = 1/ 2d espace + temps 2/ 3d espace 
    291 ;     3: BEGIN  
    292 ;         ;; si varinq.dim contient la dim infinie (no 3) -> temps 
    293 ;         dim_3 = '3d' 
    294 ;         IF nb_time GE 1 THEN dim_3 = '2d' 
    295 ;         IF jpt GE 1 THEN dim_3 = '2d' 
    296 ;         IF dim_3 EQ '2d' THEN BEGIN  
    297 ;         IF jpt EQ 1 THEN BEGIN 
    298 ; 
    299 ;            print, '    Found ', var_name, ' (3D data)', '   from ', file_name 
    300 ;            ;@read_ncdf_varget 
    301 ;            ;lec_data = res 
    302 ;            data_direc = 'xyz' 
    303 ; 
    304 ;            IF time_2 EQ time_1 THEN BEGIN  
    305 ;               print, '    Found ', var_name, ' (2D data - index ', strcompress(string(time_1)),') from ', file_name 
    306 ;               ;@read_ncdf_varget 
    307 ;               ;lec_data = res 
    308 ;               data_direc = 'xy' 
    309 ;            ENDIF ELSE BEGIN  
    310 ;               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 
    311 ;               ;IF debug_w THEN print,  '     x_offset et y_offset et time_1 :',  x_offset,  y_offset,  time_1 
    312 ;               ;@read_ncdf_varget 
    313 ;               ;lec_data = res 
    314 ;               data_direc = 'xyt'  
    315 ;            ENDELSE  
    316 ; 
    317 ;         ENDIF ELSE IF jpt GT 1 THEN BEGIN 
    318 ;         ENDELSE 
    319 ;      END  
    320 ; 
    321 ;      ;; fichier 4d : volume + temps 
    322 ;      4: BEGIN 
    323 ; 
    324 ;         IF time_2 EQ time_1 THEN BEGIN  
    325 ;            print, '    Found ', var_name, ' (3D data - index ', strcompress(string(time_1)),') from ', file_name 
    326 ;            ; read vertical levels 
    327 ;            IF debug_w THEN print, '   mesh_type= ',mesh_type  
    328 ;            IF mesh_type EQ 'atm' THEN BEGIN 
    329 ;               ncdf_diminq,cdfid,2, name_level, nb_level 
    330 ;               IF debug_w THEN print, 'name_level, nb_level = ', name_level, nb_level 
    331 ; 
    332 ;               dimidl = ncdf_dimid(cdfid, name_level) 
    333 ;               ncdf_diminq,cdfid,dimidl, name_level, nb_level 
    334 ; 
    335 ;               varidl = ncdf_varid(cdfid, name_level) 
    336 ;               ; make sure name_level is in hPa 
    337 ;               ncdf_attget, cdfid, varidl, 'units', val_unit 
    338 ;               IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar'  AND string(val_unit) NE 'mbar' THEN BEGIN  
    339 ;                  print, '  vertical levels coordinate not obvious = ', name_level 
    340 ;                  print, '  ... using <levels> ...' 
    341 ;                  varidl = ncdf_varid(cdfid, 'levels') 
    342 ;               ENDIF  
    343   
    344 ;               ncdf_varget, cdfid, varidl, gdept  
    345 ;               gdept = gdept 
    346 ;               gdepw = gdept 
    347 ;               e3t = shift(gdept, 1)-gdept 
    348 ;               e3t[0] = e3t[1] 
    349 ;               e3w = e3t 
    350 ;               jpk = nb_level 
    351 ;               tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) 
    352 ;            ENDIF  
    353 ;            IF debug_w THEN print,  ' going into @read_ncdf_varget '  
    354 ;            @read_ncdf_varget 
    355 ;            lec_data = res 
    356 ;            lec_data =  reform(lec_data,x_count, y_count, z_count)  
    357 ;            data_direc = 'xyz' 
    358 ;            IF debug_w THEN help, lec_data 
    359  
    360 ;           vertical average ? 
    361 ;            IF vert_switch ge 1 AND NOT keyword_set(no_mean) THEN BEGIN  
    362 ;               old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    363 ;               print, '      Average in vertical domain ', vert_type, vert_mean 
    364 ;               IF mesh_type EQ 'atm' THEN BEGIN 
    365 ;                  CASE atmos_msk OF  
    366 ;                     0: print, '          [atmos grid : take all points] ' ; take all points 
    367 ;                     1: BEGIN 
    368 ;                        ; take ocean points only 
    369 ;                        idx = where(tmask EQ 0) 
    370 ;                        lec_data(idx) = 1.e20 
    371 ;                        print, '          [atmos grid : take ocean points only] ' 
    372 ;                     END 
    373 ;                     2: BEGIN 
    374 ;                        ; take land points only 
    375 ;                        idx = where(tmask GT 0) 
    376 ;                        lec_data(idx) = 1.e20 
    377 ;                        print, '          [atmos grid : take land points only] ' 
    378 ;                     END 
    379 ;                  ENDCASE  
    380 ;                  CASE vert_type OF 
    381 ;                     'z': BEGIN ;  depth range 
    382 ;                        IF debug_w THEN print, ' performing average' 
    383 ;                        zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
    384 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    385 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
    386 ;                        ENDIF ELSE BEGIN 
    387 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    388 ;                        ENDELSE  
    389 ;                     END  
    390 ;                     ELSE: BEGIN ; levels range 
    391 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    392 ;                           zmean = lec_data(*, *, vert_mean[0]) 
    393 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
    394 ;                        ENDIF ELSE BEGIN 
    395 ;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) 
    396 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    397 ;                        ENDELSE  
    398 ;                     END  
    399 ;                  ENDCASE  
    400 ;;                  tmask = reform(tmask(*, *, 0), jpi, jpj) 
    401 ;               ENDIF ELSE BEGIN 
    402 ;;                  ; ocean = always mask 
    403 ;;                  idx = where(tmask EQ 0) 
    404 ;;                  lec_data(idx) = 1.e20 
    405 ;                  CASE vert_type OF 
    406 ;                     'z': BEGIN 
    407 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    408 ;                           zmean = lec_data ; right depth already selected 
    409 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    410 ;                        ENDIF ELSE BEGIN 
    411 ;                           zmean = moyenne(lec_data, 'z', boite = vert_mean) 
    412 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    413 ;                        ENDELSE  
    414 ;                     END  
    415 ;                     ELSE: BEGIN ; case level (zindex) 
    416 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    417 ;                           zmean = lec_data ; right depth already selected 
    418 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    419 ;                        ENDIF ELSE BEGIN 
    420 ;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) 
    421 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    422 ;                        ENDELSE  
    423 ;                     END  
    424 ;                  ENDCASE  
    425 ;               ENDELSE  
    426 ;               lec_data = zmean 
    427 ;               domdef, old_boite 
    428 ;               field_dim = field_dim - 1 
    429 ;               data_direc = 'xy' 
    430 ;               nzt = 1 
    431 ;               firstz = 0 
    432 ;               lastz = 0 
    433 ;            ENDIF  
    434 ; 
    435 ;         ENDIF ELSE BEGIN  
    436 ;            print, '    Reading ', var_name, ' (3D data time serie)',  strcompress(string(time_1)),'-', strtrim(string(time_2), 2), '   from ', file_name 
    437 ;            ; read vertical levels 
    438 ;            IF mesh_type EQ 'atm' THEN BEGIN 
    439 ;               ncdf_diminq,cdfid,2, name_level, nb_level 
    440  
    441                ; get name of vertical level from metadata 
    442 ;               file_grid_config = hom_def+'grid_config.def' 
    443 ;               spawn, 'grep -i "\ '+meshlec_type+' " '+file_grid_config, line 
    444 ;               line = strcompress(strtrim(line[0], 2))  
    445 ;               length = strlen(line) 
    446     
    447 ;               IF length EQ 0 THEN BEGIN 
    448 ;                  print, '    *** nc_read : define name_level from grid ', meshlec_type, $ 
    449 ;                   ' in file ', file_grid_config 
    450 ;                  name_level = '' 
    451 ;                  return, -1 
    452 ;               ENDIF ELSE BEGIN  
    453 ;                  argvar = strsplit(line, ' ', /EXTRACT)  
    454 ;                  name_level = argvar[1] 
    455 ;               ENDELSE 
    456 ;               dimidl = ncdf_dimid(cdfid, name_level) 
    457 ;               ncdf_diminq,cdfid,dimidl, name_level, nb_level 
    458 ; 
    459 ;               varidl = ncdf_varid(cdfid, name_level) 
    460 ;               ; make sure name_level is in hPa 
    461 ;               ncdf_attget, cdfid, varidl, 'units', val_unit 
    462 ;               IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar' AND string(val_unit) NE 'mbar' THEN BEGIN  
    463 ;                  print, '  vertical levels coordinate not obvious' 
    464 ;                  print, '  ... using <levels> ...' 
    465 ;                  varidl = ncdf_varid(cdfid, 'levels') 
    466 ;               ENDIF  
    467 ;               ncdf_varget, cdfid, varidl, gdept  
    468 ;               gdepw = gdept 
    469 ;               e3t = shift(gdept, 1)-gdept 
    470 ;               e3t[0] = e3t[1] 
    471 ;               e3w = e3t 
    472 ;               jpk = nb_level 
    473 ;               tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) 
    474 ;            ENDIF  
    475 ;            @read_ncdf_varget 
    476 ;            lec_data = res 
    477 ;            data_direc = 'xyzt' 
    478 ; 
    479 ;;           vertical average ? 
    480 ;            IF vert_switch ge 1  AND NOT keyword_set(no_mean) THEN BEGIN  
    481 ;               old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    482 ;               print, '      Average in vertical domain ', vert_type, vert_mean 
    483 ;               IF mesh_type EQ 'atm' THEN BEGIN 
    484 ;                  CASE atmos_msk OF  
    485 ;                     0: print, '          [take all points] ' ; take all points 
    486 ;                     1: BEGIN 
    487 ;                        ; take ocean points only 
    488 ;                        idx = where(tmask EQ 0) 
    489 ;                        lec_data(idx) = 1.e20 
    490 ;                        print, '          [take ocean points only] ' 
    491 ;                     END 
    492 ;                     2: BEGIN 
    493 ;                        ; take land points only 
    494 ;                        idx = where(tmask GT 0) 
    495 ;                        lec_data(idx) = 1.e20 
    496 ;                        print, '          [take land points only] ' 
    497 ;                     END 
    498 ;                  ENDCASE  
    499 ;                  CASE vert_type OF 
    500 ;                     'z': BEGIN ; levels 
    501 ;                         
    502 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    503 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' 
    504 ;                        ENDIF ELSE BEGIN 
    505 ;                           zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
    506 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean(0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']' 
    507 ;                        ENDELSE  
    508 ;                     END  
    509 ;                     ELSE: BEGIN ; levels 
    510 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    511 ;                           zmean = lec_data(*, *, vert_mean[0]) 
    512 ;                           name_suff = ' at '+vert_type+' '+strtrim(string(long(vert_mean(0))), 2)+' ['+strtrim(string(gdept(vert_mean(0))), 2)+' hPa]' 
    513 ;                        ENDIF ELSE BEGIN 
    514 ;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) 
    515 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    516 ;                        ENDELSE  
    517 ;                     END  
    518 ;                  ENDCASE  
    519 ;;                  tmask = reform(tmask(*, *, 0), jpi, jpj) 
    520 ;               ENDIF ELSE BEGIN 
    521 ;                  ; ocean = always mask 
    522 ;                  ; idx = where(tmask EQ 0) 
    523 ;                  ; lec_data(idx) = 1.e20 
    524 ;                  CASE vert_type OF 
    525 ;                     'z': BEGIN 
    526 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    527 ;                           zmean = lec_data 
    528 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    529 ;                        ENDIF ELSE BEGIN 
    530 ;                           zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) 
    531 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    532 ;                        ENDELSE  
    533 ;                     END  
    534 ;                     ELSE: BEGIN 
    535 ;                        IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 
    536 ;                           zmean = lec_data 
    537 ;                           name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' 
    538 ;                        ENDIF ELSE BEGIN 
    539 ;                           zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) 
    540 ;                           name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' 
    541 ;                        ENDELSE  
    542 ;                     END  
    543 ;                  ENDCASE  
    544 ;               ENDELSE  
    545 ;               lec_data = zmean 
    546 ;               domdef, old_boite 
    547 ;               field_dim = field_dim - 1 
    548 ;               data_direc = 'xyt' 
    549 ;               nzt = 1 
    550 ;               firstz = 0 
    551 ;               lastz = 0 
    552 ;            ENDIF  
    553 ;         ENDELSE  
    554 ;      END  
    555 ;     ;; erreur si dim > 4 
    556 ;      ELSE: BEGIN 
    557 ;         err_mess =  '  *** nc_read : ERROR dimension > 4' 
    558 ;         lec_data = -1.0 
    559 ;      END 
    560 ;   ENDCASE  
    561  
    562 ; Attribut du champ 
    563  
     84; Field attributes 
    56485   field = {name: '', data: lec_data, legend: '', units: '', origin: '', dim: 0, direc: data_direc} 
    56586   field.name = var_name 
    56687   field.origin = directory+file_name 
    567     
    568 ; 
    569 ; get long name  
    570 ;   result = '???' 
    571 ;   FOR i = 0, varinq.natts-1 DO BEGIN 
    572 ;      att_txt = ncdf_attname(cdfid, varid, i) 
    573 ;      IF att_txt EQ 'long_name' OR att_txt EQ 'title' THEN ncdf_attget, cdfid, varid, att_txt, result 
    574 ;   ENDFOR  
    575     
    576 ;   field.legend = strtrim(string(result),1)+name_suff 
    577  
    578 ;;;;   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 
    57988   field.legend = long_name+name_suff 
    580  
    581 ; scaling ? 
    582 ;   FOR i = 0, varinq.natts-1 DO BEGIN 
    583 ;      att_txt = ncdf_attname(cdfid, varid, i) 
    584 ;      IF att_txt EQ 'scale_factor' THEN BEGIN 
    585 ;         ncdf_attget, cdfid, varid, att_txt, valscale 
    586 ;         lec_data = lec_data*valscale 
    587 ;      ENDIF  
    588 ;   ENDFOR  
    589  
    590  
    591 ; get unit 
    592 ; if it exists 
    593 ;  isunits = ncdf_attinq(cdfid,  varid,  'units') 
    594 ;  IF isunits.datatype NE 'UNKNOWN' THEN BEGIN 
    595 ;     ncdf_attget, cdfid, varid, 'units', result 
    596 ;     field.units = strtrim(string(result),1) 
    597 ;  ENDIF ELSE BEGIN 
    598 ;     print,  'No units for the variable ',  field.name 
    599 ;     print,  ' Or units were forgotten in the file !' 
    600 ;  ENDELSE 
    60189   field.units = units 
    60290   field.dim = field_dim  
    60391 
    604 ;  get valmask (might need valmask = float(string(valmask)) 
     92; get valmask (might need valmask = float(string(valmask)) 
    60593 
    606 ;   valmask = 1.e20 
    607 ;   FOR i = 0, varinq.natts-1 DO BEGIN 
    608 ;      att_txt = ncdf_attname(cdfid, varid, i) 
    609 ;      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 
    610 ;   ENDFOR  
    611  
     94; valmask = 1.e20 
    61295   IF size(missing_value,  /TYPE) EQ 4 OR size(missing_value,  /TYPE) EQ 5 THEN BEGIN 
    61396      valmask = missing_value 
    … …  
    640123   print, '      = ',field.legend, '    [min/max = ',minf , maxf,'  ', field.units,' - masked values = ',valmask, chardim, ']' 
    641124 
    642 ;   ncdf_close, cdfid 
    643  
    644    IF keyword_set(key_yreverse) THEN print, '    key_yreverse active : ',  key_yreverse 
    645  
    646125   IF debug_w THEN print,  '  ...EXIT nc_read' 
    647126   IF debug_w THEN print,  '  ' 

  • branches/procs-read/update_data.pro

    r83 r84  
    77   IF debug_w THEN print, '   ENTER update_data...' 
    88    
     9   ; Perform mean over vertical levels or extract one vertical level ? 
    910   dim = (size(tab))[0] 
    10    ; Perform mean over vertical levels or extract one vertical level ? 
    1111   perform_mean = vert_switch GE 1 AND zinvar EQ 1 AND NOT keyword_set(no_mean) 
    1212 
    … …  
    2222         print, '    Found ', vname, ' (2D data) from file'        
    2323         d_direc = 'xy' 
    24          ;;; niveau = 1 ;;; Ca sert à quoi cette variable ? 
    2524         IF perform_mean EQ 1 THEN BEGIN 
    2625            ; vert_mean[0] can be different from vert_mean[1] to select one level or layer 
    … …  
    5857               ENDCASE 
    5958               tab = temporary(zmean) 
    60                ;domdef, old_boite ;; a quoi ca sert ??? 
    6159               d_direc = 'xy' 
    6260            ENDIF 
    … …  
    7472            d_direc = 'xyzt' 
    7573            IF perform_mean EQ 1 THEN BEGIN 
    76                ;old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    7774               print, '      Average in vertical domain ', vert_type, vert_mean 
    78                ;IF vert_mean[0] EQ vert_mean[1] THEN stop, report('vert_mean[0] cannot be equal to vert_mean[1]') 
    7975               CASE vert_type OF 
    8076                  'z': BEGIN    
    … …  
    8884               ENDCASE 
    8985               tab = temporary(zmean) 
    90                ;domdef, old_boite 
    9186               d_direc = 'xyt' 
    9287            ENDIF 
Note: See TracChangeset for help on using the changeset viewer.