function [ pert_grid, NRANDSEED, N_force_pert, N_progn_pert,   ...
           Pert_rseed, force_pert_ntrmdt, progn_pert_ntrmdt ]  ...
    = read_pert_ldas_rst( fname )

% matlab reader for LDASsa perturbation restart files
%
%  Rolf Reichle,    21 Mar 2015
%
% ------------------------------------------------------------------

int_precision   = 'int32';      % precision of fortran tag
float_precision = 'float32';    % precision of data in input file

disp(['reading from ', fname])

ifp      = fopen( fname, 'r', 'b' );


% first record: perturbations grid 
  
fortran_tag         = fread( ifp,  1, int_precision );
pert_grid.gridtype  = fread( ifp, 40, 'uint8=>char' ); 
pert_grid.ind_base  = fread( ifp,  1, int_precision );
pert_grid.i_dir     = fread( ifp,  1, int_precision );    
pert_grid.j_dir     = fread( ifp,  1, int_precision );    
pert_grid.N_lon     = fread( ifp,  1, int_precision );    
pert_grid.N_lat     = fread( ifp,  1, int_precision );     
pert_grid.i_offg    = fread( ifp,  1, int_precision );   
pert_grid.j_offg    = fread( ifp,  1, int_precision );  
pert_grid.ll_lon    = fread( ifp,  1, float_precision );   
pert_grid.ll_lat    = fread( ifp,  1, float_precision );  
pert_grid.ur_lon    = fread( ifp,  1, float_precision );   
pert_grid.ur_lat    = fread( ifp,  1, float_precision );  
pert_grid.dlon      = fread( ifp,  1, float_precision );     
pert_grid.dlat      = fread( ifp,  1, float_precision );         
fortran_tag         = fread( ifp,  1, int_precision );

pert_grid.gridtype  = deblank(pert_grid.gridtype');
 

% second record: assorted integers

fortran_tag         = fread( ifp,  1, int_precision );
NRANDSEED           = fread( ifp,  1, int_precision );
N_force_pert        = fread( ifp,  1, int_precision );
N_progn_pert        = fread( ifp,  1, int_precision );
fortran_tag         = fread( ifp,  1, int_precision );


% third record: random seed

fortran_tag         = fread( ifp,  1, int_precision  );
Pert_rseed          = fread( ifp, NRANDSEED, int_precision );
fortran_tag         = fread( ifp,  1, int_precision  );


% forcing perturbations (one record per gridded perturbation field)

force_pert_ntrmdt   = NaN( pert_grid.N_lon, pert_grid.N_lat, N_force_pert );

for i=1:N_force_pert
  
  fortran_tag       = fread( ifp,  1, int_precision  );
  dummy             = fread( ifp, pert_grid.N_lat*pert_grid.N_lon, float_precision );
  fotran_tag        = fread( ifp,  1, int_precision  );

  force_pert_ntrmdt(:,:,i) = reshape( dummy, pert_grid.N_lon, pert_grid.N_lat );

end

% prognostics perturbations (one record per gridded perturbation field)

progn_pert_ntrmdt   = NaN( pert_grid.N_lon, pert_grid.N_lat, N_progn_pert );

for i=1:N_progn_pert
  
  fortran_tag       = fread( ifp,  1, int_precision  );
  dummy             = fread( ifp, pert_grid.N_lat*pert_grid.N_lon, float_precision );
  fotran_tag        = fread( ifp,  1, int_precision  );

  progn_pert_ntrmdt(:,:,i) = reshape( dummy, pert_grid.N_lon, pert_grid.N_lat );

end

% ======= EOF ==========================================================

