rename from libGM.pl
rename to libGM76.pl
--- a/libGM.pl
+++ b/libGM76.pl
@@ -1,9 +1,9 @@
#======================================================================
-# L I B G M . P L
+# L I B G M 7 6 . P L
# doc: Sun Feb 20 14:43:47 2011
-# dlm: Tue Nov 18 12:42:30 2014
+# dlm: Sat Feb 2 02:33:11 2019
# (c) 2011 A.M. Thurnherr
-# uE-Info: 22 53 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 144 43 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -20,6 +20,14 @@
# Dec 28, 2012: - added allowance for small roundoff error to Sw()
# Oct 6, 2014: - made omega optional in Sw()
# Nov 18, 2014: - made b & jstar mandatory for Sw()
+# Feb 2, 2019: - renamed from libGM.pl to libGM76.pl
+# - replaced beta with m
+# - replaced old code based on Gregg + Kunze formalism with
+# expressions from Thurnherr et al. (GRL 2015)
+# - added GM_strain
+# - BUG: Sw usage message had wrong parameter order
+# - renamed Sw => GM_VKE; Su_z => GM_shear
+
require "$ANTS/libEOS83.pl";
@@ -38,7 +46,7 @@
# Vertical velocity spectral density
#
# Units: K.E. per frequency per wavenumber [m^2/s^2*1/s*1/m = m^3/s]
-# Version: GM79?
+# Version: GM76
#
# E. Kunze (email, Feb 2011): The GM vertical velocity w spectrum is described by
#
@@ -89,11 +97,10 @@
return 2 / $pi * $f / $omega / sqrt($omega**2 - $f**2);
}
-
-sub Sw(@)
+sub GM_VKE(@)
{
my($omega,$m,$lat,$b,$jstar,$N) =
- &antsFunUsage(-5,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]> <b[m]> <j*>',@_);
+ &antsFunUsage(-5,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
if (defined($N)) { # Sw(omega,m)
local($f) = abs(&f($lat));
@@ -102,7 +109,7 @@
return nan if ($omega < $f || $omega > $N);
my($mstar) = &m($jstar,$N,$omega);
return $GM_E0 * $b * 2 * $f**2/$omega**2/B($omega) * $jstar / ($m+$mstar)**2;
- } else { # Sw(m)
+ } else { # Sw(m), i.e. integrated over all omega; as in Thurnherr et al., GRL 2015
$N = $lat; # shift arguments to account for missing omega
$lat = $m;
local($f) = abs(&f($lat));
@@ -114,22 +121,46 @@
}
#----------------------------------------------------------------------
-# GM76, as per Gregg and Kunze (JGR 1991)
-# - beta is vertical wavenumber (m above)
+# Shear and Strain m-spectra (i.e. integrated over f)
+# - spectral density
+# - from Thurnherr et al. (GRL 2015)
#----------------------------------------------------------------------
-sub Su($$)
+sub GM_shear(@)
{
- my($beta,$N) = @_;
+ my($m,$lat,$b,$jstar,$N) =
+ &antsFunUsage(5,'fffff','<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+ local($f) = abs(&f($lat));
+ my($mstar) = &m($jstar,$N);
+ return 3 * $pi/2 * $GM_E0 * $b * $jstar * $m**2 / ($m+$mstar)**2;
+}
- my($beta_star) = &m($GM_jstar,$N); # A3
- return 3*$GM_E0*$GM_b**3*$GM_N0**2 / (2*$GM_jstar*$pi) / (1+$beta/$beta_star)**2; # A2
+sub GM_strain(@)
+{
+ my($m,$lat,$b,$jstar,$N) =
+ &antsFunUsage(5,'fffff','<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+ local($f) = abs(&f($lat));
+ my($mstar) = &m($jstar,$N);
+ return $pi/2 * $GM_E0 * $b * $jstar * $m**2 / ($m+$mstar)**2;
}
-sub Su_z($$)
-{
- my($beta,$N) = &antsFunUsage(2,'ff','<vertical wavenumber[rad/m]> <N[rad/s]>',@_);
- return $beta**2 * &Su($beta,$N);
-}
+##----------------------------------------------------------------------
+## GM76, as per Gregg and Kunze (JGR 1991)
+## - beta is vertical wavenumber
+##----------------------------------------------------------------------
+#
+#sub Su($$)
+#{
+# my($beta,$N) = @_;
+#
+# my($beta_star) = &m($GM_jstar,$N); # A3
+# return 3*$GM_E0*$GM_b**3*$GM_N0**2 / (2*$GM_jstar*$pi) / (1+$beta/$beta_star)**2; # A2
+#}
+#
+#sub Su_z($$)
+#{
+# my($beta,$N) = &antsFunUsage(2,'ff','<vertical wavenumber[rad/m]> <N[rad/s]>',@_);
+# return $beta**2 * &Su($beta,$N);
+#}
1;
new file mode 100644
--- /dev/null
+++ b/libTEOS10.pl
@@ -0,0 +1,367 @@
+#======================================================================
+# L I B T E O S 1 0 . P L
+# doc: Fri Jan 18 11:09:22 2019
+# dlm: Fri Jan 18 23:23:27 2019
+# (c) 2019 A.M. Thurnherr
+# uE-Info: 144 30 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# perl wrapper for GSW TEOS-10 libary
+
+# Notes:
+# - requires GSW shared library (currently GSW-C-3.05.0-4)
+# - requires home-grown perl-interface, implementing GSW::
+# - small subset of functions implemented
+# - no temperature scale assumed; set PARAM ITS=90|68
+# - T90*1.00024=T68
+# - check values calculated with T68
+
+# HISTORY:
+# Jan 18, 2019: - created
+
+require "$ANTS/libvec.pl";
+
+use strict;
+use GSW;
+
+#----------------------------------------------------------------------
+# ITS-68/-90 conversion
+# - call once return conversion factor
+# - divide input temperatures by conversion factor to turn into ITS-90
+# - multiply output temperatures by conv factor to return to original ITS
+# - optimized for repeat calls
+#----------------------------------------------------------------------
+
+{ # BEGIN STATIC SCOPE
+ my($TCONV);
+ sub TCONV90()
+ {
+ return $TCONV if defined($TCONV);
+ my($ITS) = &antsRequireParam('ITS');
+ if ($ITS == 68) { $TCONV = 1.00024;
+ } elsif ($ITS == 90) { $TCONV = 1;
+ } else {
+ croak("$0: illegal PARAM-value ITS=$ITS\n");
+ }
+ }
+} # END STATIC SCOPE
+
+#----------------------------------------------------------------------
+# practical salinity
+# - input: cond temp press
+# - use units from %cond.unit
+# - IDENTICAL TO EOS83 VERIFIED
+#----------------------------------------------------------------------
+
+{ my(@fc,$cscale);
+ sub salin(@)
+ {
+ my($C,$T,$P) = &antsFunUsage(3,'...','[cond, temp, press]',
+ \@fc,'cond','temp','press',@_);
+ unless (defined($cscale)) {
+ my($cu) = &antsRequireParam('cond.unit');
+ if ($cu eq 'S/m') { $cscale = 10; }
+ elsif ($cu eq 'mS/cm') { $cscale = 1; }
+ else { croak("$0: illegal PARAM-value cond.unit=$cu\n"); }
+ }
+ return numbersp($C,$T,$P) ?
+ GSW::gsw_sp_from_c($C*$cscale,$T,$P) :
+ 'nan';
+ }
+}
+
+#----------------------------------------------------------------------
+# conductivity
+# - input: salin temp press
+# - take units from %cond.unit
+# - COND(SALIN(COND)) = COND VERIFIED
+#----------------------------------------------------------------------
+
+{ my(@fc,$cscale);
+ sub cond(@)
+ {
+ my($S,$T,$P) = &antsFunUsage(3,'...','[salin, temp, press]',
+ \@fc,'salin','temp','press',@_);
+ unless (defined($cscale)) {
+ my($cu) = &antsRequireParam('cond.unit');
+ if ($cu eq 'S/m') { $cscale = 10; }
+ elsif ($cu eq 'mS/cm') { $cscale = 1; }
+ else { croak("$0: illegal PARAM-value cond.unit=$cu\n"); }
+ }
+ return numbersp($S,$T,$P) ?
+ GSW::gsw_c_from_sp($S,$T,$P)/$cscale :
+ 'nan';
+ }
+}
+
+#----------------------------------------------------------------------
+# absolute salinity
+# - input salin, press, lon, lat
+# - if %lat/%lon are available, they are used
+# - otherwise, they must be supplied as arguments
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub asalin(@)
+ {
+ my($LON,$LAT,$SP,$P);
+ $LON = &antsParam('lon');
+ $LAT = &antsParam('lat');
+ if (numbersp($LON,$LAT)) {
+ ($SP,$P) = &antsFunUsage(2,'..','[salin, press]',
+ \@fc,'salin','press',@_);
+ } else {
+ ($SP,$P,$LON,$LAT) =
+ &antsFunUsage(4,'....','[salin, press, lon, lat]',
+ \@fc,'salin','press','lon','lat',@_);
+ }
+ return numbersp($LON,$LAT,$SP,$P) ?
+ GSW::gsw_sa_from_sp($SP,$P,$LON,$LAT) :
+ 'nan';
+ }
+}
+
+#----------------------------------------------------------------------
+# potential temperature
+# - DIFFERENCES < 0.0005 degC IN FULL DEPTH PROFILE VIZ [libEOS83.pl]
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub theta(@)
+ {
+ my($SA,$T,$P,$Pref) =
+ &antsFunUsage(4,'....','[asalin, temp, press,] refpress',
+ \@fc,'asalin','temp','press',undef,@_);
+ return 'nan' unless numbersp($SA,$T,$P,$Pref);
+ return $T if ($P == $Pref);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_pt_from_t($SA,$T/$TCONV,$P,$Pref)*$TCONV;
+ }
+}
+
+#----------------------------------------------------------------------
+# conservative temperature
+# - input: asalin temp press
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub ctemp(@)
+ {
+ my($SA,$T,$P) =
+ &antsFunUsage(3,'...','[asalin, temp, press]',
+ \@fc,'asalin','temp','press',@_);
+ return 'nan' unless numbersp($SA,$T,$P);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_ct_from_t($SA,$T/$TCONV,$P)*$TCONV;
+ }
+}
+
+#----------------------------------------------------------------------
+# depth
+# - input: press lat
+# - if %lat is available, it is used
+# - CHECKED AGAINS [libEOS83.pl] => DIFF < 0.12m FOR FULL DEPTH PROF
+# - PRESS(DEPTH(PRESS)) = PRESS VERIFIED
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub depth(@)
+ {
+ my($P,$LAT);
+ $LAT = &antsParam('lat');
+ if (numberp($LAT)) {
+ ($P) = &antsFunUsage(1,'.','[press]',\@fc,'press',@_);
+ } else {
+ ($P,$LAT) = &antsFunUsage(2,'..','[press, lat]',\@fc,'press','lat',@_);
+ }
+ return numbersp($P,$LAT) ? -1*GSW::gsw_z_from_p($P,$LAT) : 'nan';
+ }
+}
+
+#----------------------------------------------------------------------
+# press
+# - input: depth lat
+# - if %lat is available, it is used
+# - PRESS(DEPTH(PRESS)) = PRESS VERIFIED
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub press(@)
+ {
+ my($D,$LAT);
+ $LAT = &antsParam('lat');
+ if (numberp($LAT)) {
+ ($D) = &antsFunUsage(1,'.','[depth]',\@fc,'depth',@_);
+ } else {
+ ($D,$LAT) = &antsFunUsage(2,'..','[depth, lat]',\@fc,'depth','lat',@_);
+ }
+ return numbersp($D,$LAT) ? GSW::gsw_p_from_z(-1*$D,$LAT) : 'nan';
+ }
+}
+
+#----------------------------------------------------------------------
+# thermal expansion coefficient
+# - input: asalin ctemp press
+# - ~ 1% higher than homegrown approximation in [libEOS83.pl]
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub alpha(@)
+ {
+ my($SA,$TC,$P) =
+ &antsFunUsage(3,'...','[asalin, ctemp, press]',
+ \@fc,'asalin','ctemp','press',@_);
+ return 'nan' unless numbersp($SA,$TC,$P);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_alpha($SA,$TC/$TCONV,$P);
+ }
+}
+
+#----------------------------------------------------------------------
+# haline contraction coefficient
+# - input: asalin ctemp press
+# - ~ 0.6% lower than homegrown approximation in [libEOS83.pl]
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub beta(@)
+ {
+ my($SA,$TC,$P) =
+ &antsFunUsage(3,'...','[asalin, ctemp, press]',
+ \@fc,'asalin','ctemp','press',@_);
+ return 'nan' unless numbersp($SA,$TC,$P);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_beta($SA,$TC/$TCONV,$P);
+ }
+}
+
+#----------------------------------------------------------------------
+# in situ density
+# - input: asalin ctemp press
+# - CHECKED RHO & ALL SIGMA AGAINST [libEOS83.pl] => CONSISTENT PATTERN
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub rho(@)
+ {
+ my($SA,$TC,$P) =
+ &antsFunUsage(3,'...','[asalin, ctemp, press]',
+ \@fc,'asalin','ctemp','press',@_);
+ return 'nan' unless numbersp($SA,$TC,$P);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_rho($SA,$TC/$TCONV,$P);
+ }
+}
+
+#----------------------------------------------------------------------
+# potential density
+# - input: asalin ctemp
+# - CHECKED RHO & ALL SIGMA AGAINST [libEOS83.pl] => CONSISTENT PATTERN
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub sigma0(@)
+ {
+ my($SA,$TC) = &antsFunUsage(2,'..','[asalin, ctemp]',\@fc,'asalin','ctemp',@_);
+ return 'nan' unless numbersp($SA,$TC);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_sigma0($SA,$TC/$TCONV);
+ }
+}
+
+{ my(@fc);
+ sub sigma1(@)
+ {
+ my($SA,$TC) = &antsFunUsage(2,'..','[asalin, ctemp]',\@fc,'asalin','ctemp',@_);
+ return 'nan' unless numbersp($SA,$TC);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_sigma1($SA,$TC/$TCONV);
+ }
+}
+
+{ my(@fc);
+ sub sigma2(@)
+ {
+ my($SA,$TC) = &antsFunUsage(2,'..','[asalin, ctemp]',\@fc,'asalin','ctemp',@_);
+ return 'nan' unless numbersp($SA,$TC);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_sigma2($SA,$TC/$TCONV);
+ }
+}
+
+{ my(@fc);
+ sub sigma3(@)
+ {
+ my($SA,$TC) = &antsFunUsage(2,'..','[asalin, ctemp]',\@fc,'asalin','ctemp',@_);
+ return 'nan' unless numbersp($SA,$TC);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_sigma3($SA,$TC/$TCONV);
+ }
+}
+
+{ my(@fc);
+ sub sigma4(@)
+ {
+ my($SA,$TC) = &antsFunUsage(2,'..','[asalin, ctemp]',\@fc,'asalin','ctemp',@_);
+ return 'nan' unless numbersp($SA,$TC);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_sigma4($SA,$TC/$TCONV);
+ }
+}
+
+#----------------------------------------------------------------------
+# sound speed
+# - input: asalin ctemp press
+# - VERIFIED AGAINST [libEOS83.pl]
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub sound_speed(@)
+ {
+ my($SA,$TC,$P) =
+ &antsFunUsage(3,'...','[asalin, ctemp, press]',
+ \@fc,'asalin','ctemp','press',@_);
+ return 'nan' unless numbersp($SA,$TC,$P);
+ my($TCONV) = &TCONV90(); # library uses ITS90
+ return GSW::gsw_sound_speed($SA,$TC/$TCONV,$P);
+ }
+}
+
+#----------------------------------------------------------------------
+# Coriolis parameter
+# - input: lat
+# - copied from [libEOS83.pl]
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub f(@)
+ {
+ my($lat) = &antsFunUsage(1,'f','[lat]',\@fc,'%lat',@_);
+ my($Omega) = 7.292e-5; # Gill (1982)
+ return 2 * $Omega * sin(rad($lat));
+ }
+}
+
+#----------------------------------------------------------------------
+# acceleration due to gravity
+# - input: press lat
+# - if %lat is available, it is used
+# - CHECKED AGAINST [libEOS83.pl] => DIFFERENCES LESS THAN 0.1%
+#----------------------------------------------------------------------
+
+{ my(@fc);
+ sub g(@)
+ {
+ my($P,$LAT);
+ $LAT = &antsParam('lat');
+ if (numberp($LAT)) {
+ ($P) = &antsFunUsage(1,'.','[press]',\@fc,'press',@_);
+ } else {
+ ($P,$LAT) = &antsFunUsage(2,'..','[press, lat]',\@fc,'press','lat',@_);
+ }
+ return numbersp($P,$LAT) ? GSW::gsw_grav($LAT,$P) : 'nan';
+ }
+}
+
+1;
+
--- a/libconv.pl
+++ b/libconv.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B C O N V . P L
# doc: Sat Dec 4 13:03:49 1999
-# dlm: Tue May 22 11:19:54 2018
+# dlm: Fri Feb 15 13:21:08 2019
# (c) 1999 A.M. Thurnherr
-# uE-Info: 68 41 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 529 24 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -66,6 +66,8 @@
# Jul 6, 2017: - BUG: date conversion routines did not parse 1/5/12 correctly
# Dec 18, 2017: - removed ambiguous-date warning
# May 22, 2018: - added NMEA2dec_time()
+# Jan 17, 2019: - added ISO_Datetime()
+# Feb 15, 3019: - added deg2lat, deg2lon
require "$ANTS/libEOS83.pl"; # &sigma()
require "$ANTS/libPOSIX.pl"; # &floor()
@@ -310,7 +312,7 @@
}
#----------------------------------------------------------------------
-# Decimal Time to Strin Conversion
+# Decimal Time to String Conversion
#----------------------------------------------------------------------
{ my(@fc);
@@ -376,6 +378,52 @@
}
}
+{ my(@fc);
+
+ sub ISO_Datetime(@) # day number -> ISO 8601
+ {
+
+ my($dnf); # find std dn field & epoch
+ if (@_ == 0) {
+ for (my($i)=0; $i<@antsLayout; $i++) {
+ next unless ($antsLayout[$i] =~ /^dn(\d\d)$/);
+ $dnf = $antsLayout[$i]; push(@_,$1);
+ last;
+ }
+ }
+
+ my($year,$fday) = &antsFunUsage(2,"cf","epoch, dayNo",\@fc,undef,$dnf,@_);
+
+ $year += ($year < 50) ? 2000 : 1900 # Y2K
+ if ($year < 100);
+
+ $day = int($fday); # prevent runover on last day of month
+ $fday -= $day;
+
+ while ($day > 365+&leapYearP($year)) { # adjust year
+ $day -= 365 + &leapYearP($year);
+ $year++;
+ }
+
+ my($month) = 1;
+ while ($day > &monthLength($year,$month)) {
+ $day -= &monthLength($year,$month);
+ $month++;
+ }
+
+ my($hour) = int(24*$fday);
+ $fday -= $hour/24;
+ my($min) = int(24*60*$fday);
+ $fday -= $min/24/60;
+ my($sec) = round(24*3600*$fday);
+ $min++,$sec=0 if ($sec == 60);
+ $hour++,$min=0 if ($min == 60);
+ $day++,$hour=0 if ($hour == 24);
+
+ return sprintf('%04d-%02d-%02dT%02d:%02d:%02d',$year,$month,$day,$hour,$min,$sec);
+ }
+}
+
#----------------------------------------------------------------------
# Other Misc Date Conversions
#----------------------------------------------------------------------
@@ -459,7 +507,7 @@
return ($b eq "") ? &dmh2d($deg,0,$a) : &dmh2d($deg,$a,$b);
}
-sub GMT2deg(@) # GMT degree format to decimal
+sub GMT2deg(@) # GMT degree format to decimal
{
my($GMT) = &antsFunUsage(1,".","GMT-degs ",@_);
return (substr($1,0,1) eq "-") ? $1-$2/60.0 : $1+$2/60.0
@@ -467,6 +515,22 @@
return $GMT;
}
+sub deg2lat(@) # decimal latitude to degrees:min.xx NS
+{
+ my($deg) = &antsFunUsage(1,'f','decimal latitude',@_);
+ return sprintf("%02d:%06.3f'%s",abs(int($deg)),
+ (abs($deg)-abs(int($deg)))*60,
+ $deg>=0 ? "N" : "S");
+}
+
+sub deg2lon(@) # decimal longitude to degrees:min.xx EW
+{
+ my($deg) = &antsFunUsage(1,'f','decimal longitude',@_);
+ return sprintf("%03d:%06.3f'%s",abs(int($deg)),
+ (abs($deg)-abs(int($deg)))*60,
+ $deg>=0 ? "E" : "W");
+}
+
#----------------------------------------------------------------------
# Temp-Scale Conversion
#----------------------------------------------------------------------