#======================================================================

#                    A N T S N C . P L 

#                    doc: Mon Jul 17 11:59:37 2006

#                    dlm: Thu Oct 13 14:27:02 2022

#                    (c) 2006 A.M. Thurnherr

#                    uE-Info: 492 55 NIL 0 0 70 0 2 4 NIL ofnI

#======================================================================

# ANTS netcdf library

# HISTORY:

#	Jul 17, 2006: - created

#	Jul 21, 2006: - documented

#				  - added NC-encoding routines

#	Jul 22: 2006: - BUG: pseudo %PARAMs were written as well

#				  -	BUG: var ATTRs were not enconded correctly

#				  - added type support

#	Jul 23, 2006: - improved type magic

#	Sep  1, 2006: - BUG: removing trainling 0s had not worked

#	Sep 23, 2006: - fiddled

#	Jul 11, 2008: - adapted to new pseudo %PARAMs

#	Jul 16, 2008: - remove \0s from strings in NC_stringify

#	Mar 20, 2008: - added progress output to NC_stringify

#	Jul 21, 2009: - allowed for suppression of %PARAMs

#	Jan 15, 2016: - BUG: %DEPS pseudo-%PARAM was encoded

#	Apr  5, 2022: - added NC_writeMDataMulti()

#				  - implemented disappeared NetCDF::DOUBLE etc.

#				  - enabled fill value handling (library bug has gone away)

#	Apr 22, 2022: - added "_coordinate" to abscissa dimension to

#					allow reading with python xr

#	Jun 18, 2022: - found better earlier solution for NetCDF::DOUBLE etc. 

#					on laptop

#	Oct 13, 2022: - added global attributes to NC_writeMDataMulti

# HISTORY END

# NOTES:

#	- multi-valued attribs are not loaded by getInfo()

#	- spaces in NC strings are replaced by underscores

# NetCDF Library Bug: NO LONGER ACTIVE AS OF APR 2022

#	The library appears to have incorrect default _FillValue types for

#	integer data types. The error appears if the "setfill" line is commented

#	out and the following command is run:

#		listNC -ct dbk100.nc | NCode -o TEMP.nc time

#	NB: The error occurs when the 1st variable value is written, NOT when

#	    the first Q_time value is written. However, when all the Q_ fields

#		are ommitted, the error disappears.

use NetCDF;

#----------------------------------------------------------------------

# NetCDF Constants

#	- in April 2022 I encoded these constants from .h file

#	- in June 2022 I found a more general earlier solution on my laptop

#----------------------------------------------------------------------

sub NetCDF::NAT 		 { return NetCDF::constant(NAT,0); }

sub NetCDF::BYTE		 { return NetCDF::constant(BYTE,0); }

sub NetCDF::CHAR		 { return NetCDF::constant(CHAR,0); }

sub NetCDF::SHORT		 { return NetCDF::constant(SHORT,0); }

sub NetCDF::INT 		 { return NetCDF::constant(INT,0); }

sub NetCDF::LONG		 { return NetCDF::constant(LONG,0); }

sub NetCDF::FLOAT		 { return NetCDF::constant(FLOAT,0); }

sub NetCDF::DOUBLE		 { return NetCDF::constant(DOUBLE,0); }

sub NetCDF::UBYTE		 { return NetCDF::constant(UBYTE,0); }

sub NetCDF::USHORT		 { return NetCDF::constant(USHORT,0); }

sub NetCDF::UINT		 { return NetCDF::constant(UINT,0); }

sub NetCDF::INT64		 { return NetCDF::constant(INT64,0); }

sub NetCDF::UINT64		 { return NetCDF::constant(UINT64,0); }

sub NetCDF::STRING		 { return NetCDF::constant(STRING,0); }

sub NetCDF::FILL_BYTE		 { return NetCDF::constant(FILL_BYTE,0); }

sub NetCDF::FILL_CHAR		 { return NetCDF::constant(FILL_CHAR,0); }

sub NetCDF::FILL_SHORT		 { return NetCDF::constant(FILL_SHORT,0); }

sub NetCDF::FILL_INT		 { return NetCDF::constant(FILL_INT,0); }

sub NetCDF::FILL_LONG		 { return NetCDF::constant(FILL_LONG,0); }

sub NetCDF::FILL_FLOAT		 { return NetCDF::constant(FILL_FLOAT,0)+36; }

sub NetCDF::FILL_DOUBLE 	 { return NetCDF::constant(FILL_DOUBLE,0)+36; }

sub NetCDF::FILL_UBYTE		 { return NetCDF::constant(FILL_UBYTE,0); }

sub NetCDF::FILL_USHORT 	 { return NetCDF::constant(FILL_USHORT,0); }

sub NetCDF::FILL_UINT		 { return NetCDF::constant(FILL_UINT,0); }

sub NetCDF::FILL_INT64		 { return NetCDF::constant(FILL_INT64,0); }

sub NetCDF::FILL_UINT64 	 { return NetCDF::constant(FILL_UINT64,0); }

sub NetCDF::FILL_STRING 	 { return NetCDF::constant(FILL_STRING,0); }

sub NetCDF::NOWRITE 	 { return NetCDF::constant(NOWRITE,0); }

sub NetCDF::CLOBBER 	 { return NetCDF::constant(CLOBBER,0); }

sub NetCDF::NOFILL		 { return NetCDF::constant(NOFILL,0); }

sub NetCDF::GLOBAL		 { return NetCDF::constant(GLOBAL,0); }

sub NetCDF::UNLIMITED	 { return NetCDF::constant(UNLIMITED,0); }

if (0) {

	sub NetCDF::NAT 		 { return 0; }

	sub NetCDF::BYTE		 { return 1; }

	sub NetCDF::CHAR		 { return 2; }

	sub NetCDF::SHORT		 { return 3; }

	sub NetCDF::INT 		 { return 4; }

	sub NetCDF::LONG		 { return 4; }

	sub NetCDF::FLOAT		 { return 5; }

	sub NetCDF::DOUBLE		 { return 6; }

	sub NetCDF::UBYTE		 { return 7; }

	sub NetCDF::USHORT		 { return 8; }

	sub NetCDF::UINT		 { return 9; }

	sub NetCDF::INT64		 { return 10; }

	sub NetCDF::UINT64		 { return 11; }

	sub NetCDF::STRING		 { return 12; }

	sub NetCDF::FILL_BYTE		 { return -127; }

	sub NetCDF::FILL_CHAR		 { return "\0"; }

	sub NetCDF::FILL_SHORT		 { return -32767; }

	sub NetCDF::FILL_INT		 { return -2147483647; }

	sub NetCDF::FILL_LONG		 { return -2147483647; }

	sub NetCDF::FILL_FLOAT		 { return 9.9692099683868690e+36; }

	sub NetCDF::FILL_DOUBLE 	 { return 9.9692099683868690e+36; }

	sub NetCDF::FILL_UBYTE		 { return 255; }

	sub NetCDF::FILL_USHORT 	 { return 65535; }

	sub NetCDF::FILL_UINT		 { return 4294967295; }

	sub NetCDF::FILL_INT64		 { return -9223372036854775806; }

	sub NetCDF::FILL_UINT64 	 { return 18446744073709551614; }

	sub NetCDF::FILL_STRING 	 { return ''; }

	sub NetCDF::NOWRITE 	 { return 0x0000; }

	sub NetCDF::CLOBBER 	 { return 0x0000; }

	sub NetCDF::NOFILL		 { return 0x100; }

	sub NetCDF::GLOBAL		 { return -1; }

	sub NetCDF::UNLIMITED    { return 0; }

}

#----------------------------------

# string representation of NC types

#----------------------------------

sub NC_typeName($)

{

	my($tp) = @_;

	return 'byte'	if ($tp == NetCDF::BYTE);

	return 'char'	if ($tp == NetCDF::CHAR);

	return 'short'	if ($tp == NetCDF::SHORT);

	return 'long'	if ($tp == NetCDF::LONG);

	return 'float'	if ($tp == NetCDF::FLOAT);

	return 'double' if ($tp == NetCDF::DOUBLE);

	croak("$0: unknown NetCDF type #$tp\n");

}

sub NC_type($)

{

	my($tn) = lc($_[0]);

	return	NetCDF::BYTE	if ($tn eq 'byte');

	return	NetCDF::CHAR	if ($tn eq 'char');

	return	NetCDF::SHORT	if ($tn eq 'short');

	return	NetCDF::LONG	if ($tn eq 'long');

	return	NetCDF::FLOAT	if ($tn eq 'float');

	return	NetCDF::DOUBLE	if ($tn eq 'double');

	croak("$0: unknown NetCDF type <$tn>\n");

}

#--------------------------------------

# test whether given NC type is numeric

#--------------------------------------

sub NC_isNumeric($)

{

	my($tp) = @_;

	return 1 if ($tp == NetCDF::BYTE);

	return 1 if ($tp == NetCDF::SHORT);

	return 1 if ($tp == NetCDF::LONG);

	return 1 if ($tp == NetCDF::FLOAT);

	return 1 if ($tp == NetCDF::DOUBLE);

	return 0;

}

#----------------------------------------

# test whether given NC type is character

#----------------------------------------

sub NC_isChar($)

{

	return $_[0] == NetCDF::CHAR;

}

#-----------------------------------

# convert character- to string array

#-----------------------------------

sub NC_stringify($@)

{

	my($len,@chars) = @_;

	my(@strings);

	my($nStrings) = @chars/$len;

	print(STDERR "$0: extracting $nStrings strings")

		if ($nStrings > 1000);

	while (@chars) {

		print(STDERR ".") if ($nStrings>1000 && $n++%1000 == 0);

		push(@strings,pack("c$len",@chars));

		$strings[$#strings] =~ s/ /_/g;

		$strings[$#strings] =~ s/\0//g;

		splice(@chars,0,$len);

	}

	print(STDERR "\n") if ($nStrings > 1000);

	return @strings;

}

#----------------------------------------------------------------------

# open netcdf file and read (most) metadata into hash

#

#	INPUT:

#		<filename>

#

#	OUTPUT:

#		$NC{id}								netcdf id

#

#		@NC{attrName}[]						names of global attrs

#		%NC{AttrType}{$aName}				types of global attrs

#		%NC{AttrLen}{$aName}				# of elts in global attrs

#		%NC{Attr}{$aName}					vals of scalar global attrs

#

#		$NC{unlim_dimId}					dim id of unlimited dim

#		@NC{dimName}[$dimId]				dim names

#		%NC{dimID}{$dName}					dim ids

#		%NC{dimLen}{$dName}					# elts in dim

#

#		@NC{varName}[$varId]				var names

#		%NC{varType}{$vName}				var types

#		%NC{varId}{$vName}					var ids

#		@%NC{varDimIDs}{$vName}[]			dims of vars, e.g. u(lon,lat)

#		@%NC{varAttrName}{$vName}[]			names of var attrs

#		%%NC{varAttrType}{$vName}{$aName}	types of var attrs

#		%%NC{varAttrLen}{$vName}{$aName}	# of elts in var attrs

#		%%NC{varAttr}{$vName}{$aName}		vals of scalar var attrs

#

#----------------------------------------------------------------------

sub NC_readMData($)

{

	my($fn) = @_;

	my(%NC);

	$NC{id} = NetCDF::ncopen($ARGV[0],NetCDF::NOWRITE);	# open

	my($nd,$nv,$nga,$udi);								# get nelts

	NetCDF::ncinquire($NC{id},$nd,$nv,$nga,$udi);

	$NC{unlim_dimId} = $udi;

	for (my($d)=0; $d<$nd; $d++) {						# dimensions

		my($dnm,$ln);

		NetCDF::ncdiminq($NC{id},$d,$dnm,$ln);

		$NC{dimName}[$d] = $dnm;

		$NC{dimId}{$dnm} = $d;

		$NC{dimLen}{$dnm} = $ln;

	}

	for (my($v)=0; $v<$nv; $v++) {						# vars & var-attribs

		my($vnm,$vtp,$nvd,$nva);

		my(@dids) = ();

		NetCDF::ncvarinq($NC{id},$v,$vnm,$vtp,$nvd,\@dids,$nva);

		$NC{varName}[$v] = $vnm;

		$NC{varId}{$vnm} = $v;

		$NC{varType}{$vnm} = $vtp;

		@{$NC{varDimIds}{$vnm}} = @dids[0..$nvd-1];

		for (my($a)=0; $a<$nva; $a++) {					# var-attribs

			my($anm,$atp,$aln);

			NetCDF::ncattname($NC{id},$v,$a,$anm);

			$NC{varAttrName}{$vnm}[$a] = $anm;

			NetCDF::ncattinq($NC{id},$v,$anm,$atp,$aln);

			$NC{varAttrType}{$vnm}{$anm} = $atp;

			$NC{varAttrLen}{$vnm}{$anm} = $aln;

			if ($atp == NetCDF::BYTE || $atp == NetCDF::CHAR || $aln == 1) {

				my($val) = "";

				NetCDF::ncattget($NC{id},$v,$anm,\$val);

				$val =~ s{\0+$}{} if ($atp == NetCDF::CHAR);	# trailing \0

				$NC{varAttr}{$vnm}{$anm} = $val;

			}		

		}

	}

	for (my($a)=0; $a<$nga; $a++) {						#  global attribs

		my($anm,$atp,$aln);

		NetCDF::ncattname($NC{id},NetCDF::GLOBAL,$a,$anm);

		$NC{attrName}[$a] = $anm;

		NetCDF::ncattinq($NC{id},NetCDF::GLOBAL,$anm,$atp,$aln);

		$NC{attrType}{$anm} = $atp;

		$NC{attrLen}{$anm} = $aln;

		if ($atp == NetCDF::BYTE || $atp == NetCDF::CHAR || $aln == 1) {

			my($val) = "";

			NetCDF::ncattget($NC{id},NetCDF::GLOBAL,$anm,\$val);

			$val =~ s{\0+$}{} if ($atp == NetCDF::CHAR);

			$NC{attr}{$anm} = $val;

		}	    

    }

	return %NC;

}

#----------------------------------------------------------------------

# create new nc file and write metadata

#

#	INPUT:

#		<filename>

#		<abscissa>			name of unlimited dimension

#		<suppress-params>	if true, don't write %PARAMs

#

#	OUTPUT:

#		<netcdf id>

#

#	NOTES:

#		- netcdf types can be set with %<var>:NC_type to

#			byte, long, short, double

#		- string types are as in old PASCAL convention (e.g. string80)

#		- default type is NetCDF::DOUBLE

#		- %<var>:NC_type are not added to ATTRIBs

#----------------------------------------------------------------------

sub NC_writeMData($$$)

{

	my($fn,$abscissa,$suppress_params) = @_;

	my(%attrDone,@slDim,@NCtype);

	my($ncId) = NetCDF::nccreate($fn,NetCDF::CLOBBER);

#	NetCDF::ncsetfill($ncId,NetCDF::NOFILL);			# NetCDF library bug

														# DIMENSIONS

	my($aid) = NetCDF::ncdimdef($ncId,$abscissa . '_coordinate',NetCDF::UNLIMITED);

	for (my($f)=0; $f<=$#antsLayout; $f++) {			# types

		my($tpa) = $antsLayout[$f] . ':NC_type';

		my($sl) = ($P{$tpa} =~ m{^string(\d+)$});

		if ($sl > 0) {									# string

			$slDim[$f] = NetCDF::ncdimdef($ncId,"$antsLayout[$f]:strlen",$sl);

			$NCtype[$f] = NetCDF::CHAR;

		} elsif (defined($P{$tpa})) {					# custom

			$NCtype[$f] = NC_type($P{$tpa});

		} else {										# default

			$NCtype[$f] = NetCDF::DOUBLE;

		}

#		printf(STDERR "type %s set to %s\n",$antsLayout[$f],NC_typeName($NCtype[$f]));

		undef($P{$tpa});								# do not add to ATTRIBs

    }

	for (my($f)=0; $f<=$#antsLayout; $f++) {			# VARIABLES

		my($vid);

		if (defined($slDim[$f])) {

			$vid = NetCDF::ncvardef($ncId,$antsLayout[$f],$NCtype[$f],[$aid,$slDim[$f]]);

		} else {

			$vid = NetCDF::ncvardef($ncId,$antsLayout[$f],$NCtype[$f],[$aid]);

		}

		croak("$0: varid != fnr (implementation restriction)")

			unless ($vid == $f);

		foreach my $anm (keys(%P)) {					# VARIABLE ATTRIBUTES

			next unless defined($P{$anm});

			my($var,$attr) = ($anm =~ m{([^:]+):(.*)});

			next unless ($var eq $antsLayout[$f]);

			$attrDone{$anm} = 1;						# mark

			if (numberp($P{$anm}) || lc($P{$anm}) eq nan) {

				NetCDF::ncattput($ncId,$f,$attr,NetCDF::DOUBLE,$P{$anm});

			} else {

				NetCDF::ncattput($ncId,$f,$attr,NetCDF::CHAR,$P{$anm});

			}

        }		                  

	}

	unless ($suppress_params) {

		foreach my $anm (keys(%P)) {					# GLOBAL ATTRIBUTES

			next unless defined($P{$anm});

			next if ($anm eq 'FILENAME' || $anm eq 'DIRNAME' || # skip pseudo 

					 $anm eq 'BASENAME' || $anm eq 'EXTN' ||

					 $anm eq 'PATHNAME' || $anm eq 'DEPS' ||

					 $anm eq 'RECNO'	|| $anm eq 'LINENO');

			next if $attrDone{$anm};

			if (numberp($P{$anm}) || lc($P{$anm}) eq nan) {

				NetCDF::ncattput($ncId,NetCDF::GLOBAL,$anm,NetCDF::DOUBLE,$P{$anm});

			} else {

				NetCDF::ncattput($ncId,NetCDF::GLOBAL,$anm,NetCDF::CHAR,$P{$anm});

			}

	    }

	}

	NetCDF::ncendef($ncId);

	return $ncId;

}

#----------------------------------------------------------------------

# create new NC "multi" file and write metadata

#

#	INPUT:

#		<filename>

#		<abscissa>			name of unlimited dimension

#		<suppress-params>	flag to suppress params

#		<file-id>			name of dimension enumerating input files (e.g. profile_number, cruise_id, etc.)

#		<n-files>			number of 1-D files to encode in multi file

#		[global-attribs]	optional reference to hash with global attributes

#

#	OUTPUT:

#		<netcdf id>

#		%nc_vid				netcdf variable ids for %PARAMs

#		%nc_varAttr			flags for variable (not global) attributes

#

#	NOTES:

#		- netcdf types can be set with %<var|param>:NC_type to

#			byte, chat, long, short, double, float (NC_type())

#			- don't use same-named var and param with different types...

#		- string types are as in old PASCAL convention (e.g. string80)

#		- default type is NetCDF::DOUBLE; other types must be set

#		  explicitly

#		- %<var>:NC_type are not added to ATTRIBs

#----------------------------------------------------------------------

sub NC_writeMDataMulti($$$$$@)

{

	my($fn,$abscissa,$suppress_params,$file_id,$n_files,$global_attribs) = @_;

	my(@slDim,@NCtype);

	my($ncId) = NetCDF::nccreate($fn,NetCDF::CLOBBER);

																				# DIMENSIONS

	my($aid) = NetCDF::ncdimdef($ncId,$abscissa . '_coordinate',NetCDF::UNLIMITED);

	my($mid) = NetCDF::ncdimdef($ncId,$file_id,$n_files);

	for (my($f)=0; $f<=$#antsLayout; $f++) {									# TYPES

		my($tpa) = $antsLayout[$f] . ':NC_type';								# defined explicitly

		my($sl) = ($P{$tpa} =~ m{^string(\d+)$});

		if ($sl > 0) {															# string

			$slDim[$f] = NetCDF::ncdimdef($ncId,"$antsLayout[$f]:strlen",$sl);

			$NCtype[$f] = NetCDF::CHAR;

		} elsif (defined($P{$tpa})) {											# other custom type 

			$NCtype[$f] = NC_type($P{$tpa});

		} else {																# default type

			$NCtype[$f] = NetCDF::DOUBLE;

		}

#		printf(STDERR "type %s set to %s\n",$antsLayout[$f],NC_typeName($NCtype[$f]));

		undef($P{$tpa});														# do not add to ATTRIBs

    }

	for (my($f)=0; $f<=$#antsLayout; $f++) {									# VARIABLES

		my($vid);

		if (defined($slDim[$f])) {

			$vid = NetCDF::ncvardef($ncId,$antsLayout[$f],$NCtype[$f],[$aid,$mid,$slDim[$f]]);

		} else {

			$vid = NetCDF::ncvardef($ncId,$antsLayout[$f],$NCtype[$f],[$aid,$mid]);

		}

		croak("$0: varid != fnr (implementation restriction)")

			unless ($vid == $f);

		foreach my $anm (keys(%P)) {											# VARIABLE ATTRIBUTES

			next unless defined($P{$anm});

			my($var,$attr) = ($anm =~ m{([^:]+):(.*)});

			next unless ($var eq $antsLayout[$f]);

			$nc_varAttr{$anm} = 1;												# mark as done

			if (numberp($P{$anm}) || lc($P{$anm}) eq nan) {

				NetCDF::ncattput($ncId,$f,$attr,NetCDF::DOUBLE,$P{$anm});

			} else {

				NetCDF::ncattput($ncId,$f,$attr,NetCDF::CHAR,$P{$anm});

			}

        }		                  

	}

	unless ($suppress_params) {													# PARAM vectors

		foreach my $anm (keys(%P)) {

			next if ($anm eq 'FILENAME' || $anm eq 'DIRNAME' || 				# skip pseudo %PARAMs

					 $anm eq 'BASENAME' || $anm eq 'EXTN' ||

					 $anm eq 'PATHNAME' || $anm eq 'DEPS' ||

					 $anm eq 'RECNO'	|| $anm eq 'LINENO');

			next if $nc_varAttr{$anm};											# variable attribute (done already)

			next if ($anm =~ /:NC_type$/);										# this is needed to remove the :NC_type params!?

			my($tpa) = $anm . ':NC_type';										# defined explicitly

			my($sl) = ($P{$tpa} =~ m{^string(\d+)$});

			if ($sl > 0) {														# string

				$slDim{$anm} = NetCDF::ncdimdef($ncId,"$anm:strlen",$sl);

				$NCtype{$anm} = NetCDF::CHAR;

			} elsif (defined($P{$tpa})) {										# other custom type 

				$NCtype{$anm} = NC_type($P{$tpa});

			} else {

				$NCtype{$anm} = NetCDF::DOUBLE;										# not a NC type %PARAM

			}

			if (defined($slDim{$anm})) {										# string type

				$nc_vid{$anm} = NetCDF::ncvardef($ncId,$anm,$NCtype{$anm},[$mid,$slDim{$anm}]);

			} else {															# other type

				$nc_vid{$anm} = NetCDF::ncvardef($ncId,$anm,$NCtype{$anm},[$mid]);

			}

	    }

	}

	if (defined($global_attribs)) {												# global attributes

		foreach my $anm (keys(%$global_attribs)) {

			croak("global attributes must be numeric (implementation restriction)")

				unless numberp($$global_attribs{$anm});

			NetCDF::ncattput($ncId,NetCDF::GLOBAL,$anm,NetCDF::DOUBLE,$$global_attribs{$anm});

	    }

	}

	NetCDF::ncendef($ncId);

	return $ncId;

}

1;