U and Ubar patterns disagreement in the results

[feroda]

I check my model output by using of Ncview. What's a really big surprise to me, the Ubar and U patterns looks totally different, which should exactly be. Are you guys surprised upon seeing the .ps file bellow?

I have checked my ocean-grd.nc carefully, the mask_r,mask_psi,mask_u,mask_v looks fine. I am doubt about why the 3-D U mask pattern in the model output does Not conform the mask_u in the grid file.

I am greatly appreciated your precious comment and valuable advice on it. Thank you very much!

[kate]

Ew, there's an indexing problem somewhere. What do you get on an ncdump of the file you are plotting? U and ubar should have the same horizontal indices. Which version of ROMS is this? Are you the first to use last week's update?

As for mask_u, it should follow directly from mask_rho. I think ROMS recomputes it to be consistent with the mask_rho you pass in and therefore doesn't need to read it from the grid file at all.

[feroda]

What do you get on an ncdump of the file you are plotting?

netcdf ocean_avg_0001 {
dimensions:
xi_rho = 362 ;
xi_u = 361 ;
xi_v = 362 ;
xi_psi = 361 ;
eta_rho = 322 ;
eta_u = 322 ;
eta_v = 321 ;
eta_psi = 321 ;
s_rho = 30 ;
s_w = 31 ;
tracer = 2 ;
boundary = 4 ;
ocean_time = UNLIMITED ; // (0 currently)
variables:
int ntimes ;
ntimes:long_name = "number of long time-steps" ;
int ndtfast ;
ndtfast:long_name = "number of short time-steps" ;
double dt ;
dt:long_name = "size of long time-steps" ;
dt:units = "second" ;
double dtfast ;
dtfast:long_name = "size of short time-steps" ;
dtfast:units = "second" ;
double dstart ;
dstart:long_name = "time stamp assigned to model initilization" ;
dstart:units = "days since 0001-01-01 00:00:00" ;
dstart:calendar = "julian" ;
int nHIS ;
nHIS:long_name = "number of time-steps between history records" ;
int ndefHIS ;
ndefHIS:long_name = "number of time-steps between the creation of history files" ;
int nRST ;
nRST:long_name = "number of time-steps between restart records" ;
nRST:cycle = "only latest two records are maintained" ;
int ntsAVG ;
ntsAVG:long_name = "starting time-step for accumulation of time-averaged fields" ;
int nAVG ;
nAVG:long_name = "number of time-steps between time-averaged records" ;
int ndefAVG ;
ndefAVG:long_name = "number of time-steps between the creation of average files" ;
int ntsDIA ;
ntsDIA:long_name = "starting time-step for accumulation of diagnostic fields" ;
int nDIA ;
nDIA:long_name = "number of time-steps between diagnostic records" ;
int ndefDIA ;
ndefDIA:long_name = "number of time-steps between the creation of diagnostic files" ;
double Falpha ;
Falpha:long_name = "Power-law shape barotropic filter parameter" ;
double Fbeta ;
Fbeta:long_name = "Power-law shape barotropic filter parameter" ;
double Fgamma ;
Fgamma:long_name = "Power-law shape barotropic filter parameter" ;
double tnu2(tracer) ;
tnu2:long_name = "Laplacian mixing coefficient for tracers" ;
tnu2:units = "meter2 second-1" ;
double visc2 ;
visc2:long_name = "Laplacian mixing coefficient for momentum" ;
visc2:units = "meter2 second-1" ;
double Akt_bak(tracer) ;
Akt_bak:long_name = "background vertical mixing coefficient for tracers" ;
Akt_bak:units = "meter2 second-1" ;
double Akv_bak ;
Akv_bak:long_name = "background vertical mixing coefficient for momentum" ;
Akv_bak:units = "meter2 second-1" ;
double rdrg ;
rdrg:long_name = "linear drag coefficient" ;
rdrg:units = "meter second-1" ;
double rdrg2 ;
rdrg2:long_name = "quadratic drag coefficient" ;
double Zob ;
Zob:long_name = "bottom roughness" ;
Zob:units = "meter" ;
double Zos ;
Zos:long_name = "surface roughness" ;
Zos:units = "meter" ;
double Znudg ;
Znudg:long_name = "free-surface nudging/relaxation inverse time scale" ;
Znudg:units = "day-1" ;
double M2nudg ;
M2nudg:long_name = "2D momentum nudging/relaxation inverse time scale" ;
M2nudg:units = "day-1" ;
double M3nudg ;
M3nudg:long_name = "3D momentum nudging/relaxation inverse time scale" ;
M3nudg:units = "day-1" ;
double Tnudg(tracer) ;
Tnudg:long_name = "Tracers nudging/relaxation inverse time scale" ;
Tnudg:units = "day-1" ;
double rho0 ;
rho0:long_name = "mean density used in Boussinesq approximation" ;
rho0:units = "kilogram meter-3" ;
double R0 ;
R0:long_name = "background density used in linear equation of state" ;
R0:units = "kilogram meter-3" ;
double Tcoef ;
Tcoef:long_name = "thermal expansion coefficient" ;
Tcoef:units = "Celsius-1" ;
double Scoef ;
Scoef:long_name = "Saline contraction coefficient" ;
Scoef:units = "PSU-1" ;
double gamma2 ;
gamma2:long_name = "slipperiness parameter" ;
char spherical ;
spherical:long_name = "grid type logical switch" ;
spherical:option_T = "spherical" ;
spherical:option_F = "Cartesian" ;
double xl ;
xl:long_name = "domain length in the XI-direction" ;
xl:units = "meter" ;
double el ;
el:long_name = "domain length in the ETA-direction" ;
el:units = "meter" ;
double theta_s ;
theta_s:long_name = "S-coordinate surface control parameter" ;
double theta_b ;
theta_b:long_name = "S-coordinate bottom control parameter" ;
double Tcline ;
Tcline:long_name = "S-coordinate surface/bottom layer width" ;
Tcline:units = "meter" ;
double hc ;
hc:long_name = "S-coordinate parameter, critical depth" ;
hc:units = "meter" ;
double s_rho(s_rho) ;
s_rho:long_name = "S-coordinate at RHO-points" ;
s_rho:valid_min = -1. ;
s_rho:valid_max = 0. ;
s_rho:standard_name = "ocean_s_coordinate" ;
s_rho:formula_terms = "s: s_rho eta: zeta depth: h a: theta_s b: theta_b depth_c: hc" ;
s_rho:field = "s_rho, scalar" ;
double s_w(s_w) ;
s_w:long_name = "S-coordinate at W-points" ;
s_w:valid_min = -1. ;
s_w:valid_max = 0. ;
s_w:standard_name = "ocean_s_coordinate" ;
s_w:formula_terms = "s: s_w eta: zeta depth: h a: theta_s b: theta_b depth_c: hc" ;
s_w:field = "s_w, scalar" ;
double Cs_r(s_rho) ;
Cs_r:long_name = "S-coordinate stretching curves at RHO-points" ;
Cs_r:valid_min = -1. ;
Cs_r:valid_max = 0. ;
Cs_r:field = "Cs_r, scalar" ;
double Cs_w(s_w) ;
Cs_w:long_name = "S-coordinate stretching curves at W-points" ;
Cs_w:valid_min = -1. ;
Cs_w:valid_max = 0. ;
Cs_w:field = "Cs_w, scalar" ;
double h(eta_rho, xi_rho) ;
h:long_name = "bathymetry at RHO-points" ;
h:units = "meter" ;
h:coordinates = "lon_rho lat_rho" ;
h:field = "bath, scalar" ;
double f(eta_rho, xi_rho) ;
f:long_name = "Coriolis parameter at RHO-points" ;
f:units = "second-1" ;
f:coordinates = "lon_rho lat_rho" ;
f:field = "coriolis, scalar" ;
double pm(eta_rho, xi_rho) ;
pm:long_name = "curvilinear coordinate metric in XI" ;
pm:units = "meter-1" ;
pm:coordinates = "lon_rho lat_rho" ;
pm:field = "pm, scalar" ;
double pn(eta_rho, xi_rho) ;
pn:long_name = "curvilinear coordinate metric in ETA" ;
pn:units = "meter-1" ;
pn:coordinates = "lon_rho lat_rho" ;
pn:field = "pn, scalar" ;
double 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" ;
double 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" ;
double lon_u(eta_u, xi_u) ;
lon_u:long_name = "longitude of U-points" ;
lon_u:units = "degree_east" ;
lon_u:field = "lon_u, scalar" ;
double lat_u(eta_u, xi_u) ;
lat_u:long_name = "latitude of U-points" ;
lat_u:units = "degree_north" ;
lat_u:field = "lat_u, scalar" ;
double lon_v(eta_v, xi_v) ;
lon_v:long_name = "longitude of V-points" ;
lon_v:units = "degree_east" ;
lon_v:field = "lon_v, scalar" ;
double lat_v(eta_v, xi_v) ;
lat_v:long_name = "latitude of V-points" ;
lat_v:units = "degree_north" ;
lat_v:field = "lat_v, scalar" ;
double lon_psi(eta_psi, xi_psi) ;
lon_psi:long_name = "longitude of PSI-points" ;
lon_psi:units = "degree_east" ;
lon_psi:field = "lon_psi, scalar" ;
double lat_psi(eta_psi, xi_psi) ;
lat_psi:long_name = "latitude of PSI-points" ;
lat_psi:units = "degree_north" ;
lat_psi:field = "lat_psi, scalar" ;
double angle(eta_rho, xi_rho) ;
angle:long_name = "angle between XI-axis and EAST" ;
angle:units = "radians" ;
angle:coordinates = "lon_rho lat_rho" ;
angle:field = "angle, scalar" ;
double mask_rho(eta_rho, xi_rho) ;
mask_rho:long_name = "mask on RHO-points" ;
mask_rho:option_0 = "land" ;
mask_rho:option_1 = "water" ;
mask_rho:coordinates = "lon_rho lat_rho" ;
double mask_u(eta_u, xi_u) ;
mask_u:long_name = "mask on U-points" ;
mask_u:option_0 = "land" ;
mask_u:option_1 = "water" ;
mask_u:coordinates = "lon_u lat_u" ;
double mask_v(eta_v, xi_v) ;
mask_v:long_name = "mask on V-points" ;
mask_v:option_0 = "land" ;
mask_v:option_1 = "water" ;
mask_v:coordinates = "lon_v lat_v" ;
double mask_psi(eta_psi, xi_psi) ;
mask_psi:long_name = "mask on psi-points" ;
mask_psi:option_0 = "land" ;
mask_psi:option_1 = "water" ;
mask_psi:coordinates = "lon_psi lat_psi" ;
double ocean_time(ocean_time) ;
ocean_time:long_name = "averaged time since initialization" ;
ocean_time:units = "seconds since 0001-01-01 00:00:00" ;
ocean_time:calendar = "julian" ;
ocean_time:field = "time, scalar, series" ;
float zeta(ocean_time, eta_rho, xi_rho) ;
zeta:long_name = "time-averaged free-surface" ;
zeta:units = "meter" ;
zeta:time = "ocean_time" ;
zeta:coordinates = "lon_rho lat_rho ocean_time" ;
zeta:field = "free-surface, scalar, series" ;
zeta:_FillValue = 1.e+37f ;
float ubar(ocean_time, eta_u, xi_u) ;
ubar:long_name = "time-averaged vertically integrated u-momentum component" ;
ubar:units = "meter second-1" ;
ubar:time = "ocean_time" ;
ubar:coordinates = "lon_u lat_u ocean_time" ;
ubar:field = "ubar-velocity, scalar, series" ;
ubar:_FillValue = 1.e+37f ;
float vbar(ocean_time, eta_v, xi_v) ;
vbar:long_name = "time-averaged vertically integrated v-momentum component" ;
vbar:units = "meter second-1" ;
vbar:time = "ocean_time" ;
vbar:coordinates = "lon_v lat_v ocean_time" ;
vbar:field = "vbar-velocity, scalar, series" ;
vbar:_FillValue = 1.e+37f ;
float u(ocean_time, s_rho, eta_u, xi_u) ;
u:long_name = "time-averaged u-momentum component" ;
u:units = "meter second-1" ;
u:time = "ocean_time" ;
u:coordinates = "lon_u lat_u s_rho ocean_time" ;
u:field = "u-velocity, scalar, series" ;
u:_FillValue = 1.e+37f ;
float v(ocean_time, s_rho, eta_v, xi_v) ;
v:long_name = "time-averaged v-momentum component" ;
v:units = "meter second-1" ;
v:time = "ocean_time" ;
v:coordinates = "lon_v lat_v s_rho ocean_time" ;
v:field = "v-velocity, scalar, series" ;
v:_FillValue = 1.e+37f ;
float omega(ocean_time, s_w, eta_rho, xi_rho) ;
omega:long_name = "time-averaged S-coordinate vertical momentum component" ;
omega:units = "meter3 second-1" ;
omega:time = "ocean_time" ;
omega:coordinates = "lon_rho lat_rho s_w ocean_time" ;
omega:field = "omega, scalar, series" ;
omega:_FillValue = 1.e+37f ;
float w(ocean_time, s_w, eta_rho, xi_rho) ;
w:long_name = "time-averaged vertical momentum component" ;
w:units = "meter second-1" ;
w:time = "ocean_time" ;
w:coordinates = "lon_rho lat_rho s_w ocean_time" ;
w:field = "w-velocity, scalar, series" ;
w:_FillValue = 1.e+37f ;
float temp(ocean_time, s_rho, eta_rho, xi_rho) ;
temp:long_name = "time-averaged potential temperature" ;
temp:units = "Celsius" ;
temp:time = "ocean_time" ;
temp:coordinates = "lon_rho lat_rho s_rho ocean_time" ;
temp:field = "temperature, scalar, series" ;
temp:_FillValue = 1.e+37f ;
float salt(ocean_time, s_rho, eta_rho, xi_rho) ;
salt:long_name = "time-averaged salinity" ;
salt:time = "ocean_time" ;
salt:coordinates = "lon_rho lat_rho s_rho ocean_time" ;
salt:field = "salinity, scalar, series" ;
salt:_FillValue = 1.e+37f ;

// global attributes:
:file = "ocean_avg_0001.nc" ;
:type = "ROMS/TOMS averages file" ;
:Conventions = "CF-1.0" ;
:title = "Wind-Driven Upwelling/Downwelling over a Periodic Channel" ;
:var_info = "./varinfo.dat" ;
:rst_file = "ocean_rst.nc" ;
:his_base = "ocean_his" ;
:avg_base = "ocean_avg" ;
:dia_file = "ocean_dia.nc" ;
:grd_file = "./Input/scs-grd.nc" ;
:script_file = "ocean_upwelling.in" ;
:svn_url = "https://www.myroms.org/svn/src/trunk" ;
:svn_rev = "291M" ;
:code_dir = "/home/dkwang/trunk" ;
:header_dir = "/home/dkwang/SCS" ;
:header_file = "upwelling.h" ;
:os = "Linux" ;
:cpu = "x86_64" ;
:compiler_system = "ifort" ;
:compiler_command = "/usr/local/mpich2/bin/mpif90" ;
:compiler_flags = " -ip -O3 -xW -free" ;
:tiling = "001x010" ;
:history = "ROMS/TOMS, Version 3.2, Monday - January 19, 2009 - 6:01:00 PM" ;
:ana_file = "ROMS/Functionals/ana_btflux.h, ROMS/Functionals/ana_initial.h, ROMS/Functionals/ana_smflux.h, ROMS/Functionals/ana_stflux.h, ROMS/Functionals/ana_vmix.h" ;
:CPP_options = "UPWELLING, ANA_BSFLUX, ANA_BTFLUX, ANA_INITIAL, ANA_SMFLUX, ANA_SSFLUX, ANA_STFLUX, ANA_VMIX, ASSUMED_SHAPE, AVERAGES, CURVGRID, DIAGNOSTICS_TS, DIAGNOSTICS_UV, DJ_GRADPS, DOUBLE_PRECISION, EW_PERIODIC, INLINE_2DIO, MASKING, MIX_S_TS, MIX_S_UV, MPI, NONLINEAR, !NONLIN_EOS, POWER_LAW, PROFILE, !RST_SINGLE, SALINITY, SOLVE3D, SPLINES, TS_U3HADVECTION, TS_C4VADVECTION, TS_DIF2, UV_ADV, UV_COR, UV_U3HADVECTION, UV_C4VADVECTION, UV_LDRAG, UV_VIS2, VAR_RHO_2D" ;
}

[feroda]

Which version of ROMS is this? Are you the first to use last week's update?

The version is
svn: $LastChangedRevision: 291 $
svn: $LastChangedDate: 2008-12-23 22:42:42 -0500 (Tue, 23 Dec 2008) $

So, it is not the latest version.

The method to calculate mask_u is from seagrd2roms.m in seagrid package.
Specifically:

umask = zeros(size(rmask));

% The xi direction (left-right):
LP = (m-1)/2; % The rho dimension.
% The eta direction (up-down):
MP = (n-1)/2; % The rho dimension.
for i = 2:LP
for j = 1:MP
umask(j, i-1) = rmask(j, i) * rmask(j, i-1);
end
end
nc{'mask_u'}(:) = umask(1:end, 1:end-1);

Hi Kate, I learn from your response that ROMS calculates mask_u from mask_rho input from the grid file. But, the mask_rho in my grid file looks fine. So, ROMS should get mask_u by itself, shouldn't it?

Would you please explain more about the "indexing problem"? How could I make sure where the error really is and solve it?

Thanks!

[kate]

I don't see anything wrong there. Ncview works for me in terms of indexing the u variables. What version of that do you have? I've got 1.93c, probably not the latest. Do you have any other way to visualize the model fields?

[feroda]

I don't see anything wrong there. Ncview works for me in terms of indexing the u variables. What version of that do you have? I've got 1.93c, probably not the latest. Do you have any other way to visualize the model fields?

Thank you very much!
The version of NCview for me is:
Ncview 1.93f David W. Pierce 12 September 2008

Yep, I visualize the model fields in matlab as well as Plotting Package(based on NCL) recommended by wikiroms. The results are with the same disagreement. So, I think it is not the ncview's problem.
I met this problem two weeks ago. I am trying my best to figure out the problem. Unfortunately, all ideas do NOT work.

[feroda]

Another great doubt for me is that, the model can run smoothly even though the the mask_u seems totally wrong and the model results for other variables are greatly reasonable. Why? Thank you!

[kate]

Your fields likely are fine inside the model. It's not that the mask for u is funky, it's that somehow it got shifted somewhere in the writing or plotting so that line 2 is off by 1, line 3 is off by 2, and so on. It's like u got put into an array that's the correct dimension for rho. I don't recall this being a recent ROMS error.

Ideas:

1. Download ROMS again and see if it changes.
2. Try a test problem or similar configuration with one box island - is there such a test problem? It looks like FLT_TEST, RIVERPLUME1 and RIVERPLUME2 use ana_mask.h. If you can reproduce it there, it would be easier for others to investigate.

[m.hadfield]

Hi Kate, I learn from your response that ROMS calculates mask_u from mask_rho input from the grid file. But, the mask_rho in my grid file looks fine. So, ROMS should get mask_u by itself, shouldn't it?

I don't think this is true. I don't see code to calculate umask anywhere except ana_mask.h. I recall that a few years ago I passed inconsistent u, v and rho masks to ROMS and the model blew up soon after initialisation.

[feroda]

Your fields likely are fine inside the model. It's not that the mask for u is funky, it's that somehow it got shifted somewhere in the writing or plotting so that line 2 is off by 1, line 3 is off by 2, and so on. It's like u got put into an array that's the correct dimension for rho. I don't recall this being a recent ROMS error.

Ideas:

1. Download ROMS again and see if it changes.

Hi Kate, I am trying the latest ROMS. There are some trouble in compile the codes.
The errors are:

/home/dkwang/Projects/SCS/Bath-ok/Build/libMODS.a -L/usr/local/netcdf-3.6.1/lib -lnetcdf
IPO link: Warning unknown option '-h'.
IPO link: can not find "eap-arrays"
ifort: error: problem during multi-file optimization compilation (code 1)
make: *** [/home/Projects/SCS/Bath-ok/oceanM] Error 1

Does it mean netcdf-4 and hdf5 is needed by the latest ROMS?
Thank you!

[kate]

No, it has nothing to do with Netcdf4 or HDF5. You somehow have a -heap-arrays argument to the loader. It's likely coming from a file in the Compilers directory. I see you are using ifort: can you fix the corresponding file?

[feroda]

No, it has nothing to do with Netcdf4 or HDF5. You somehow have a -heap-arrays argument to the loader. It's likely coming from a file in the Compilers directory. I see you are using ifort: can you fix the corresponding file?

Thank you, Kate. Yep, I use ifort to compile ROMS.
Do you mean
"FFLAGS := -heap-arrays"
in the Linux-ifort.mk should be commented out?

[feroda]

I try the latest version of ROMS to check the problem stated as the title.
Even worse problem appears:

STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME

0 0 00:00:00 0.000000E+00 3.967632E+03 3.967632E+03 5.390653E+14
DEF_HIS - creating history file: ocean_his_0001.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
1 0 00:01:30 NaN NaN NaN NaN

Elapsed CPU time (seconds):

It means that the mask on u, v ,psi and rho are NOT coordinate with each other.
I check the mask_rho, mask_u and mask_v by using of Ncview, they looks well. But
mask_psi meets the same pattern disagreement problem!

So, it come to the conclusion that:
My model setup
goes well in ROMS with the version of 177;
goes well in ROMS with the version of 291, but get wrong u-mask in the model output which does NOT affect the running of the model;
goes bad in ROMS with the version of 304(latest) while reading the mask_psi from the grid input file.

Where is the way to go ahead?

Thank you very much!

[m.hadfield]

Your skewed u (but not ubar) field looks very like something I encountered a week ago. I filed a bug report here

https://www.myroms.org/projects/src/ticket/268

If you can't see that page, here is what it says

When DISTRIBUTE and INLINE_2DIO are defined (which normally happens only on the UNICOS_SN platform) there is a misalignment problem when 3D u-grid variables are written to a netCDF file: successive rows of data are out of step by 1, leading to a skewing of the fields. The same error occurs (I think) when data are read from a restart file.

I have traced this to functions mp_gather2d and mp_scatter2d in distribute.F. These have code to adjust the offset for U and V grids:

IF ((MyType.eq.p2dvar).or.(MyType.eq.u2dvar)) Io=1
IF ((MyType.eq.p2dvar).or.(MyType.eq.v2dvar)) Jo=1

However this doesn't anticipate the possibility that mp_gather2d and mp_scatter2d will be called for 3D variables, as happens when INLINE_2DIO is defined. To cover this case, the code should be changed to the following

IF ((MyType.eq.p2dvar).or.(MyType.eq.u2dvar).or. &
& (MyType.eq.p3dvar).or.(MyType.eq.u3dvar)) Io=1
IF ((MyType.eq.p2dvar).or.(MyType.eq.v2dvar).or. &
& (MyType.eq.p3dvar).or.(MyType.eq.v3dvar)) Jo=1

A modified version of distribute.F is attached.

Now this applies only when DISTRIBUTE and INLINE_2DIO are defined, and by default this occurs only on the UNICOS_SN platform (ie. Cray T3E, pretty rare now). So it probably isn't directly relevant in your case (though it is strikingly similar!). However, given the recent major changes to the I/O system you might try either updating to the latest version (304) or reverting to a version before the I/O changes.

RE the NaNs in the initial model fields, I also saw this recently, when I constructed an initialisation file from a history or averages file. As of a couple of months ago, history and averages file have NaNs in masked areas for easier visualisation. This has had the side effect of makeing them unsuitable for model initialisation: you need to replace the NaNs with finite values.

[feroda]

I just want to let ROMSers know that the latest version 310 is ok for my application.
It means the issue described in the topic has been settled by Hernan.
Thanks for his great job.

Enjoy