Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Option to output H #169

Merged
merged 1 commit into from
Jan 29, 2024
Merged

Option to output H #169

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/default_parameters.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@

# OUTPUT OPTIONS
output::Bool=false # netcdf output?
output_vars::Array{String,1}=["u","v","η","sst"] # which variables to output? "du","dv","dη" also allowed
output_vars::Array{String,1}=["u","v","η","sst"] # which variables to output? "du","dv","dη","H","ζ" also allowed
output_dt::Real=24 # output time step [hours]
outpath::String=pwd() # path to output folder
compression_level::Int=3 # compression level
Expand Down
55 changes: 36 additions & 19 deletions src/output.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,12 @@ struct NcFiles
du::Union{NcFile,Nothing} # tendency of u [m^2/s^2]
dv::Union{NcFile,Nothing} # tendency of v [m^2/s^2]
dη::Union{NcFile,Nothing} # tendency of η [m^2/s]
H::Union{NcFile,Nothing} # layer thickness at rest [m]
iout::Array{Integer,1} # output index, stored in array for mutability
end

"""Generator function for "empty" NcFiles struct."""
NcFiles(x::Nothing) = NcFiles(x,x,x,x,x,x,x,x,x,[0])
NcFiles(x::Nothing) = NcFiles(x,x,x,x,x,x,x,x,x,x,[0])

"""Generator function for NcFiles struct, creating the underlying netCDF files."""
function NcFiles(feedback::Feedback,S::ModelSetup)
Expand All @@ -35,16 +36,22 @@ function NcFiles(feedback::Feedback,S::ModelSetup)
ncdu = if "du" in output_vars nc_create(x_u,y_u,"du",runpath,"m^2/s^2","zonal velocity tendency",P) else nothing end
ncdv = if "dv" in output_vars nc_create(x_v,y_v,"dv",runpath,"m^2/s^2","meridional velocity tendency",P) else nothing end
ncdη = if "dη" in output_vars nc_create(x_T,y_T,"deta",runpath,"m^2/s","sea surface height tendency",P) else nothing end
nH = if "H" in output_vars nc_create(x_T,y_T,"H",runpath,"m","undisturbed layer thickness",P,false) else nothing end

for nc in (ncu,ncv,ncη,ncsst,ncq,ncζ,ncdu,ncdv,ncdη)
if nc != nothing
if !isnothing(nc)
NetCDF.putatt(nc,"t",Dict("units"=>"s","long_name"=>"time"))
NetCDF.putatt(nc,"x",Dict("units"=>"m","long_name"=>"zonal coordinate"))
NetCDF.putatt(nc,"y",Dict("units"=>"m","long_name"=>"meridional coordinate"))
end
end

return NcFiles(ncu,ncv,ncη,ncsst,ncq,ncζ,ncdu,ncdv,ncdη,[0])
if !isnothing(nH) # output constant fields
NetCDF.putatt(nH,"x",Dict("units"=>"m","long_name"=>"zonal coordinate"))
NetCDF.putatt(nH,"y",Dict("units"=>"m","long_name"=>"meridional coordinate"))
end

return NcFiles(ncu,ncv,ncη,ncsst,ncq,ncζ,ncdu,ncdv,ncdη,nH,[0])
else
return NcFiles(nothing)
end
Expand All @@ -57,18 +64,24 @@ function nc_create( x::Array{T,1},
path::String,
unit::String,
long_name::String,
P::Parameter) where {T<:Real}
P::Parameter,
output_tdim::Bool=true) where {T<:Real}

@unpack compression_level = P

xdim = NcDim("x",length(x),values=x)
ydim = NcDim("y",length(y),values=y)
tdim = NcDim("t",0,unlimited=true)

if output_tdim
tdim = NcDim("t",0,unlimited=true)
var = NcVar(name,[xdim,ydim,tdim],t=Float32,compress=compression_level)
tvar = NcVar("t",tdim,t=Int32)
nc = NetCDF.create(joinpath(path,name*".nc"),[var,tvar],mode=NC_NETCDF4)
else
var = NcVar(name,[xdim,ydim],t=Float32,compress=compression_level)
nc = NetCDF.create(joinpath(path,name*".nc"),[var],mode=NC_NETCDF4)
end

var = NcVar(name,[xdim,ydim,tdim],t=Float32,compress=compression_level)
tvar = NcVar("t",tdim,t=Int32)

nc = NetCDF.create(joinpath(path,name*".nc"),[var,tvar],mode=NC_NETCDF4)
# add missing_value although irrelevant for ncview compatibility
NetCDF.putatt(nc,name,Dict("units"=>unit,"long_name"=>long_name,"missing_value"=>-999999f0))
return nc
Expand Down Expand Up @@ -96,51 +109,55 @@ function output_nc!(i::Int,
# As output is before copyto!(u,u0), take u0,v0,η0
# Tendencies calculate from the last time step, du = u_n+1-u_n etc
# WRITING THE VARIABLES
if ncs.u != nothing
if !isnothing(ncs.u)
@views u = Float32.(scale_inv*Diag.RungeKutta.u0[halo+1:end-halo,halo+1:end-halo])
NetCDF.putvar(ncs.u,"u",u,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.v != nothing
if !isnothing(ncs.v)
@views v = Float32.(scale_inv*Diag.RungeKutta.v0[halo+1:end-halo,halo+1:end-halo])
NetCDF.putvar(ncs.v,"v",v,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.η != nothing
if !isnothing(ncs.η)
@views η = Float32.(Diag.RungeKutta.η0[haloη+1:end-haloη,haloη+1:end-haloη])
NetCDF.putvar(ncs.η,"eta",η,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.sst != nothing
if !isnothing(ncs.sst)
@views sst = Float32.(Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty]/scale_sst)
NetCDF.putvar(ncs.sst,"sst",sst,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.q != nothing
if !isnothing(ncs.q)
@views q = Float32.(scale_inv*Diag.Vorticity.q[haloη+1:end-haloη,haloη+1:end-haloη])
NetCDF.putvar(ncs.q,"q",q,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.ζ != nothing
if !isnothing(ncs.ζ)
@unpack dvdx,dudy = Diag.Vorticity
@unpack f_q = S.grid
@views ζ = Float32.((dvdx[2:end-1,2:end-1]-dudy[2+ep:end-1,2:end-1])./abs.(f_q))
NetCDF.putvar(ncs.ζ,"relvort",ζ,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.du != nothing
if !isnothing(ncs.du)
@views u = Float32.(scale_inv*Diag.RungeKutta.u0[halo+1:end-halo,halo+1:end-halo])
@views du = u-Float32.(scale_inv*Prog.u[halo+1:end-halo,halo+1:end-halo])
NetCDF.putvar(ncs.du,"du",du,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.dv != nothing
if !isnothing(ncs.dv)
@views v = Float32.(scale_inv*Diag.RungeKutta.v0[halo+1:end-halo,halo+1:end-halo])
@views dv = v-Float32.(scale_inv*Prog.v[halo+1:end-halo,halo+1:end-halo])
NetCDF.putvar(ncs.dv,"dv",dv,start=[1,1,iout],count=[-1,-1,1])
end
if ncs.dη != nothing
if !isnothing(ncs.dη)
@views η = Float32.(Diag.RungeKutta.η0[haloη+1:end-haloη,haloη+1:end-haloη])
@views dη = η-Float32.(Prog.η[haloη+1:end-haloη,haloη+1:end-haloη])
NetCDF.putvar(ncs.dη,"deta",dη,start=[1,1,iout],count=[-1,-1,1])
end
if i==0 && !isnothing(ncs.H) # constant field, output only once, before the timestepping
@views H = Float32.(S.forcing.H[haloη+1:end-haloη,haloη+1:end-haloη])
NetCDF.putvar(ncs.H,"H",H,start=[1,1],count=[-1,-1])
end

# WRITING THE TIME
for nc in (ncs.u,ncs.v,ncs.η,ncs.sst,ncs.q,ncs.ζ,ncs.du,ncs.dv,ncs.dη)
if nc != nothing
if !isnothing(nc)
NetCDF.putvar(nc,"t",Int64[i*dtint],start=[iout])
NetCDF.sync(nc) # sync to view netcdf while model is still running
end
Expand Down
Loading