diff --git a/src/Grids.jl b/src/Grids.jl index 0f2a6a32..0e8a2289 100644 --- a/src/Grids.jl +++ b/src/Grids.jl @@ -2,17 +2,19 @@ Dict_to_NamedTuple(tmp::Dict) = (; zip(Symbol.(keys(tmp)), values(tmp))...) """ - land_mask(Γ) + land_mask(m::MeshArray) -Define land mask from `Γ.hFacC[:,1]` +Define land mask from `m` (1 if m>0; NaN if otherwise). """ -function land_mask(Γ) - μ =Γ.hFacC[:,1] +function land_mask(m::MeshArray) + μ=m μ[findall(μ.>0.0)].=1.0 μ[findall(μ.==0.0)].=NaN μ end +land_mask(Γ::NamedTuple)=land_mask(Γ.hFacC[:,1]) + """ UnitGrid(γ::gcmgrid;option="minimal") @@ -29,6 +31,10 @@ function UnitGrid(γ::gcmgrid;option="minimal") list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","hFacC","hFacS","hFacW"); list_u=(u"m",u"m",u"m",u"m",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0,1.0) list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5]) + elseif option=="light" + list_n=("XC","XG","YC","YG","RAC","DXC","DXG","DYC","DYG","Depth"); + list_u=(u"m",u"m",u"m",u"m",u"m^2",u"m",u"m",u"m",u"m",u"m") + list_p=(pc,pg,pc,pg,pc,pu,pv,pv,pu,pc) else list_n=("XC","YC"); list_u=(u"°",u"°") @@ -231,6 +237,12 @@ function GridLoad(γ::gcmgrid;option="minimal") end list_n=(list_n...,"DRC","DRF","RC","RF") list_n=(list_n...,"hFacC","hFacS","hFacW"); + elseif option=="light" + list_n=("XC","XG","YC","YG","RAC","DXC","DXG","DYC","DYG","Depth"); + if (!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path)))) + list_n=(list_n...,"AngleCS","AngleSN"); + end + list_n=(list_n...,"DRC","DRF","RC","RF") else list_n=("XC","YC"); end diff --git a/src/Interpolation.jl b/src/Interpolation.jl index 45f56b52..266e01cd 100644 --- a/src/Interpolation.jl +++ b/src/Interpolation.jl @@ -193,8 +193,8 @@ InterpolationFactors(Γ,lon::Number,lat::Number) = InterpolationFactors(Γ,[lon] Read e.g. `interp_coeffs_halfdeg.jld2` ``` -file_int=MeshArrays.interpolation_setup() -λ=MeshArrays.interpolation_setup(file_int) +fil=joinpath(tempdir(),"interp_coeffs_halfdeg.jld2") +λ=MeshArrays.interpolation_setup(fil) ``` """ function interpolation_setup(fil::String) @@ -203,21 +203,30 @@ function interpolation_setup(fil::String) end """ - interpolation_setup(;path=tempdir()) + interpolation_setup(;Γ,lon,lat,path,url) -Download e.g. `interp_coeffs_halfdeg.jld2` +Download or recompute interpolation coefficients (e.g. `interp_coeffs_halfdeg.jld2`). -``` -file_int=MeshArrays.interpolation_setup() -λ=MeshArrays.interpolation_setup(file_int) -``` +- `λ=interpolation_setup()` to download from `url` to `path` +- `λ=interpolation_setup(Γ=Γ)` to recompute at `lon,lat` """ -function interpolation_setup(; +function interpolation_setup(;Γ=NamedTuple(), + lon=[i for i=-179.:2.0:179., j=-89.:2.0:89.], + lat=[j for i=-179.:2.0:179., j=-89.:2.0:89.], path=tempdir(), url="https://zenodo.org/record/5784905/files/interp_coeffs_halfdeg.jld2") - fil=joinpath(tempdir(),basename(url)) - !isfile(fil) ? MeshArrays.download_file(url, fil) : nothing - fil + + if isempty(Γ) + fil=joinpath(path,basename(url)) + !isfile(fil) ? MeshArrays.download_file(url, fil) : nothing + else + lon= + lat= + (f,i,j,w)=InterpolationFactors(Γ,vec(lon),vec(lat)) + fil=tempname()*"_interp_coeffs.jld2" + MeshArrays.write_JLD2(fil; lon=lon, lat=lat, f=f, i=i, j=j, w=w) + end + interpolation_setup(fil) end ## diff --git a/src/Operations.jl b/src/Operations.jl index 9fab322a..807a934c 100644 --- a/src/Operations.jl +++ b/src/Operations.jl @@ -79,7 +79,7 @@ function curl(u::MeshArray,v::MeshArray,Γ::NamedTuple) uvcurl=similar(Γ.XC) fac=exchange(1.0 ./Γ.RAZ,1) (U,V)=exchange(u,v,1) - (DXC,DYC)=exchange(Γ.DXC,Γ.DYC,1) + (DXC,DYC)=exchange(Γ.DXC,Γ.DYC,1) [DXC[i].=abs.(DXC[i]) for i in eachindex(U)] [DYC[i].=abs.(DYC[i]) for i in eachindex(V)] diff --git a/test/runtests.jl b/test/runtests.jl index ef266c87..0ef550cd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,20 +54,23 @@ end γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90) Tx=γ.read(MeshArrays.GRID_LLC90*"TrspX.bin",MeshArray(γ,Float32)) Ty=γ.read(MeshArrays.GRID_LLC90*"TrspY.bin",MeshArray(γ,Float32)) - Γ=GridLoad(γ;option="full") + Γ=GridLoad(γ;option="light") - μ =Γ.hFacC[:,1] - μ[findall(μ.>0.0)].=1.0 - μ[findall(μ.==0.0)].=NaN + hFacC=GridLoadVar("hFacC",γ) + μ=land_mask(hFacC[:,1]) lons=[-68 -63]; lats=[-54 -66]; name="Drake Passage" Trsct=Transect(name,lons,lats,Γ) #Various vector operations - U=0*Γ.hFacW; V=0*Γ.hFacS; + hFacW=GridLoadVar("hFacW",γ) + hFacS=GridLoadVar("hFacS",γ) + RAZ=GridLoadVar("RAZ",γ) + + U=0*hFacW; V=0*hFacS; UVtoTransport(U,V,Γ) UVtoUEVN(U[:,1],V[:,1],Γ) - curl(U[:,1],V[:,1],Γ) + curl(U[:,1],V[:,1], merge(Γ,(hFacW=hFacW,hFacS=hFacS,RAZ=RAZ,)) ) #Meridional transport integral uv=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"]) @@ -90,7 +93,7 @@ end TxR = Tx-TxD TyR = Ty-TyD #vector Potential - TrspPsi=VectorPotential(TxR,TyR,Γ); + TrspPsi=VectorPotential(TxR,TyR, merge(Γ,(hFacC=hFacC,hFacW=hFacW,hFacS=hFacS)) ); tst=extrema(write(TrspPsi))