function unit = vs_trim2nc_test
%vs_trim2nc_test  unit test for vs_trim2nc, incl dry point mask
%
% Checks conversion of primary variables from NEFIS to netCDF using 
% Delft3D input matrices. The test valyes are hard-coded to match 
% those in the mdf and associated files.
%
% nan         `---+----@----+----@----+----@----+----@        
%        n=5   ::::::::|:::::::::|:::::::::|:::::::::|        
%   nan        :::x::::+::::x::::+::::x::::+::::x::::|        
%       - - -  ::::::::|:::::::::|:::::::::|:::::::::|        
% 40          `---+----o----+----o----+----o----+----@        
%        n=4   ::::::::|         |         |:::::::::|        
%   35         :::x::::+    x    +    x    +::::x::::|        
%       - - -  ::::::::|         |         |:::::::::|        
% 30          `---+----o----+----o----+----o----+----@        
%        n=3   ::::::::|         |         |:::::::::|        
%   25         :::x::::+    x    +    x    +::::x::::|        
%       - - -  ::::::::|         |         |:::::::::|        
% 20          `---+----o----+----o----+----o----+----@        
%        n=2   ::::::::|         |:::::::::|:::::::::|        
%   15         :::x::::+    x    +::::x::::+::::x::::|        
%       - - -  ::::::::|         |:::::::::|:::::::::|        
% 10          `---+----o----+----o----+----o----+----@        
%        n=1   ::::::::|:::::::::|:::::::::|:::::::::|        
%   nan        :::x::::+::::x::::+::::x::::+::::x::::|        
% nan   - - -  ::::::::|:::::::::|:::::::::|:::::::::|        
%                      `      `       `       `       ‘      
% =      y     .       .       .       .       .      
% = grid_y     .  m=1  .  m=2  .  m=3  .  m=6  .      
%              .       .       .       .       .      
%  grid_x = nan       10        30        50        nan
%       x =      nan        20        40       nan
%%
%See also: VS_TRIM2NC, NETCDF, NCWRITESCHEMA, DELFT3D, DELFT3D_GRID_IMAGE
%          nc_cf_bounds2cor, nc_cf_cor2bounds

% TO DO: test also for missing variables in trim file (keywords SMhydr, SMderv, SMproc not all Y)

rehash toolbox % update textpad vs_trim2nc

OPT = vs_trim2nc;

RUNIDs = {'2x3','2x3spherical','2x3z'};
epsg   = [28992 nan 28992];
OK     = [];

for iRUNID=1:3

   RUNID = RUNIDs{iRUNID};
   
   clear D
      vs_trim2nc([pwd,filesep,'trim-',RUNID,'.dat'],'epsg',epsg(iRUNID),'var',{OPT.var{:},'pea','dpeadt'},'Format','classic') % 64bit
   D = nc2struct([pwd,filesep,'trim-',RUNID,'.nc' ],'rename',{{'velocity_x','velocity_y'},{'u','v'}});
   %structshow(D)
   
   fid = fopen(['trim-',RUNID,'.cdl'],'w');
   nc_dump(    ['trim-',RUNID,'.nc' ],fid);
   fclose(fid);
   
   n = 2:3;
   m = 2:4;
   
   D0.x     = [                           nan nan nan nan nan;   nan 20 20 20 nan;    nan NaN 40 40 nan;    nan nan nan nan nan];
   D0.y     = [                           nan nan nan nan nan;   nan 15 25 35 nan;    nan NaN 25 35 nan;    nan nan nan nan nan];
   % for xcor we added TWO rows and columns (at the start and at the end)
   % as nc_cf_cor2bounds removes one, so then it is same size as xcen + 2*dummy row/cols
   D0.xcor = addrowcol([ 10  10  10  10; 30 30 30 30;nan  50 50 50],[-1 1],[-1 1],nan);
   D0.ycor = addrowcol([ 10  20  30  40; 10 20 30 40;nan  20 30 40],[-1 1],[-1 1],nan);

   if iRUNID==1 | iRUNID==3
     D0.x    = D0.x*10;    
     D0.y    = D0.y*10;    
     D0.xcor = D0.xcor*10;    
     D0.ycor = D0.ycor*10;    
   elseif iRUNID==2
     D.x      = D.longitude   ; 
     D.y      = D.latitude    ;
     D.grid_x = D.grid_longitude;
     D.grid_y = D.grid_latitude ;
   end
   
   OK(end+1) = all(D.k==[1 2 3 4]');
   if iRUNID<3
   OK(end+1) = all((D.Layer       - [-0.05 -0.2 -0.45 -0.8   ]')< 1e-7);
   OK(end+1) = all((D.LayerInterf - [ 0    -0.1 -0.3  -0.6 -1]')< 1e-7);
   OK(end+1) = isequalwithequalnans(D.time,[datenum(1970,0,0,0,0:100:600,0) - datenum(1970,0,0)]');
   else
   OK(end+1) = all(D.Layer'       + corner2center((-(10-cumsum([0 1 2 3 4]))*12))< 1e-7);
   OK(end+1) = all(D.LayerInterf' + (-(10-cumsum([0 1 2 3 4]))*12)               < 1e-7);
   OK(end+1) = isequalwithequalnans(D.time,[datenum(1970,0,0,0,0,0) - datenum(1970,0,0)]');
   end
   
   OK(end+1) = isequalwithequalnans(D.x,D0.x);
   OK(end+1) = isequalwithequalnans(D.y,D0.y);
   
   OK(end+1) = isequalwithequalnans(D.area,grid_area(nc_cf_bounds2cor(D.grid_x),nc_cf_bounds2cor(D.grid_y))); % note NEFIS is wrong here
   
   OK(end+1) = isequalwithequalnans(nc_cf_cor2bounds(D0.xcor),D.grid_x);
   OK(end+1) = isequalwithequalnans(nc_cf_cor2bounds(D0.ycor),D.grid_y);
   OK(end+1) = isequalwithequalnans(D0.xcor ,nc_cf_bounds2cor(D.grid_x));
   OK(end+1) = isequalwithequalnans(D0.ycor ,nc_cf_bounds2cor(D.grid_y));
   if iRUNID<3
   OK(end+1) = all(all(abs(D.u(1,1,n,m)-.06  )<1e-4 | isnan(D.u(1,1,n,m))));
   OK(end+1) = all(all(abs(D.u(1,2,n,m)-.05  )<1e-4 | isnan(D.u(1,2,n,m))));
   OK(end+1) = all(all(abs(D.u(1,3,n,m)-.04  )<1e-4 | isnan(D.u(1,3,n,m))));
   OK(end+1) = all(all(abs(D.u(1,4,n,m)-.03  )<1e-4 | isnan(D.u(1,4,n,m))));
   
   OK(end+1) = all(all(abs(D.v(1,1,n,m)-.08  )<1e-4 | isnan(D.v(1,1,n,m))));
   OK(end+1) = all(all(abs(D.v(1,2,n,m)-.0667)<1e-3 | isnan(D.v(1,2,n,m))));
   OK(end+1) = all(all(abs(D.v(1,3,n,m)-.0533)<1e-3 | isnan(D.v(1,3,n,m))));
   OK(end+1) = all(all(abs(D.v(1,4,n,m)-.04  )<1e-4 | isnan(D.v(1,4,n,m))));
   
   for k=1:D.k(end)
   OK(end+1) = isequalwithequalnans(permute(D.temperature(1,  k,n,m),[  3 4 1 2]),[ 10  10  10;NaN  10  10]+k); % mixes to 11
   OK(end+1) = isequalwithequalnans(permute(D.salinity   (1,  k,n,m),[  3 4 1 2]),[ 30  30  30;NaN  30  30]+k); % mixes to 32.3,  not (.1*31 +.2*32+.3*33+.4*34) = 33
   OK(end+1) = isequalwithequalnans(permute(D.suspsedconc(1,1,k,n,m),[4 5 1 2 3]),[100 100 100;NaN 100 100]+k);
   OK(end+1) = isequalwithequalnans(permute(D.suspsedconc(1,2,k,n,m),[4 5 1 2 3]),[200 200 200;NaN 200 200]+k);
   end
   
   %% PEA
   %  with full mixing, dpeadt at t=2 (central difference from timestep 1 to 3)
   %  is identical to PEA at t=1 dvided by 2*dt
   
   % pea_simpson_et_al_1990(z,permute(D.density(1,:,:,:),[2 3 4 1]))
   OK(end+1) = all(all(abs(permute(D.dpeadt(2,:,:),[2 3 1]) - permute(-D.pea(1,:,:),[2 3 1])./200./60)< 1e-5 | isnan(permute(D.dpeadt(2,:,:),[2 3 1]))));
   
   % clc;disp('salinity');for t=1:7;for k=1:4;disp(permute(D.salinity(t,k,:,:),[3 4 1 2]));end;disp('-----------');pausedisp;end
   % clc;disp('pea')     ;for t=1:7;          disp(permute(D.dpeadt  (t,:,:)  ,[2 3 1  ]));    disp('-----------');pausedisp;end
   end

end

unit = (OK);
   
