convert sigma coordinate to depth coordinate

[fdaryabor]

Hi Everybody
ROMS is in xi,eta,sigma
coordiantes. I want to convert sigma co-ordinates to
depth coordinates. I will get the output nc file in
the following form. Everybody could help me for guide or have script for converting thanks so much.
netcdf roms_avg_Y1M1 {
dimensions:
xi_rho = 103 ;
xi_u = 102 ;
eta_rho = 104 ;
eta_v = 103 ;
s_rho = 20 ;
s_w = 21 ;
time = UNLIMITED ; // (30 currently)
auxil = 4 ;
variables:
char spherical ;
spherical:long_name = "grid type logical switch" ;
spherical:option_T = "spherical" ;
spherical:option_F = "cartesian" ;
float xl ;
xl:long_name = "domain length in the XI-direction" ;
xl:units = "meter" ;
float el ;
el:long_name = "domain length in the ETA-direction" ;
el:units = "meter" ;
float h(eta_rho, xi_rho) ;
h:long_name = "bathymetry at RHO-points" ;
h:units = "meter" ;
h:field = "bath, scalar" ;
float f(eta_rho, xi_rho) ;
f:long_name = "Coriolis parameter at RHO-points" ;
f:units = "second-1" ;
f:field = "coriolis, scalar" ;
float pm(eta_rho, xi_rho) ;
pm:long_name = "curvilinear coordinate metric in XI" ;
pm:units = "meter-1" ;
pm:field = "pm, scalar" ;
float pn(eta_rho, xi_rho) ;
pn:long_name = "curvilinear coordinate metric in ETA" ;
pn:units = "meter-1" ;
pn:field = "pn, scalar" ;
float lon_rho(eta_rho, xi_rho) ;
lon_rho:long_name = "longitude of RHO-points" ;
lon_rho:units = "degree_east" ;
lon_rho:field = "lon_rho, scalar" ;
float lat_rho(eta_rho, xi_rho) ;
lat_rho:long_name = "latitude of RHO-points" ;
lat_rho:units = "degree_north" ;
lat_rho:field = "lat_rho, scalar" ;
float angle(eta_rho, xi_rho) ;
angle:long_name = "angle between XI-axis and EAST" ;
angle:units = "radians" ;
angle:field = "angle, scalar" ;
float mask_rho(eta_rho, xi_rho) ;
mask_rho:long_name = "mask on RHO-points" ;
mask_rho:option_0 = "land" ;
mask_rho:option_1 = "water" ;
int time_step(time, auxil) ;
time_step:long_name = "time step and record numbers from initialization" ;
float scrum_time(time) ;
scrum_time:long_name = "averaged time since initialization" ;
scrum_time:units = "second" ;
scrum_time:field = "time, scalar, series" ;
float zeta(time, eta_rho, xi_rho) ;
zeta:long_name = "averaged free-surface" ;
zeta:units = "meter" ;
zeta:field = "free-surface, scalar, series" ;
float ubar(time, eta_rho, xi_u) ;
ubar:long_name = "averaged vertically integrated u-momentum component" ;
ubar:units = "meter second-1" ;
ubar:field = "ubar-velocity, scalar, series" ;
float vbar(time, eta_v, xi_rho) ;
vbar:long_name = "averaged vertically integrated v-momentum component" ;
vbar:units = "meter second-1" ;
vbar:field = "vbar-velocity, scalar, series" ;
float u(time, s_rho, eta_rho, xi_u) ;
u:long_name = "averaged u-momentum component" ;
u:units = "meter second-1" ;
u:field = "u-velocity, scalar, series" ;
float v(time, s_rho, eta_v, xi_rho) ;
v:long_name = "averaged v-momentum component" ;
v:units = "meter second-1" ;
v:field = "v-velocity, scalar, series" ;
float temp(time, s_rho, eta_rho, xi_rho) ;
temp:long_name = "averaged potential temperature" ;
temp:units = "Celsius" ;
temp:field = "temperature, scalar, series" ;
float salt(time, s_rho, eta_rho, xi_rho) ;
salt:long_name = "averaged salinity" ;
salt:units = "PSU" ;
salt:field = "salinity, scalar, series" ;
float omega(time, s_w, eta_rho, xi_rho) ;
omega:long_name = "averaged S-coordinate vertical momentum component" ;
omega:units = "meter second-1" ;
omega:field = "omega, scalar, series" ;
float w(time, s_rho, eta_rho, xi_rho) ;
w:long_name = "averaged vertical momentum component" ;
w:units = "meter second-1" ;
w:field = "w-velocity, scalar, series" ;
float AKt(time, s_w, eta_rho, xi_rho) ;
AKt:long_name = "averaged temperature vertical diffusion coefficient" ;
AKt:units = "meter2 second-1" ;
AKt:field = "AKt, scalar, series" ;
float hbl(time, eta_rho, xi_rho) ;
hbl:long_name = "averaged depth of planetary boundary layer" ;
hbl:units = "meter" ;
hbl:field = "hbl, scalar, series" ;
float hbbl(time, eta_rho, xi_rho) ;
hbbl:long_name = "averaged depth of bottom boundary layer" ;
hbbl:units = "meter" ;
hbbl:field = "hbbl, scalar, series" ;
float bostr(time, eta_rho, xi_rho) ;
bostr:long_name = "averaged Kinematic bottom stress" ;
bostr:units = "N/m2" ;

// global attributes:
:type = "ROMS averages file" ;
:title = "Southern Benguela" ;
:date = ;
:rst_file = "ROMS_FILES/roms_rst.nc" ;
:his_file = "ROMS_FILES/roms_his.0000.nc" ;
:avg_file = "ROMS_FILES/roms_avg.0000.nc" ;
:grd_file = "ROMS_FILES/roms_grd.nc" ;
:ini_file = "ROMS_FILES/roms_ini.nc" ;
:frc_file = "ROMS_FILES/roms_frc.nc" ;
:theta_s = 5.f ;
:theta_s_expl = "S-coordinate surface control parameter" ;
:theta_b = 0.4f ;
:theta_b_expl = "S-coordinate bottom control parameter" ;
:Tcline = 10.f ;
:Tcline_expl = "S-coordinate surface/bottom layer width" ;
:Tcline_units = "meter" ;
:hc = 10.f ;
:hc_expl = "S-coordinate parameter, critical depth" ;
:hc_units = "meter" ;
:sc_w = -1.f, -0.95f, -0.9f, -0.85f, -0.8f, -0.75f, -0.7f, -0.65f, -0.6f, -0.55f, -0.5f, -0.45f, -0.4f, -0.35f, -0.3f, -0.25f, -0.2f, -0.15f, -0.1f, -0.05f, 0.f ;
:sc_w_expl = "S-coordinate at W-points" ;
:Cs_w = -1.f, -0.8655258f, -0.7593114f, -0.6742046f, -0.6041494f, -0.5437741f, -0.4881475f, -0.4328656f, -0.3746809f, -0.3126323f, -0.2489214f, -0.188284f, -0.1356491f, -0.09380978f, -0.06283176f, -0.04099445f, -0.02601683f, -0.01581968f, -0.008792158f, -0.003783539f, -2.220446e-17f ;
:Cs_w_expl = "S-coordinate stretching curves at W-points" ;
:sc_r = -0.975f, -0.925f, -0.875f, -0.825f, -0.775f, -0.725f, -0.675f, -0.625f, -0.575f, -0.525f, -0.475f, -0.425f, -0.375f, -0.325f, -0.275f, -0.225f, -0.175f, -0.125f, -0.075f, -0.025f ;
:sc_r_expl = "S-coordinate at W-points" ;
:Cs_r = -0.9287272f, -0.8093643f, -0.7145184f, -0.6376433f, -0.5730692f, -0.5156592f, -0.4607052f, -0.4042644f, -0.3440788f, -0.2807264f, -0.2178818f, -0.1607245f, -0.1133185f, -0.07704455f, -0.05091174f, -0.0327879f, -0.02043305f, -0.01199137f, -0.006092094f, -0.001777454f ;
:Cs_r_expl = "S-coordinate stretching curves at RHO-points" ;
:ntimes = 216000 ;
:ndtfast = 48 ;
:dt = 720.f ;
:dtfast = 15.f ;
:nwrt = 120 ;
:ntsavg = 1 ;
:ntsavg_expl = "starting time-step for accumulation of time-averaged fields" ;
:navg = 120 ;
:navg_expl = "number of time-steps between time-averaged records" ;
:visc2 = 0.f ;
:visc2_expl = "Laplacian mixing coefficient for momentum" ;
:visc2_units = "meter2 second-1" ;
:tnu2 = 0.f ;
:tnu2_expl = "Laplacian mixing coefficient for tracers" ;
:tnu2_units = "meter2 second-1" ;
:Zob = 0.01f ;
:Zob_expl = "VonKarman/Prandtl log layer : roughness scale" ;
:Zob_units = "meter" ;
:Cdb_max = 0.1f ;
:Cdb_min = 0.0001f ;
:Cdb_expl = "Range of quadratic drag coefficient\000M" ;
:rho0 = 1025.f ;
:rho0_expl = "Mean density used in Boussinesq approximation" ;
:rho0_units = "kilogram meter-3" ;
:gamma2 = 1.f ;
:gamma2_expl = "Slipperiness parameter" ;
:x_sponge = 200000.f ;
:v_sponge = 800.f ;
:sponge_expl = "Sponge parameters : extent (m) & viscosity (m2.s-1)" ;
:SRCS = "main.F step.F read_inp.F timers.F init_scalars.F init_arrays.F set_weights.F set_scoord.F ana_grid.F setup_grid1.F setup_grid2.F set_nudgcof.F ana_initial.F analytical.F setup_kwds.F check_kwds.F check_srcs.F check_switches1.F check_switches2.F zoom.F set_nudgcof_fine.F step2d.F u2dbc.F v2dbc.F zetabc.F obc_volcons.F pre_step3d.F step3d_t.F step3d_uv1.F step3d_uv2.F prsgrd.F rhs3d.F set_depth.F omega.F uv3dmix.F t3dmix.F smagorinsky.F u3dbc.F v3dbc.F t3dbc.F rho_eos.F ab_ratio.F alfabeta.F ana_vmix.F bvf_mix.F lmd_vmix.F lmd_skpp.F lmd_bkpp.F lmd_swfrac.F lmd_wscale.F diag.F wvlcty.F checkdims.F grid_stiffness.F bio_diag.F get_date.F lenstr.F closecdf.F ana_initracer.F insert_node.F nf_fread.F get_grid.F get_initial.F def_grid.F def_his.F def_rst.F set_avg.F wrt_grid.F wrt_his.F wrt_rst.F wrt_avg.F def_diags.F wrt_diags.F output.F wrt_diags_avg.F set_diags_avg.F put_global_atts.F get_vbc.F set_cycle.F get_wwave.F get_tclima.F get_uclima.F get_ssh.F get_sss.F get_smflux.F get_stflux.F get_srflux.F get_sst.F get_tides.F clm_tides.F get_bulk.F bulk_flux.F get_bry.F nf_read_bry.F MPI_Setup.F MessPass2D.F MessPass3D.F exchange.F biology.F sediment.F bbl.F init_floats.F wrt_floats.F step_floats.F rhs_floats.F interp_rho.F def_floats.F init_arrays_floats.F random_walk.F get_initial_floats.F init_sta.F wrt_sta.F step_sta.F interp_sta.F def_sta.F init_arrays_sta.F AMRDIR AGRIFZOOM/AGRIF_YOURFILES" ;
:CPP-options = "REGIONAL BENGUELA OPENMP TIDES OBC_EAST OBC_NORTH OBC_SOUTH SOLVE3D UV_COR UV_ADV SSH_TIDES UV_TIDES TIDERAMP CURVGRID SPHERICAL MASKING AVERAGES AVERAGES_K SALINITY NONLIN_EOS SPLIT_EOS QCORRECTION SFLX_CORR DIURNAL_SRFLUX BULK_FLUX BULK_EP BULK_SMFLUX BULK_WVEC BULK_WSTR SPONGE CLIMATOLOGY ZCLIMATOLOGY M2CLIMATOLOGY M3CLIMATOLOGY TCLIMATOLOGY ZNUDGING M2NUDGING M3NUDGING TNUDGING FRC_BRY Z_FRC_BRY M2_FRC_BRY M3_FRC_BRY T_FRC_BRY ANA_BSFLUX ANA_BTFLUX UV_VIS2 MIX_GP_UV TS_DIF2 MIX_GP_TS CLIMAT_TS_MIXH LMD_MIXING LMD_SKPP LMD_BKPP LMD_RIMIX LMD_CONVEC OBC_M2FLATHER OBC_M3ORLANSKI OBC_TORLANSKI M2FILTER_FLAT" ;
}

[staalstrom]

Hallo

Here's a matlab script to convert from sigma to z coordinates.

% USAGE: z=zdepth_roms(h,hc,theta,b,N)
% z z-coordinate center of layer (m)
% h depth (m)
% hc = 20;
% theta = 5;
% b = 0.6;
% N number of layers
function z=zdepth_roms(h,hc,theta,b,N)
% Center of layers
ds=1/N;
s = -1+ds/2:ds:0-ds/2;
A = sinh(theta*s)/sinh(theta);
B = (tanh(theta*(s+.5))-tanh(theta*.5))/(2*tanh(theta*.5));
C = (1-b)*A + b*B;
z = hc*s + (h-hc)*C;

Andre Staalstrom

[fdaryabor]

Hi my friend
Thanks so much for your guide.
I have model output to following form for 5 years.
roms_avg_Y1M1,roms_avg_Y1M2,...,roms_avg_Y1M12
roms_avg_Y2M1,roms_avg_Y2M2,...,roms_avg_Y2M12
roms_avg_Y3M1,roms_avg_Y3M2,...,roms_avg_Y3M12
roms_avg_Y4M1,roms_avg_Y4M2,...,roms_avg_Y4M12
roms_avg_Y5M1,roms_avg_Y5M2,...,roms_avg_Y5M12
for example size(roms_avg_Y1M1)

ans =

30 20 104 102

time record for first month from model of 1 year=30(30 day)
N number of layers
eta_rho = 104 and xi_u = 102

Totally for per file we have 13 variables and have following metrics

function [n3dvars,varcell,L,M,N]=get_3dvars(nc)

n3dvars =

13


varcell =

'zeta'
'ubar'
'vbar'
'u'
'v'
'temp'
'salt'
'omega'
'w'
'AKt'
'hbl'
'hbbl'
'bostr'


L =

103 102 103 102 103 103 103 103 103 103 103 103 103


M =

104 104 103 104 103 104 104 104 104 104 104 104 104


N =

1 1 1 20 20 20 20 21 20 21 1 1 1

I want computing mean monthly climatology for first month from 1 year until 5 year,second month from 1 year until 5 year,...,12th month from 1 year until 5 year. could you please guide or help me for a matlab script for computing it.
Best Regards
F.Daryabor

[magaesco]

Hallo

Here's a matlab script to convert from sigma to z coordinates.

% USAGE: z=zdepth_roms(h,hc,theta,b,N)
% z z-coordinate center of layer (m)
% h depth (m)
% hc = 20;
% theta = 5;
% b = 0.6;
% N number of layers
function z=zdepth_roms(h,hc,theta,b,N)
% Center of layers
ds=1/N;
s = -1+ds/2:ds:0-ds/2;
A = sinh(theta*s)/sinh(theta);
B = (tanh(theta*(s+.5))-tanh(theta*.5))/(2*tanh(theta*.5));
C = (1-b)*A + b*B;
z = hc*s + (h-hc)*C;

Andre Staalstrom

maga O teu comentario está á espera de ser moderado. Marzo 11, 2011 ás 23:54 | Responder

Hi,
I understand what parameter I need to use, but Could you tell me please which depth is “h”¿¿??

Because, from roms configuration I get the following:

N = 32; vertical levels
theta_s = 6.; theta surface control parameter
theta_b = 0.; theta b bottom control parameter
hc = 0.9;

[nacholibre]

Could you tell me please which depth is “h”¿¿??

That is the depth at each grid cell. It is saved in the grid file and all history files with the name "h".

[wilkin]

That code is not correct for sea level zeta not equal to zero.

It also is not applicable to coordinate stretching options other than the old default ocean-s-coordinate.

The coordinate stretching is explained in detail here:
https://www.myroms.org/wiki/index.php/V ... coordinate

There are directions to accurate and up-to-date Matlab code (Utility/scoord.m) for computing the z coordinate at this post:
viewtopic.php?f=13&t=1240&p=4334&hilit= ... oord#p4334