Unnamed Fossil Project

    if (FIG)
        using Colors
    end
    include("/home/ezaron/opt/julia/edzlib.jl")
    
    nodes=128
    # S2+ all other constituents are in here:
    workdir = "/scratch/ezaron/workdir/"
    # M2 is in here, incorrectly named output files are Xout_Zhao.nc:
#    workdir = "/scratch/ezaron/workdir12.2/"
    rearth = 6.371e3
    dtr = pi/180
    rtd = 180/pi

#    fn = "/Xout_Zhao.nc"
#    fn = "/XoutM2.nc"
    fn = "/XoutS2.nc"
#    fn = "/XoutK1.nc"
#    fn = "/XoutO1.nc"
#    fn = "/XoutN2.nc"
#    fn = "/XoutMA2.nc"
#    fn = "/XoutMB2.nc"
end

# coarse grain scalars from each solve:
dlonc = 1.0
dlatc = 1.0
# NEXT TIME YOU DO A COMPLETE SET, FIX THIS BUG:
#      We should only go to 360-dlonc+epsilon():
#clon = collect(0:dlonc:360.0)
clon = collect(0:dlonc:360.0-dlonc+eps()) # Fixed!
ind = find( x -> x > 180.0 , clon )
clonx= copy(clon)
clonx[ind] = clonx[ind] .- 360.0
clat = collect(-66:dlatc:66)
nxc = length(clon)
nyc = length(clat)
is = ["kwoa",

fout = workdir * "HRET12.4.nc" # K1
fout = workdir * "HRET12.5.nc" # O1
fout = workdir * "HRET12.6.nc" # N2
fout = workdir * "HRET12.7.nc" # MA2
fout = workdir * "HRET12.8.nc" # MB2

fout = workdir * "HRET13.1.nc" # M2 -- no testing for small HA, and use bootstrap for explained variance.
fout = workdir * "HRET13.2.nc" # S2

run(`rm -f $fout`)
nccreate(fout,"res12","longitude",glon,"latitude",glat,"nc",nc,
         atts=Dict("longname" => "tide from HA tuned to GM",
                   "units"    => "m"))
nccreate(fout,"res11","longitude","latitude","nc",
         atts=Dict("longname" => "tide from HA tuned to HA",

# solutions -- for computing specific mode solutions, or solutions with wavenumebers all in a specific
# direction or whatever.

@sync @everywhere begin
    using Distributed
    using SharedArrays
    using Dates
    using DSP
    using GSL
    using FFTW
    using LinearAlgebra
    using LinearMaps
    using IterativeSolvers
    using CRC32c
    using SpecialFunctions
    using ColorSchemes

    # To get graphics from root proc when running interactively:
    #FIG=true
    FIG=false
    if (FIG)
        using Colors
    end
    include("/home/ezaron/opt/julia/edzlib.jl")
    
    nodes=128
    # S2+ all other constituents are in here:
#    workdir = "/scratch/ezaron/workdir/"
    # M2 is in here, incorrectly named output files are Xout_Zhao.nc:
    workdir = "/scratch/ezaron/workdir12.2/"
    rearth = 6.371e3
    dtr = pi/180
    rtd = 180/pi

    nc = 2
    nmodes = 2
    cid = "M2"
#    fn = "/Xout_Zhao.nc"
    fn = "/Xout" * cid * ".nc"
#    fn = "/XoutS2.nc"
#    fn = "/XoutK1.nc"
#    fn = "/XoutO1.nc"

# coarse grain scalars from each solve:
dlonc = 1.0
dlatc = 1.0
# NEXT TIME YOU DO A COMPLETE SET, FIX THIS BUG:
#      We should only go to 360-dlonc+epsilon():
#clon = collect(0:dlonc:360.0)
clon = collect(0:dlonc:360.0-dlonc+eps()) # Fixed!
ind = find( x -> x > 180.0 , clon )
clonx= copy(clon)
clonx[ind] = clonx[ind] .- 360.0
clat = collect(-66:dlatc:66)
nxc = length(clon)
nyc = length(clat)
# Kludgy:

# coarse scalars:
csca  = SharedArray{Float64}(nxc,nyc,ns,init = S -> S[localindices(S)] .= 0.0)
# coarse weight:
cwgt  = SharedArray{Float64}(nxc,nyc,   init = S -> S[localindices(S)] .= 0.0)

@sync @everywhere function doer!(myj,
                                 glon,glonx,glat,gwgt,
                                 gres12,gres32,gresx,
                                 clon,clonx,clat,cwgt,csca,is)
    
#    println("doer has been called with myid, myj =",myid()," ",myj)

    fin = myj * "/lonlatL.nc"
    lon0 = ncread(fin,"lon0")[1]
    lat0 = ncread(fin,"lat0")[1]
    Ldata0 = ncread(fin,"Ldata")[1]
    fin = myj * fn
    Ldata = ncgetatt(fin,"global","Ldata")
    dx    = ncgetatt(fin,"global","dx")
    Lplane= ncgetatt(fin,"global","Lplane")
    cidv  = ncgetatt(fin,"Global","cidv")
    nc    = length(cidv)
    if ( Ldata == nothing )
        println("ERROR: empty Ldata in Xout.nc on ",myj)
        # oplot([lon0],[lat0],"kx",markersize=0.25)
        return "rx"
    end

    rdx = 1.0 ./dx
    n,ix,iy,x1,y1,ilon,ilat,dlon,dlat,lon_lim,lat_lim,lon0,lat0 = geography_init(lon0,lat0,Lplane,dx)

    ilon = wrapto(ilon,lon0)

    function xy2ll(x,y)
        ilon = lon0 .+ x/(rearth*dtr*cos(lat0*dtr))
        ilat = lat0 .+ y/(rearth*dtr)

             :cori => cori,
             :mbasis => mbasis, :mbas => mbas,
             :ffttmp1 => ffttmp1, :ffttmp2 => ffttmp2, :fftp1 => fftp1, :fftq1 => fftq1)

    kwoa = ncgetatt(fin,"Global","kwoa")
    #    ampvec = ncgetatt(fin,"Global","ampvec")
    bwvec  = ncgetatt(fin,"Global","bwvec")
    acoef  = ncgetatt(fin,"Global","acoef")

    if (nmodes > length(kwoa))
        println("ERROR: too many modes sought by patchmerge_special.")
        println("ABORTING")
        stophere
    end

    B2 = ncvarget(fin,"B2")
    B2ind = []
    for k=1:round(Int,nc/2)
        b2 = view(B2,:,:,k)
        ind = find( x -> x != 0.0, b2)
        push!(B2ind,ind)
    end

    bvec12 = squeeze(ncvarget(fin,"bvec",start=[1,1,2],count=[-1,1,1]))
    bvec32 = squeeze(ncvarget(fin,"bvec",start=[1,3,2],count=[-1,1,1]))
    cof12 = bvec2cof(bvec12,B2ind,G)
    cof32 = bvec2cof(bvec32,B2ind,G)

    for nmod=1:nmodes
        
    kx,lx = psd2dkl(x1,y1,0.0)
    k2 = ones(iy,1)*kx' ;
    l2 = lx*ones(1,ix);
    kk2 = sqrt.( k2.^2 + l2.^2 );
    # fftshift these 2d arrays into the same rotation as the B2 array so that
    # ind will be able to access into either B2 or the cof arrays:
    k2 = ifftshift(k2)
    l2 = ifftshift(l2)

        cof12t[:,:,p,1] = cof12[:,:,p,1].*B2m[:,:]
        cof12t[:,:,p,2] = cof12[:,:,p,2].*B2m[:,:]
        cof32t[:,:,p,1] = cof32[:,:,p,1].*B2m[:,:]
        cof32t[:,:,p,2] = cof32[:,:,p,2].*B2m[:,:]
    end
    res12 = innerfunb(cof12t,G)
    res32 = innerfunb(cof32t,G)
    resx  = acoef*res12 + (1.0 - acoef)*res32

    # Need to permute these like the original file:
    res12 = permutedims(res12,(2,1,3))
    res32 = permutedims(res32,(2,1,3))
    resx  = permutedims(resx ,(2,1,3))
    n1,m1,nc = size(res12)