.
#======================================================================
# 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;