Using gsw_nsquared

[mukund-gupta]

Hello, I am trying to use gsw_nsquared as follows:

using GibbsSeaWater

# Constants
g = 9.81
Nz = 64
Δz = 2.5
Lz = Δz*Nz
z = -((Δz/2):Δz:(Lz-Δz/2))
Stop = 34.1 
Sbot = 34.7
Ttop = -0.5
Tbot = 0
ztop = -40
zbot = -Lz
dSdz = (Sbot - Stop)/(zbot - ztop)
dTdz = (Tbot - Ttop)/(zbot - ztop)

# Profiles
S = Stop .+ dSdz*(z .< ztop).*(z .- ztop)
T = Ttop .+ dTdz*(z .< ztop).*(z .- ztop)
Pres = -rho_O*g*z*1e-4
N2 = float(Nz-1)
Pmid = float(Nz-1)
a = gsw_nsquared.(S, T, Pres, 0, Nz, N2, Pmid)

But get the following error:

 ERROR: LoadError: MethodError: no method matching unsafe_convert(::Type{Ptr{Float64}}, ::Float64)
Closest candidates are:
  unsafe_convert(::Type{T}, ::T) where T<:Ptr at essentials.jl:391
  unsafe_convert(::Type{Ptr{T}}, ::Ptr{Tuple{Vararg{T,N}}}) where {N, T} at refpointer.jl:136
  unsafe_convert(::Type{P}, ::Ptr) where P<:Ptr at essentials.jl:392
  ...
Stacktrace:
 [1] gsw_nsquared(::Float64, ::Float64, ::Float64, ::Int64, ::Int64, ::Float64, ::Float64) at /Users/guptam/.julia/packages/GibbsSeaWater/vq3GE/src/gen_gswteos10.jl:370
 [2] _broadcast_getindex_evalf at ./broadcast.jl:648 [inlined]
 [3] _broadcast_getindex at ./broadcast.jl:621 [inlined]
 [4] getindex at ./broadcast.jl:575 [inlined]
 [5] copy at ./broadcast.jl:876 [inlined]
 [6] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(gsw_nsquared),Tuple{Array{Float64,1},Array{Float64,1},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64,Int64,Float64,Float64}}) at ./broadcast.jl:837
 [7] top-level scope at compare_nsquared.jl:137
 [8] include(::String) at ./client.jl:457
 [9] top-level scope at REPL[14]:1
in expression starting at compare_nsquared.jl:137

Could you advise how I should use gsw_nsquared instead? Many thanks.

[kouketsu]

Now we did not have wrapper functions easily to use some subroutines including outputs in the arguments (ex. gsw_nsquared). The functions are confusing. Sorry @mukund-gupta.

In the released version, please try below.

using GibbsSeaWater

# Constants
g = 9.81
Nz = 64
Δz = 2.5
Lz = Δz*Nz
z = collect(-((Δz/2):Δz:(Lz-Δz/2))) # to Array{Float64} for Pres
Stop = 34.1 
Sbot = 34.7
Ttop = -0.5
Tbot = 0
ztop = -40
zbot = -Lz
dSdz = (Sbot - Stop)/(zbot - ztop)
dTdz = (Tbot - Ttop)/(zbot - ztop)
rho_O = 1000.0 # !!! I miss this value in the original script
Lat0 = 0.0 # latitude where the profile is measured

# Profiles
S = Stop .+ dSdz*(z .< ztop).*(z .- ztop)
T = Ttop .+ dTdz*(z .< ztop).*(z .- ztop)
Pres = -rho_O*g*z*1e-4
N2    = zeros(Nz-1) # Array for output of Brunt-Vaisalla values @ Pmid) 
Pmid = zeros(Nz-1) # Array for outputs of pressure for N2
Lats = fill(Lat0, Nz) 
a = gsw_nsquared(S, T, Pres, Lats, Nz, N2, Pmid)

Or

In the master repository (pkg> add GibbsSeaWater#master), you can use

N2, Pmid = gsw_nsquared(S, T, Pres, Lats)

[mukund-gupta]

Thank you very much! Both methods work for me. I'll just point out that when using

N2, Pmid = gsw_nsquared(S, T, Pres, Lats),

N2 and Pmid have the same dimensions as T and S, and that the last element of those vectors are zeros. I don't know whether this is intentional or not.

[kouketsu]

Sorry, the lengths of outputs are wrong. I fixed them. Thank you for pointing out @mukund-gupta :).

[mukund-gupta]

Great! Thanks for fixing that.

I would personally have Lats as an optional argument, since that would make it consistent with the Matlab and python versions of gsw_nsquared.

Thanks for the help!