From 4186d7a50a4b3deecec3c80348846cac2409421f Mon Sep 17 00:00:00 2001 From: Max Date: Mon, 29 Jan 2024 16:39:45 +0100 Subject: [PATCH] option to output H --- src/default_parameters.jl | 2 +- src/output.jl | 55 +++++++++++++++++++++++++-------------- 2 files changed, 37 insertions(+), 20 deletions(-) diff --git a/src/default_parameters.jl b/src/default_parameters.jl index 570faa5..4ef41c3 100644 --- a/src/default_parameters.jl +++ b/src/default_parameters.jl @@ -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 diff --git a/src/output.jl b/src/output.jl index 3d271a4..a30589d 100644 --- a/src/output.jl +++ b/src/output.jl @@ -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) @@ -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 @@ -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 @@ -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