DoMORE-2 draft
authorA.M. Thurnherr <athurnherr@yahoo.com>
Fri, 04 Sep 2015 06:37:43 +0000
changeset 17 f5a63c03d9c8
parent 16 54126cc9bb4a
child 18 6e7c8d592f7f
DoMORE-2
HISTORY
begin_processing_step.m
edit_data.m
getinv.m
ladcp2cdf.m
loadctd.m
loadrdi.m
loadsadcp.m
plotraw.m
process_cast.m
--- a/HISTORY	Wed Apr 15 09:55:34 2015 +0000
+++ b/HISTORY	Fri Sep 04 06:37:43 2015 +0000
@@ -1,9 +1,9 @@
 ======================================================================
                     H I S T O R Y 
                     doc: Tue Aug 16 11:41:46 2005
-                    dlm: Wed Apr 15 08:56:15 2015
+                    dlm: Mon Jul 27 15:50:13 2015
                     (c) 2005 A.M. Thurnherr
-                    uE-Info: 254 49 NIL 0 0 72 3 2 4 NIL ofnI
+                    uE-Info: 276 37 NIL 0 0 72 3 2 4 NIL ofnI
 ======================================================================
 
 NB: CHANGE VERSION IN [default.m] BEFORE UPLOADING
@@ -252,3 +252,25 @@
 	Apr 15, 2015:
 		- modified ambiguity-velocity warning message in [loadrdi.m]
 		  following a suggestion of Diana Cardoso
+		- implemented fixes for bugs reported by Achim Randelhoff via
+  		  yahoo email on 04/23/2014 [getinv.m] [plotraw.m]
+	Apr 22, 2015:
+		- modified [begin_processing_step.m] to evaluate the expression
+		  provided as an argument before re-loading [set_cast_params.m]
+	May  8, 2015:
+		- folded in DoMORE-1 [getinv.m] with changes to dragfac
+		  and smallfac
+	May 12, 2015:
+		- finally removed step 16 from [process_cast.m]
+	May 27, 2015:
+		- updated diagnostic mesages in [loadctd.m] [loadrdi.m]
+	May 28, 2015:
+		- added error message to [loadctd.m]
+	Jun 12, 2015:
+		- set p.lat and p.lon if only SADCP gps data are available
+	Jul 21, 2015:
+		- [edit_data.m] made bin masking more permissive
+	Jul 27, 2015:
+		- replaced [ladcp2cdf.m] with homegrown alternative that is
+		  *not* consistent with old format; Diana Cardoso and EF are
+		  working on replacement code
--- a/begin_processing_step.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/begin_processing_step.m	Fri Sep 04 06:37:43 2015 +0000
@@ -1,9 +1,9 @@
 %======================================================================
 %                    B E G I N _ P R O C E S S I N G _ S T E P . M 
 %                    doc: Fri Jun 25 16:13:41 2004
-%                    dlm: Thu Jun 26 13:47:02 2008
+%                    dlm: Wed Apr 22 08:10:18 2015
 %                    (c) 2004 ladcp@
-%                    uE-Info: 13 50 NIL 0 0 72 2 2 8 NIL ofnI
+%                    uE-Info: 15 58 NIL 0 0 72 2 2 8 NIL ofnI
 %======================================================================
 
 % start new processing step (in [process_cast.m])
@@ -11,6 +11,8 @@
 % HISTORY:
 %   Jun 25, 2004: - created
 %   Jun 26, 2008: - BUG: typo related to eval_expr
+%   Apr 22, 2015: - added evaluation of eval_expr before re-loading set_cast_params.m
+%		    to allow setting of processing_version
 
 msg = sprintf('#################### step %d: %s ',pcs.cur_step,pcs.step_name);
 while length(msg)<70, msg = [msg '#']; end
@@ -19,10 +21,14 @@
   save_pcs = pcs; % save state
   disp(sprintf('LOADING CHECKPOINT %s_%d',f.checkpoints,pcs.cur_step-1));
   load(sprintf('%s_%d',f.checkpoints,pcs.cur_step-1));
+  if ~isempty(save_pcs.eval_expr)
+    disp(sprintf('EVALUATING EXPRESSION <%s>...',save_pcs.eval_expr));
+    eval(save_pcs.eval_expr);
+  end
   disp('RE-LOADING PER-CAST PARAMETERS');
   set_cast_params;
   if ~isempty(save_pcs.eval_expr)
-    disp(sprintf('EVALUATING EXPRESSION <%s>...',save_pcs.eval_expr));
+    disp(sprintf('RE-EVALUATING EXPRESSION <%s>...',save_pcs.eval_expr));
     eval(save_pcs.eval_expr);
   end
   pcs = save_pcs; % restore state
--- a/edit_data.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/edit_data.m	Fri Sep 04 06:37:43 2015 +0000
@@ -19,13 +19,14 @@
 %  Jan 23, 2015: - BUG: automatic zero blanking editing had a typo
 %		 - BUG: automatic zero blanking editing did not work with DL-only data
 %		 - added p.edit_{dn,up}_bad_ensembles
+%  Jul 21, 2015: - made bin-masking more permissive (allow indices > #bins)
 
 %======================================================================
 %                    E D I T _ D A T A . M 
 %                    doc: Sat Jul  3 17:13:05 2004
-%                    dlm: Fri Jan 23 21:40:14 2015
+%                    dlm: Tue Jul 21 11:39:18 2015
 %                    (c) 2004 A.M. Thurnherr
-%                    uE-Info: 93 23 NIL 0 0 72 0 2 8 NIL ofnI
+%                    uE-Info: 126 23 NIL 0 0 72 0 2 8 NIL ofnI
 %======================================================================
 
 %----------------------------------------------------------------------
@@ -113,16 +114,20 @@
   nbad = 0;
   if length(d.zu) > 0
     for bi=1:length(p.edit_mask_up_bins)
-      bn = length(d.zu)+1 - p.edit_mask_up_bins(bi);
-      nbad = nbad + length(find(isfinite(d.weight(bn,:))));
-      d.weight(bn,:) = NaN; d.ts_edited(bn,:) = NaN;
+      if p.edit_mask_up_bins(bi)<=p.nbin_u
+         bn = length(d.zu)+1 - p.edit_mask_up_bins(bi);
+	 nbad = nbad + length(find(isfinite(d.weight(bn,:))));
+	 d.weight(bn,:) = NaN; d.ts_edited(bn,:) = NaN;
+      end
     end
   end
   if length(d.zd) > 0
     for bi=1:length(p.edit_mask_dn_bins)
-      bn = length(d.zu) + p.edit_mask_dn_bins(bi);
-      nbad = nbad + length(find(isfinite(d.weight(bn,:))));
-      d.weight(bn,:) = NaN; d.ts_edited(bn,:) = NaN;
+      if p.edit_mask_dn_bins(bi)<=p.nbin_d
+        bn = length(d.zu) + p.edit_mask_dn_bins(bi);
+	nbad = nbad + length(find(isfinite(d.weight(bn,:))));
+	d.weight(bn,:) = NaN; d.ts_edited(bn,:) = NaN;
+      end
     end
   end
 
--- a/getinv.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/getinv.m	Fri Sep 04 06:37:43 2015 +0000
@@ -8,9 +8,9 @@
 %======================================================================
 %                    G E T I N V . M 
 %                    doc: Thu Jun 17 15:36:21 2004
-%                    dlm: Fri Jan  6 11:57:45 2012
+%                    dlm: Thu Jul 23 14:06:09 2015
 %                    (c) 2004 ladcp@
-%                    uE-Info: 31 50 NIL 0 0 72 2 2 8 NIL ofnI
+%                    uE-Info: 292 1 NIL 0 0 72 0 2 8 NIL ofnI
 %======================================================================
 
 % CHANGE HISTORY:
@@ -29,6 +29,9 @@
 % Jan  6, 2012: - lesqchol removed because it does not deal with
 %		  nearly singular matrices gracefully
 %		- BUG: replaced all imag() by new imagnan()
+% Jul 28, 2014: - modified how to specify smallfac
+% Aug  9, 2014: - modified how to specify dragfac (+ve for fixed; -ve for scaled Martin's default)
+% Jul 23, 2015: - commented on bug below (#!#)
 
 if nargin<5, iplot=0; end
 
@@ -46,16 +49,31 @@
 % how strong do you believe in the simple wire dynamics (not very good yet)
 %ps=setdefv(ps,'dragfac',tanh(p.maxdepth/3000)*0.2);
 ps=setdefv(ps,'dragfac',0);
+
+% negative values weigh Martin's original default
+if (ps.dragfac < 0)
+  ps.dragfac = -ps.dragfac * tanh(p.maxdepth/3000)*0.2;
+end
+
 % how much do you want to smooth the ocean profile
 ps=setdefv(ps,'smoofac',0);
-% how much do you want to smooth the ocean profile
-imax=ceil(p.maxdepth/500);
-smallfac=[0 0];
-for i=1:imax
- smallfac(i,1)=i;
- smallfac(i,2)=0.02/(1+abs(i-(imax/2)))*tanh(p.maxdepth/3000);
-end
-ps=setdefv(ps,'smallfac',smallfac);
+
+% how much do you want to restrict large shears on short vertical wavelengths
+%  	ps.smallfac(1) is the vertical-wavelength filter cutoff
+%	ps.smallfac(2) is the filter strength (larger implies more filtering)
+
+%ps=setdefv(ps,'smallfac',[500 0.02]);		% Martin's defaults
+%ps=setdefv(ps,'smallfac',[200 0.02]);		% works reasonably well for DoMORE1 data
+
+ps=setdefv(ps,'smallfac',[99999 0]);		% default: disable
+
+    imax = ceil(p.maxdepth/ps.smallfac(1));
+    smallfac = [0 0];
+    for i=1:imax
+	smallfac(i,1) = i;
+	smallfac(i,2) = ps.smallfac(2)/(1+abs(i-(imax/2)))*tanh(p.maxdepth/3000);
+    end
+    
 % do you want up/down cast seperately then set to 1
 ps=setdefv(ps,'down_up',1);
 % decide how to weight data
@@ -269,6 +287,9 @@
 
 [jmax,jbott]=max(izv);
 
+%#!# the following line is wrong, as izv>(ps.dz*(jmax-1) is never true;
+%#!# the code is supposed to remove the "last data bin" (whatever that means)
+%#!# but it does not do anything
 ii=find(izv<(ps.dz) | izv>(ps.dz*(jmax-1)) | izv>(zbottom-ps.dz));
 wm(ii)=0;
 
@@ -467,8 +488,8 @@
 
 
 %### small deep ocean velocity
-if sum(ps.smallfac(:,2))>0
- [A2,A1,d]=lainsmal(A2,A1,d,ps.smallfac);
+if sum(smallfac(:,2))>0
+ [A2,A1,d]=lainsmal(A2,A1,d,smallfac);
 end
 
 de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];
--- a/ladcp2cdf.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/ladcp2cdf.m	Fri Sep 04 06:37:43 2015 +0000
@@ -1,176 +1,250 @@
-function [] = ladcp2cdf(fname,dr_struct,da,p,ps,f,att);
-% function [] = ladcp2cdf(fname,dr_struct,da,p,ps,f,att);
-%
-% function to save LADCP data into a netcdf file for MatLab version 2012a
-%
-% input  :	fname		- output filename
-%			dr_struct	- main inversion results (velocity profiles)
-%			da..att	- arbitrary metadata structures
-%
-% Created By:   Diana Cardoso, Bedford Institute of Oceangraphy
-%               Diana.Cardoso@dfo-mpo.gc.ca
-% Description:  Based on LDEO software to Process LADCP, version IX.8,
-%               script ladcp2cdf.m version 0.1	last change 08.03.2002. 
-%               maintained by A.M. Thurnherr and downloaded from:
-%       http://www.ldeo.columbia.edu/cgi-bin/ladcp-cgi-bin/hgwebdir.cgi
-%       The function ladcp2cdf was changed to run with the the Matlab
-%       version 2012, which now supports netcdf.
-
 %======================================================================
-%                    L A D C P 2 C D F . M 
-%                    doc: Thu Aug 15 10:52:55 2013
-%                    dlm: Wed Aug 28 12:31:16 2013
-%                    (c) 2013 A.M. Thurnherr, from code contributed by D. Cardoso
-%                    uE-Info: 99 0 NIL 0 0 72 0 2 8 NIL ofnI
+%                    L A D C P 2 C D F _ V 2 . M 
+%                    doc: Fri Jul 24 09:31:43 2015
+%                    dlm: Mon Jul 27 15:47:29 2015
+%                    (c) 2015 A.M. Thurnherr
+%                    uE-Info: 21 0 NIL 0 0 72 2 2 8 NIL ofnI
 %======================================================================
 
-% NOTES:
-%	- This version creates slightly different files than the original version
-%	  created by Visbeck/Krahmann. In the original version, the contents of the
-%	  dr structure end up as top-level variables and the contents of
-%	  the da, p, ps, f and att structures end of as global attributes. In 
-%	  the new version, the latter are saved as sub-structures, with _struct appended
-%	  to the internal names to avoid conflicts.
+% HISTORY:
+%  Jul 24, 2015: - created, losely based on code provided by D. Cardoso
+%  Jul 27, 2015: - reverted to Martin's original dimensions as requested
+%		   by E. Firing
 
-% HISTORY:
-%   Aug 15, 2013: - incorporated this code, supplied Diana Cardoso, into IX_10beta
-%		  - modified doc in header
-%		  - renamded struct variable to dr_struct
-%		  - removed 'cd' in and out of results directory (pathnames work just fine)
-%		  - delete netcdfile before it is written to (old 'clobber' option)
-%		  - removed 'l' suffix from all dims
-%		  - replaced yes/no logical vals by true/false
-%		  - renamed substructures from st2..st6 to internal names (da,p,ps,f,att)
-%   Aug 28, 2013: - incorporated bug fix provided by Diana Cardoso to prevent lat,lon,name and 
-%		    date to be stored 2cd in the nc file, which can make the code
-%		    bomb if the length of any other var is 6 (or equal to the length of name?)
+function [] = ladcp2cdf(fname,dr,da,p,ps,f,att);
+    netcdfile = deblank(fname);
+    if exist(netcdfile,'file'), delete(netcdfile); end
 
-% check arguments
-if nargin<2
-  error('need two input arguments')
-end
-if ~isstruct(dr_struct)
-  error('second argument must be a dr structure')
-end
-
-netcdfile = deblank(fname); %remove any blanks from string end
-if exist(netcdfile,'file')
-	delete(netcdfile)
+    add_dr_struct(netcdfile,'dr',dr,att)
+    add_struct(netcdfile,'da',da)
+    add_struct(netcdfile,'p',p)
+    add_struct(netcdfile,'ps',ps)
+    add_struct(netcdfile,'f',f)
 end
 
-%Determine dimensions of variables
-lbot = 0;
-lz = 0;
-ltim = 0;
-lsadcp = 0;
+%----------------------------------------------------------------------
+
+function [] = addatts(ncf,fnm,att)
+    as = getfield(att,fnm);
+    an = fieldnames(as);
+    for i=1:length(an)
+    	ncwriteatt(ncf,sprintf('dr.%s',fnm),char(an(i)),char(getfield(as,char(an(i)))));
+    end
+end   
+
+function [] = add_dr_struct(ncf,snm,dr,att)
 
-if isfield(dr_struct,'z');
-  lz = length(getfield(dr_struct,'z'));
-end  
-if isfield(dr_struct,'tim');
-  ltim = length(getfield(dr_struct,'tim'));
-end  
-if isfield(dr_struct,'zbot');
-  lbot = length(getfield(dr_struct,'zbot'));
-end  
-if isfield(dr_struct,'z_sadcp');
-  lsadcp = length(getfield(dr_struct,'z_sadcp'));
-end  
+    % scalars & misc dims
+    nccreate(ncf,'dr.name','Datatype','char','Dimensions',{'dr.name' length(dr.name)});
+	ncwrite(ncf,'dr.name',dr.name,[1]);
+	addatts(ncf,'name',att);
+    nccreate(ncf,'dr.date','Datatype','int32','Dimensions',{'dr.date' length(dr.date)});
+	ncwrite(ncf,'dr.date',dr.date,[1]);
+	addatts(ncf,'date',att);
+    nccreate(ncf,'dr.lat','Datatype','double');
+    	ncwrite(ncf,'dr.lat',dr.lat);
+    	addatts(ncf,'lat',att);
+    nccreate(ncf,'dr.lon','Datatype','double');
+    	ncwrite(ncf,'dr.lon',dr.lon);
+    	addatts(ncf,'lon',att);
+    nccreate(ncf,'dr.ubar','Datatype','double');
+    	ncwrite(ncf,'dr.ubar',dr.ubar);
+    	addatts(ncf,'ubar',att);
+    nccreate(ncf,'dr.vbar','Datatype','double');
+    	ncwrite(ncf,'dr.vbar',dr.vbar);
+    	addatts(ncf,'vbar',att);
 
-% % define dimensions in netcdf file and standard variables
-nccreate(netcdfile,'lat','Dimensions',{'lat' 1},'Datatype','single');
-nccreate(netcdfile,'lon','Dimensions',{'lon' 1},'Datatype','single');
-nccreate(netcdfile,'date','Dimensions',{'date' 6},'Datatype','int32');
-nccreate(netcdfile,'name','Dimensions',{'name' length(getfield(dr_struct,'name'))},'Datatype','char');
+    % zbot dim
+    nccreate(ncf,'dr.zbot','Datatype','double','Dimensions',{'zbot' length(dr.zbot)});
+	ncwrite(ncf,'dr.zbot',dr.zbot,[1]);
+    	addatts(ncf,'zbot',att);
+    nccreate(ncf,'dr.ubot','Datatype','double','Dimensions',{'zbot'});
+	ncwrite(ncf,'dr.ubot',dr.ubot,[1]);
+    	addatts(ncf,'ubot',att);
+    nccreate(ncf,'dr.vbot','Datatype','double','Dimensions',{'zbot'});
+	ncwrite(ncf,'dr.vbot',dr.vbot,[1]);
+    	addatts(ncf,'vbot',att);
+    nccreate(ncf,'dr.uerrbot','Datatype','double','Dimensions',{'zbot'});
+	ncwrite(ncf,'dr.uerrbot',dr.uerrbot,[1]);
+    	addatts(ncf,'uerrbot',att);
 
-% store standard variables
-ncwrite(netcdfile,'lat',dr_struct.lat);
-ncwrite(netcdfile,'lon',dr_struct.lon);
-ncwrite(netcdfile,'date',dr_struct.date);
-ncwrite(netcdfile,'name',dr_struct.name);
+    % z_sadcp dim
+    nccreate(ncf,'dr.z_sadcp','Datatype','double','Dimensions',{'z_sadcp' length(dr.z_sadcp)});
+	ncwrite(ncf,'dr.z_sadcp',dr.z_sadcp,[1]);
+    	addatts(ncf,'z_sadcp',att);
+    nccreate(ncf,'dr.u_sadcp','Datatype','double','Dimensions',{'z_sadcp'});
+	ncwrite(ncf,'dr.u_sadcp',dr.u_sadcp,[1]);
+    	addatts(ncf,'u_sadcp',att);
+    nccreate(ncf,'dr.v_sadcp','Datatype','double','Dimensions',{'z_sadcp'});
+	ncwrite(ncf,'dr.v_sadcp',dr.v_sadcp,[1]);
+    	addatts(ncf,'v_sadcp',att);
+    nccreate(ncf,'dr.uerr_sadcp','Datatype','double','Dimensions',{'z_sadcp'});
+	ncwrite(ncf,'dr.uerr_sadcp',dr.uerr_sadcp,[1]);
+    	addatts(ncf,'uerr_sadcp',att);
 
-% parse fieldnames, define the proper variable and store it
-fnames = fieldnames(dr_struct);
-nn=strncmp('name',fnames,6); nla=strncmp('lat',fnames,3); 	% find, name, lat,lon,date from fnames
-nlo=strncmp('lon',fnames,3); nda=strncmp('date',fnames,4);
-ntot=[nn+nla+nlo+nda]; Ktot = logical(ntot);
-fnames(Ktot,:)=[];						% remove , name, lat,lon,date from fnames
+    % z dim
+    nccreate(ncf,'dr.z','Datatype','double','Dimensions',{'z' length(dr.z)});
+	ncwrite(ncf,'dr.z',dr.z,[1]);
+    	addatts(ncf,'z',att);
+    nccreate(ncf,'dr.u','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.u',dr.u,[1]);
+    	addatts(ncf,'u',att);
+    nccreate(ncf,'dr.v','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.v',dr.v,[1]);
+    	addatts(ncf,'v',att);
+    nccreate(ncf,'dr.nvel','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.nvel',dr.nvel,[1]);
+    	addatts(ncf,'nvel',att);
+    nccreate(ncf,'dr.uerr','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.uerr',dr.uerr,[1]);
+    	addatts(ncf,'uerr',att);
+    nccreate(ncf,'dr.range','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.range',dr.range,[1]);
+    	addatts(ncf,'range',att);
+    nccreate(ncf,'dr.range_do','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.range_do',dr.range_do,[1]);
+    	addatts(ncf,'range_do',att);
+    nccreate(ncf,'dr.range_up','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.range_up',dr.range_up,[1]);
+    	addatts(ncf,'range_up',att);
+    nccreate(ncf,'dr.ts','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.ts',dr.ts,[1]);
+    	addatts(ncf,'ts',att);
+    nccreate(ncf,'dr.ts_out','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.ts_out',dr.ts_out,[1]);
+    	addatts(ncf,'ts_out',att);
+    nccreate(ncf,'dr.p','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.p',dr.p,[1]);
+    	addatts(ncf,'p',att);
+    nccreate(ncf,'dr.ctd_t','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.ctd_t',dr.ctd_t,[1]);
+    	addatts(ncf,'ctd_t',att);
+    nccreate(ncf,'dr.ctd_s','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.ctd_s',dr.ctd_s,[1]);
+    	addatts(ncf,'ctd_s',att);
+    nccreate(ncf,'dr.u_do','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.u_do',dr.u_do,[1]);
+    	addatts(ncf,'u_do',att);
+    nccreate(ncf,'dr.v_do','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.v_do',dr.v_do,[1]);
+    	addatts(ncf,'v_do',att);
+    nccreate(ncf,'dr.u_up','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.u_up',dr.u_up,[1]);
+    	addatts(ncf,'u_up',att);
+    nccreate(ncf,'dr.v_up','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.v_up',dr.v_up,[1]);
+    	addatts(ncf,'v_up',att);
+    nccreate(ncf,'dr.ensemble_vel_err','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.ensemble_vel_err',dr.ensemble_vel_err,[1]);
+    	addatts(ncf,'ensemble_vel_err',att);
+    nccreate(ncf,'dr.u_shear_method','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.u_shear_method',dr.u_shear_method,[1]);
+    	addatts(ncf,'u_shear_method',att);
+    nccreate(ncf,'dr.v_shear_method','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.v_shear_method',dr.v_shear_method,[1]);
+    	addatts(ncf,'v_shear_method',att);
+    nccreate(ncf,'dr.w_shear_method','Datatype','double','Dimensions',{'z'});
+	ncwrite(ncf,'dr.w_shear_method',dr.w_shear_method,[1]);
+    if existf(dr,'ctd_ss')
+	nccreate(ncf,'dr.ctd_ss','Datatype','double','Dimensions',{'z'});
+	    ncwrite(ncf,'dr.ctd_ss',dr.ctd_ss,[1]);
+	    addatts(ncf,'ctd_ss',att);
+    end
+    if existf(dr,'ctd_N2')
+	nccreate(ncf,'dr.ctd_N2','Datatype','double','Dimensions',{'z'});
+	    ncwrite(ncf,'dr.ctd_N2',dr.ctd_N2,[1]);
+	    addatts(ncf,'ctd_N2',att);
+    end
 
-for n=1:size(fnames,1)
-  dummy = getfield(dr_struct,fnames{n});
-  if length(dummy)==lz
-    nccreate(netcdfile,fnames{n},'Dimensions',{fnames{n} fix(lz)},'Datatype','single');
-    ncwrite(netcdfile,fnames{n},dummy);
-  end
-  if length(dummy)==ltim
-    nccreate(netcdfile,fnames{n},'Dimensions',{fnames{n} fix(ltim)},'Datatype','single');
-    ncwrite(netcdfile,fnames{n},dummy);
-  end
-  if length(dummy)==lbot
-    nccreate(netcdfile,fnames{n},'Dimensions',{fnames{n} fix(lbot)},'Datatype','single');
-    ncwrite(netcdfile,fnames{n},dummy);
-  end
-  if length(dummy)==lsadcp
-    nccreate(netcdfile,fnames{n},'Dimensions',{fnames{n} fix(lsadcp)},'Datatype','single');
-    ncwrite(netcdfile,fnames{n},dummy);
-  end
-end
-
-add_struct(netcdfile,'da_struct',da)
-add_struct(netcdfile,'p_struct',p)
-add_struct(netcdfile,'ps_struct',ps)
-add_struct(netcdfile,'f_struct',f)
-add_struct(netcdfile,'att_struct',att)
-
+    % tim dim
+    nccreate(ncf,'dr.tim','Datatype','double','Dimensions',{'tim' length(dr.tim)});
+	ncwrite(ncf,'dr.tim',dr.tim,[1]);
+    	addatts(ncf,'tim',att);
+    nccreate(ncf,'dr.tim_hour','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.tim_hour',dr.tim_hour,[1]);
+    	addatts(ncf,'tim_hour',att);
+    nccreate(ncf,'dr.shiplon','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.shiplon',dr.shiplon,[1]);
+    	addatts(ncf,'shiplon',att);
+    nccreate(ncf,'dr.shiplat','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.shiplat',dr.shiplat,[1]);
+    	addatts(ncf,'shiplat',att);
+    nccreate(ncf,'dr.xship','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.xship',dr.xship,[1]);
+    	addatts(ncf,'xship',att);
+    nccreate(ncf,'dr.yship','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.yship',dr.yship,[1]);
+    	addatts(ncf,'yship',att);
+    nccreate(ncf,'dr.uship','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.uship',dr.uship,[1]);
+    	addatts(ncf,'uship',att);
+    nccreate(ncf,'dr.vship','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.vship',dr.vship,[1]);
+    	addatts(ncf,'vship',att);
+    nccreate(ncf,'dr.zctd','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.zctd',dr.zctd,[1]);
+    	addatts(ncf,'zctd',att);
+    nccreate(ncf,'dr.wctd','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.wctd',dr.wctd,[1]);
+    	addatts(ncf,'wctd',att);
+    nccreate(ncf,'dr.uctd','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.uctd',dr.uctd,[1]);
+    	addatts(ncf,'uctd',att);
+    nccreate(ncf,'dr.vctd','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.vctd',dr.vctd,[1]);
+    	addatts(ncf,'vctd',att);
+    nccreate(ncf,'dr.xctd','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.xctd',dr.xctd,[1]);
+    	addatts(ncf,'xctd',att);
+    nccreate(ncf,'dr.yctd','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.yctd',dr.yctd,[1]);
+    	addatts(ncf,'yctd',att);
+    nccreate(ncf,'dr.uctderr','Datatype','double','Dimensions',{'tim'});
+	ncwrite(ncf,'dr.uctderr',dr.uctderr,[1]);
+    	addatts(ncf,'uctderr',att);
 end % function
 
 %----------------------------------------------------------------------
 
-function [] = add_struct(ncf,snm,a)
-  
-   fnames = fieldnames(a);
-   if isstruct(a)
-      if ~isstruct(eval(['a.' fnames{1}])) % No SubStructure
-	nccreate(ncf,snm,'Datatype','char');
+function [] = add_struct(ncf,snm,struct)
+
+    fname = fieldnames(struct);
+    for i=1:length(fname)
+	fns = char(fname(i));
+	vnm = sprintf('%s.%s',snm,fns);
+    
+	f = getfield(struct,fns);
+	if isstruct(f)
+	    error(sprintf('ladcp2cdf:add_struct(%s) substructures are not allowed',vnm));
+	end
+    
+	if ischar(f),		type = 'char';					% define data type
+	elseif isnumeric(f),	type = 'double';
+	elseif islogical(f),	type = 'int8'; if f, f=1; else, f=0; end;	% logical -> int
+	else, error(sprintf('ladcp2cdf:create_var(%s) unsupported type',vnm));
+	end
 
-	for n = 1:size(fnames,1)
-		dummy = getfield(a,fnames{n});
-	    	if size(dummy,1)==1
-	       		if isstr(dummy)
-                		ncwriteatt(ncf,snm,fnames{n},dummy);
-	       		elseif islogical(dummy)
-		        	if dummy, dummy='true';
-		               	else, 	  dummy='false';
-		               	end
-				ncwriteatt(ncf,snm,fnames{n},dummy);
-			else
-		                ncwriteatt(ncf,snm,fnames{n},dummy(:));
-			end    
-		end
-	end % for n
-      else % SubStructures -> Variable Attributes
-	for n = 1:size(fnames,1)
-		atts = eval(['fieldnames(a.' fnames{n} ');']);
-	        finfo = ncinfo(ncf);
-        	FieldNames = {finfo. Variables.Name};
-	        existField=strmatch(fnames{n}, FieldNames);
-        	if isempty(existField)
-			nccreate(ncf,fnames{n},'Datatype','char');
-	        end
-		for j = 1:size(atts,1)
-			dummy = eval(['a.' fnames{n} '.' atts{j} ';']);
-			if size(dummy,1) == 1
-				if ischar(dummy)
-		                    ncwriteatt(ncf,fnames{n},atts{j},dummy);
-                		else
-		                    ncwriteatt(ncf,fnames{n},atts{j},dummy(:));
-                		end
-			end
-		end	       
-	 end % for n
-      end % else (substructures or not)
-   else % if issstruct(a)
-      disp(' not structure')
-   end
+	[nr nc] = size(f);							% define variable
+	if (nr*nc < 2)								% scalar
+	    nccreate(ncf,vnm,'Datatype',type);
+	elseif (nr==1 || nc==1) 						% vector
+	    nccreate(ncf,vnm,'Datatype',type,'Dimensions',{vnm nr*nc});
+	else									% matrix
+	    nccreate(ncf,vnm,'Datatype',type,'Dimensions',{sprintf('%s_c',vnm) nc sprintf('%s_r',vnm) nr});
+	end
+
+	if (nr*nc == 1) 							% write data: scalar
+	    %disp(sprintf('Writing one %s value to %s...',type,vnm));
+	    ncwrite(ncf,vnm,f);
+	elseif (nr==1 || nc==1) 						% vector
+	    %disp(sprintf('Writing %d %s values to %s...',nr*nc,type,vnm));
+	    ncwrite(ncf,vnm,f,[1]);
+	elseif (nr*nc > 1)							% matrix
+	    %disp(sprintf('Writing %dx%d %s values to %s...',nr,nc,type,vnm));
+	    ncwrite(ncf,vnm,f',[1 1]);
+        end % if
+    end % for
 end % function
 
+
+
--- a/loadctd.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/loadctd.m	Fri Sep 04 06:37:43 2015 +0000
@@ -1,9 +1,9 @@
 %======================================================================
 %                    L O A D C T D . M 
 %                    doc: Sat Jun 26 15:56:43 2004
-%                    dlm: Fri Mar 21 11:01:39 2014
+%                    dlm: Thu May 28 07:51:39 2015
 %                    (c) 2004 M. Visbeck & A. Thurnherr
-%                    uE-Info: 52 1 NIL 0 0 72 0 2 8 NIL ofnI
+%                    uE-Info: 169 24 NIL 0 0 72 0 2 8 NIL ofnI
 %======================================================================
 
 function [d,p]=loadctd(f,d,p)
@@ -83,6 +83,8 @@
 %   Jan  7, 2009: - tightened use of exist()
 %   Jun 16, 2009: - BUG: patching short nav time series did not work correctly
 %   Mar 21, 2014: - BUG: f.ctd_time_base used p.ctd_time_base set as default
+%   May 27, 2015: - removed confusing diagnostic message regarding adjusting NAV time
+%   May 28, 2015: - added error message when there are no valid vertical velocities
 
 % read SEABIRD ctd timeseries file
 disp(['LOADCTD: load CTD time series ',f.ctd])
@@ -164,6 +166,9 @@
 
 % calc LADCP depth
 w=meannan(d.rw);
+if sum(isfinite(w)) == 0
+    error('No valid vertical velocities --- aborting');
+end
 ii=find(~isfinite(w));
 w(ii)=0;
 dt=diff(d.time_jul)*24*3600;
@@ -204,8 +209,10 @@
   timctd = timctd + p.ctdtimoff;
 end
 
-if p.navdata & f.nav_time_base == 0
-  disp(sprintf(' adjusting GPS time to CTD time (%+d seconds)',floor(p.ctdtimoff*24*3600)))
+if p.navdata && f.nav_time_base == 0
+  if ~strcmp(f.ctd,f.nav)
+    disp(sprintf(' adjusting GPS time to CTD time (%+d seconds)',floor(p.ctdtimoff*24*3600)))
+  end
   d.navtime_jul = d.navtime_jul + p.ctdtimoff;
 end
 
--- a/loadrdi.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/loadrdi.m	Fri Sep 04 06:37:43 2015 +0000
@@ -13,9 +13,9 @@
 %======================================================================
 %                    L O A D R D I . M 
 %                    doc: Fri Jun 18 18:21:56 2004
-%                    dlm: Wed Apr 15 09:53:54 2015
+%                    dlm: Wed May 27 08:17:42 2015
 %                    (c) 2004 ladcp@
-%                    uE-Info: 260 45 NIL 0 0 72 2 2 8 NIL ofnI
+%                    uE-Info: 51 41 NIL 0 0 72 2 2 8 NIL ofnI
 %======================================================================
 
 % CHANGES BY ANT
@@ -48,6 +48,7 @@
 %		 - added separate nbin, blnk, dist for DL/UL to p struct
 %  Jan 23, 2015: - made updown() bomb when UL file is not found
 %  Apr 15, 2015: - modified ambiguity-velocity warning as suggested by Diana Cardoso
+%  May 27, 2015: - clarified time-related warnings
 
 % p=setdefv(p,'pg_save',[1 2 3 4]);
 % Default =3 for loadctd_whoi.
@@ -257,7 +258,7 @@
 jj = find(j2>skipnens & j2<size(l.u,2)-skipnens);
 
 if length(jj)>10
- warn = sprintf('** found %d (%.1f%% of total) relative horizontal velocities > %g m/s',...
+ warn = sprintf('** found %d (%.1f%% of total) velocity measurements > %g m/s',...
 		length(jj),length(jj)/length(vel)*100,p.vlim);
  disp(warn);
  if length(jj)>100
@@ -888,27 +889,25 @@
   if pingdiff | timediff | pingconst
    disp(' WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING ')
    if pingconst
-    warn=(' warning ping rate not equal in down instrument');
+    warn=(' Warning: non-constant ping rate in downlooker data (staggered pinging?)');
     disp(warn)
     disp(['  min down ping rate :',num2str(-24*3600*maxnan(-diff(timd))),...
          '  max down ping rate :',num2str(24*3600*maxnan(diff(timd)))])
    end
    if pingdiff
-    warn=(' warning mean ping rate not equal between both instruments ');
+    warn=(' Warning: mean ping rates differ in downlooker/uplooker data ');
     disp(warn)
     l.warn(size(l.warn,1)+1,1:length(warn))=warn;
     disp(['  mean down ping rate :',num2str(24*3600*meannan(diff(timd))),...
          '  mean up ping rate :',num2str(24*3600*meannan(diff(timu)))])
    end
    if timediff
-    warn=(' warning instruments have no fixed ping rate ');
+    warn=(' Warning: cast duration differs in downlooker/uplooker data ');
     disp(warn)
     l.warn(size(l.warn,1)+1,1:length(warn))=warn;
     disp(['  down dt for common ping number:',num2str((timd(ii)-timd(1))*24),...
           '  up dt :',num2str((timu(ii)-timu(1))*24),' hours '])
    end
-   %timd=timd-timd(1);
-   %timu=timu-timu(1);
    iu=1:length(timd);
    ii=find(iu>length(timu));
    iu(ii)=length(timu);
--- a/loadsadcp.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/loadsadcp.m	Fri Sep 04 06:37:43 2015 +0000
@@ -1,9 +1,9 @@
 %======================================================================
 %                    L O A D S A D C P . M 
 %                    doc: Sun Jun 27 23:42:04 2004
-%                    dlm: Tue May 22 11:01:21 2012
+%                    dlm: Fri Jun 12 11:35:40 2015
 %                    (c) 2004 ladcp@
-%                    uE-Info: 109 4 NIL 0 0 72 0 2 8 NIL ofnI
+%                    uE-Info: 97 33 NIL 0 0 72 0 2 8 NIL ofnI
 %======================================================================
 
 % CHANGES BY ANT:
@@ -15,6 +15,8 @@
 %		   because these data are unlikely to be accurate enough
 %		   for the ship-drift constraint; if they are, the user
 %		   should verify and make a GPS file during pre-processing
+%  Jun 12, 2015: - added code to set p.lat and p.lon when only SADCP 
+%		   nav data are available
 
 function   [di,p]=loadsadcp(f,di,p)
 % function   [di,p]=loadsadcp(f,di,p)
@@ -89,28 +91,34 @@
   di.sadcp_lat=interp1(tim_sadcp,lat_sadcp(:),di.time_jul);
 
   % set position from SADCP nav
+  % NB: this is only used to set lat lon in the output files; magdev is NOT calculated from these
   if abs(p.lon)+abs(p.lat)==0
-   error('as of version IX_9, using GPS info from SADCP data stream is no longer supported');
-   slat=di.sadcp_lat(1);
-   slon=di.sadcp_lon(1);
-   elat=di.sadcp_lat(end);
-   elon=di.sadcp_lon(end);
-   p.poss=[fix(slat), (slat-fix(slat))*60, fix(slon), (slon-fix(slon))*60];
-   p.pose=[fix(elat), (slat-fix(elat))*60, fix(elon), (slon-fix(elon))*60];
+    p.lon = medianan(di.sadcp_lon);
+    p.lat = medianan(di.sadcp_lat);
   end
+
+%  if abs(p.lon)+abs(p.lat)==0
+%   error('as of version IX_9, using GPS info from SADCP data stream is no longer supported');
+%   slat=di.sadcp_lat(1);
+%   slon=di.sadcp_lon(1);
+%   elat=di.sadcp_lat(end);
+%   elon=di.sadcp_lon(end);
+%   p.poss=[fix(slat), (slat-fix(slat))*60, fix(slon), (slon-fix(slon))*60];
+%   p.pose=[fix(elat), (slat-fix(elat))*60, fix(elon), (slon-fix(elon))*60];
+%  end
  
-% if no other ship navigation exists, use SADCP navigation
-  if existf(di,'slon')==0
-   error('as of version IX_9, using GPS info from SADCP data stream is no longer supported');
-   di.slon=di.sadcp_lon;
-   di.slat=di.sadcp_lat;
-  else
-   if sum(isfinite(di.slon+di.slat))==0
-    error('as of version IX_9, using GPS info from SADCP data stream is no longer supported');
-    di.slon=di.sadcp_lon;
-    di.slat=di.sadcp_lat;
-   end
-  end
+%% if no other ship navigation exists, use SADCP navigation
+%  if existf(di,'slon')==0
+%   error('as of version IX_9, using GPS info from SADCP data stream is no longer supported');
+%   di.slon=di.sadcp_lon;
+%   di.slat=di.sadcp_lat;
+%  else
+%   if sum(isfinite(di.slon+di.slat))==0
+%    error('as of version IX_9, using GPS info from SADCP data stream is no longer supported');
+%    di.slon=di.sadcp_lon;
+%    di.slat=di.sadcp_lat;
+%   end
+%  end
 
   if length(find(isfinite(u_sadcp(:,ii))))<1, 
    disp(' no finite SADCP data found ')
--- a/plotraw.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/plotraw.m	Fri Sep 04 06:37:43 2015 +0000
@@ -6,15 +6,17 @@
 %======================================================================
 %                    P L O T R A W . M 
 %                    doc: Fri Jan  5 15:38:43 2007
-%                    dlm: Mon Jun  3 14:09:21 2013
+%                    dlm: Wed Apr 15 11:30:50 2015
 %                    (c) 2007 M. Visbeck with contribs from A. Thurnherr
-%                    uE-Info: 79 0 NIL 0 0 72 0 2 4 NIL ofnI
+%                    uE-Info: 19 76 NIL 0 0 72 0 2 4 NIL ofnI
 %======================================================================
 
 % MODIFICATIONS BY ANT:
 %	Jan  5, 2007: - fixed checkbeam() as suggested by B. Huber
 %	Jun  3, 2013: - BUG: top panel of Fig. 2 was wrong for dual-headed
 %						 LADCPs with different UL/DL bin sizes
+%	Apr 15, 2015: - BUG: ax() is always used but only sometimes defined()
+%						 reported by Achim Randelhoff (yahoo email 04/23/14)
 
 pmax=200;
 orient tall
@@ -74,6 +76,8 @@
 	ax(4) = -zz(end);
 	axis(ax);
 	pcolorn(ii,zz(length(d.izu)+1:end),rw(length(d.izu)+1:end,:))
+else
+	ax = axis;	% alternative solution to bug reported by Achim Randelhoff yahoo email 04/23/2014
 end
 
 hold on
--- a/process_cast.m	Wed Apr 15 09:55:34 2015 +0000
+++ b/process_cast.m	Fri Sep 04 06:37:43 2015 +0000
@@ -34,16 +34,15 @@
 % 13: (RE-)LOAD SADCP DATA
 % 14: CALCULATE INVERSE SOLUTION
 % 15: CALCULATE SHEAR SOLUTION
-% 16: CALCULATE DIFFUSIVITY PROFILE
-% 17: PLOT RESULTS & SHOW WARNINGS
-% 18: SAVE OUTPUT
+% 16: PLOT RESULTS & SHOW WARNINGS
+% 17: SAVE OUTPUT
 
 %======================================================================
 %                    P R O C E S S _ C A S T . M 
 %                    doc: Thu Jun 24 16:54:23 2004
-%                    dlm: Fri Sep 26 15:57:32 2014
+%                    dlm: Tue May 12 08:36:53 2015
 %                    (c) 2004 A.M. Thurnherr
-%                    uE-Info: 83 80 NIL 0 0 72 2 2 8 NIL ofnI
+%                    uE-Info: 378 0 NIL 0 0 72 2 2 8 NIL ofnI
 %======================================================================
 
 % NOTES:
@@ -81,6 +80,7 @@
 %			data were loaded
 %  Apr 26, 2012: - finally removed finestructure kz code
 %  Sep 26, 2014: - added support for p.orig in [saveres.m] (patch by Dan Torres)
+%  May 12, 2015: - finally removed entire step 16 (diffusivity)
 
 %----------------------------------------------------------------------
 % STEP 0: EXECUTE ALWAYS
@@ -368,6 +368,13 @@
      dino=di;
      lanarrow
      diary on
+     if existf(d,'ctdprof_p')
+       dr.ctd_t=interp1q(d.ctdprof_z,d.ctdprof_t,dr.z);
+       dr.ctd_s=interp1q(d.ctdprof_z,d.ctdprof_s,dr.z);
+     end
+     if existf(d,'ctdprof_ss')
+       dr.ctd_ss=interp1q(d.ctdprof_z,d.ctdprof_ss,dr.z);
+     end
   end
 
   end_processing_step;
@@ -459,30 +466,7 @@
 end % OF STEP 15: CALCULATE SHEAR SOLUTION
 
 %----------------------------------------------------------------------
-% STEP 16: CALCULATE DIFFUSIVITY PROFILE
-%----------------------------------------------------------------------
-
-pcs.cur_step = pcs.cur_step + 1;
-if pcs.begin_step <= pcs.cur_step
-  pcs.step_name = 'CALCULATE DIFFUSIVITY PROFILE'; begin_processing_step;
-
-  % 
-  %  prepare ctd data for output
-  %
-  if existf(d,'ctdprof_p')
-   % first save CTD data with profile
-   dr.ctd_t=interp1q(d.ctdprof_z,d.ctdprof_t,dr.z);
-   dr.ctd_s=interp1q(d.ctdprof_z,d.ctdprof_s,dr.z);
-  end
-  if existf(d,'ctdprof_ss')
-   dr.ctd_ss=interp1q(d.ctdprof_z,d.ctdprof_ss,dr.z);
-  end
-  
-  end_processing_step;
-end % OF STEP 16: CALCULATE DIFFUSIVITY PROFILE
-
-%----------------------------------------------------------------------
-% STEP 17: PLOT RESULTS & SHOW WARNINGS
+% STEP 16: PLOT RESULTS & SHOW WARNINGS
 %----------------------------------------------------------------------
 
 pcs.cur_step = pcs.cur_step + 1;
@@ -529,10 +513,10 @@
   pause(0.01)
 
   end_processing_step;
-end % OF STEP 17: PLOT RESULTS & SHOW WARNINGS
+end % OF STEP 16: PLOT RESULTS & SHOW WARNINGS
 
 %----------------------------------------------------------------------
-% STEP 18: SAVE OUTPUT
+% STEP 17: SAVE OUTPUT
 %----------------------------------------------------------------------
 
 pcs.cur_step = pcs.cur_step + 1;
@@ -584,7 +568,7 @@
   end
     
   end_processing_step;
-end % OF STEP 18: SAVE OUTPUT
+end % OF STEP 17: SAVE OUTPUT
 
 %----------------------------------------------------------------------
 % FINAL STEP: CLEAN UP