  function tp2_new
%      *   program tp2_new.m
%      *   tp2_new -- "NetCDF Toolbox for Matlab-5" code.
%      *   tp2_new (no argument) is a code  that builds cdf file, lists
%      *   itself, then displays its variables.
 
% ---------------------------- DEFINE THE FILE --------------------------- %
addpath /usr/people/valery/netcdf/
addpath /usr/people/valery/netcdf/ncfiles
addpath /usr/people/valery/netcdf/nctype
addpath /usr/people/valery/netcdf/ncutility
addpath /whoinet/deltorchio/p2/valery/TP_netCDF

ncquiet                                              

nc = netcdf('tp2_new.nc', 'clobber');                
%% Global attributes:
nc.description = 'The output of the QG model configured in T/P diamond'; 
nc.authors = 'D. McGillicuddy & V. Kosnyrev';
nc.date = '20 July, 2001';
%% Dimensions:
% define dimensions
nc('X_BATS') = 31;
nc('Y_BATS') = 25;
nc('TIME') = 2151;%% (record dimension)
nc('DEPTH') = 6;
%% Variables and attributes:
%  define variables
nc{'lat'} = ncfloat('Y_BATS','X_BATS');                      
nc{'lat'}.missing_value = ncfloat(-1.e+34);
nc{'lat'}.y_axis_rotation = 'rotated 64 deg cckw';
nc{'lat'}.units = 'degrees_north';
nc{'lat'}.point_spacing = 'even';
nc{'lat'}.long_name = 'latitude';

nc{'lon'} = ncfloat('Y_BATS','X_BATS');
nc{'lon'}.missing_value = ncfloat(-1.e+34);
nc{'lon'}.x_axis_rotation = 'rotated 64 deg cckw';
nc{'lon'}.units = 'degrees_east';                    
nc{'lon'}.point_spacing = 'even';                        
nc{'lon'}.long_name = 'longitude';

nc{'streamf_sim'} = ncfloat('TIME','DEPTH','Y_BATS','X_BATS');
nc{'streamf_sim'}.missing_value = ncfloat(-1.e+34);
nc{'streamf_sim'}.units = 'meters^2/second';
nc{'streamf_sim'}.long_name = 'predicted streamfunction';

nc{'vort_sim'} = ncfloat('TIME','DEPTH','Y_BATS','X_BATS');
nc{'vort_sim'}.missing_value = ncfloat(-1.0e+34);
nc{'vort_sim'}.units = 'second^-1';
nc{'vort_sim'}.long_name = 'predicted vorticity';

nc{'DEPTH'} = ncdouble('DEPTH');
nc{'DEPTH'}.positive = 'down';
nc{'DEPTH'}.points_spacing = 'uneven';
nc{'DEPTH'}.units = 'meters';
nc{'DEPTH'}.long_name = 'depth streamfunction defined';

nc{'TIME'}=ncdouble('TIME');
nc{'TIME'}.time_origin = '26-OCT-1992';
nc{'TIME'}.units = 'DAYS since 1992-10-27 00:00:00';

nc{'lat'} = 'latitude';                        
nc{'lon'} = 'longitude';
nc{'streamf_sim'} = 'simf';
nc{'vort_sim'} = 'simv';
nc{'TIME'} = 'time';
nc{'DEPTH'} = 'depth';
% ---------------------------- LOAD DATA ----------------------------------%
% ---------------------------- STORE DATA ---------------------------- %
ti = 2151;
de = 6;
xb=31;
yb=25;

% prepare large scale grid  
      [x3,y3]=meshgrid(.85:.57627:34.85,1:.68085:33);
%transform grid to the geographycal coordinate system 
      x3=-70+(x3-.85)*10/34;
      y3=28.25+(y3-1)*8/32;
% rotation angle
      f=pi*64.0/180;
% rotated coordinate system origin
      yn=10.2728;
      xn=23.4591;
% build meshgrid
      xi=ones(25,31);
      yi=ones(25,31);
      for k=1:1:25
      for l=1:1:31
      xi(k,l)=(xn)+(-(k-1)*sin(f)+(l-1)*cos(f))*0.366;
      yi(k,l)=(yn)+((k-1)*cos(f)+(l-1)*sin(f))*0.366;
      end
      end
% T/P meshgrid     
      longitude = -70+(xi-.85)*10/34;
      latitude=28.25+(yi-1)*8/32;
% load simulated  streamfunction

im=ti*de;

fid =fopen('dh3a15b5f_2001.vbc','r')
a=fscanf(fid,'%f',[775 im])';
fclose(fid);
for i=[1:ti]
kk=(i-1)*6;
for k=[1:de]
kk=kk+1;
simf(i,k,:,:)=reshape(a(kk,:),31,25)';
end
end



% load simulated  vorticity
fid =fopen('dh3a15b5f_2001.voc','r')
a=fscanf(fid,'%f',[775 im])';
fclose(fid);
for i=[1:ti]
kk=(i-1)*de;
for k=[1:de]
kk=kk+1;
simv(i,k,:,:)=reshape(a(kk,:),31,25)';
end
end
clear a

    time = 1:2151;
    depth=[ 100 400  800 1200 2000 3550];   

% load variables
nc{'lat'}(:) = latitude;
nc{'lon'}(:) = longitude;
% rescale to the dimensional units (m**2/sec ** 10^-3)
nc{'streamf_sim'}(:) = simf;
% rescale to the dimensional units (1/sec **10^6)
nc{'vort_sim'}(:) = simv;
nc{'DEPTH'}(:) = depth;
nc{'TIME'}(:)  = time;

nc = close(nc);                                      % Close the file.

% ---------------------------- RECALL THE DATA --------------------------- %

nc = netcdf('tp2_new.nc', 'nowrite');                 % Open NetCDF file.
description = nc.description(:)                      % Global attribute.
variables = var(nc)                                 % Get variable data.
for i = 1:length(variables)
   disp([name(variables{i}) ' =']), disp(' ')
   disp(variables{i})
   disp([size(variables{i})])
   disp([att(variables{i})])
   disp([dim(variables{i})])
end
nc = close(nc);                                      % Close the file.

% --------------------------------- DONE --------------------------------- %
