libGM.pl
changeset 4 ff72b00b4342
parent 3 55a8c407d38e
child 7 f7b4692ad805
child 14 4097040ad31c
equal deleted inserted replaced
3:55a8c407d38e 4:ff72b00b4342
     1 #======================================================================
     1 #======================================================================
     2 #                    L I B G M . P L 
     2 #                    L I B G M . P L 
     3 #                    doc: Sun Feb 20 14:43:47 2011
     3 #                    doc: Sun Feb 20 14:43:47 2011
     4 #                    dlm: Fri Sep 14 12:35:40 2012
     4 #                    dlm: Sat Apr 20 20:05:51 2013
     5 #                    (c) 2011 A.M. Thurnherr
     5 #                    (c) 2011 A.M. Thurnherr
     6 #                    uE-Info: 61 0 NIL 0 0 72 2 2 4 NIL ofnI
     6 #                    uE-Info: 62 0 NIL 0 0 72 2 2 4 NIL ofnI
     7 #======================================================================
     7 #======================================================================
     8 
     8 
     9 # HISTORY:
     9 # HISTORY:
    10 #	Feb 20, 2011: - created
    10 #	Feb 20, 2011: - created
    11 #	Feb 28, 2011: - cosmetics
    11 #	Feb 28, 2011: - cosmetics
    15 #					affects only long wavelengths
    15 #					affects only long wavelengths
    16 #				  - return nan for omega outside internal-wave band
    16 #				  - return nan for omega outside internal-wave band
    17 #	Mar 29, 2012: - re-wrote using definition of B(omega) from Munk (1981)
    17 #	Mar 29, 2012: - re-wrote using definition of B(omega) from Munk (1981)
    18 #	Aug 23, 2012: - cosmetics?
    18 #	Aug 23, 2012: - cosmetics?
    19 #	Sep  7, 2012: - made N0, E0, b, jstar global
    19 #	Sep  7, 2012: - made N0, E0, b, jstar global
       
    20 #	Dec 28, 2012: - added allowance for small roundoff error to Sw()
    20 
    21 
    21 require "$ANTS/libEOS83.pl";
    22 require "$ANTS/libEOS83.pl";
    22 
    23 
    23 my($pi) = 3.14159265358979;
    24 my($pi) = 3.14159265358979;
    24 
    25 
    56 
    57 
    57 sub m($$)	# vertical wavenumber as a function of mode number & stratification params
    58 sub m($$)	# vertical wavenumber as a function of mode number & stratification params
    58 {
    59 {
    59 	my($j,$N,$omega) = @_;
    60 	my($j,$N,$omega) = @_;
    60 
    61 
    61 	return defined($omega)
    62 	return defined($omega) && ($omega <= $GM_N0)
    62 		   ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j
    63 		   ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j
    63 		   : $pi * $j * $N / ($GM_b * $GM_N0);			# valid, except in vicinity of buoyancy turning frequency (p. 285)
    64 		   : $pi * $j * $N / ($GM_b * $GM_N0);			# valid, except in vicinity of buoyancy turning frequency (p. 285)
    64 }
    65 }
    65 
    66 
    66 sub B($)											# structure function (omega dependence)
    67 sub B($)											# structure function (omega dependence)
    75 sub Sw($$$$)
    76 sub Sw($$$$)
    76 {
    77 {
    77 	my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_);
    78 	my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_);
    78 
    79 
    79 	local($f) = abs(&f($lat));
    80 	local($f) = abs(&f($lat));
       
    81 	$omega += $PRACTICALLY_ZERO if ($omega < $f);
       
    82 	$omega -= $PRACTICALLY_ZERO if ($omega > $N);
    80 	return nan if ($omega < $f || $omega > $N);
    83 	return nan if ($omega < $f || $omega > $N);
    81 
    84 
    82 	my($GM_b) = 1300; #m								# pycnocline lengthscale
    85 	my($GM_b) = 1300; #m								# pycnocline lengthscale
    83 
    86 
    84 	my($mstar) = &m($GM_jstar,$N,$omega);
    87 	my($mstar) = &m($GM_jstar,$N,$omega);