Merge pull request #114 from gaelforget/V0p3p1b

V0p3p1b

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

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
 
 ##

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)]
 

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))