How can make sure proper river point position in ROMS

[CBian]

Actually this is the same problem I have post in ROMS Bugs: Minus sediment concentration in river sediment case.
But I think this is not a sediment concentration problem any more. So I make a new post so that more people can have a look.

These days I tried to figure out how the river point was added to the ROMS. There are 10 rivers in my model domain. I tried to put all the rivers into the model. After lots of tests, I added 7 of them to the model successfully. However, other 3 rivers cannot be accepted near their location. These 3 rivers always have abnormal salt value (salt>100), or disappear.

I read a message from kate:

If a source goes east to west, its indices will be at a place with zero mask_u, since that is the land-sea boundary. In fact, any source should be where the mask_u or mask_v is zero, with land on one side and water on the other.

So I test the mask_rho, mask_u, mask_v in the river point.
The results are very interesting, all the rivers can be successfully added in the model have follow characters:
mask_rho(lat,lon-1)==0 & mask_rho(lat,lon)==0 & mask_rho(lat,lon+1)==1
mask_u(lat,lon-1)==0 & mask_u(lat,lon)==0 & mask_u(lat,lon+1)==1
mask_v(lat,lon-1)==0 & mask_v(lat,lon)==0 & mask_v(lat,lon+1)==1

I output all the points in my domain which have these characters and test some of them. All of them can be added as a river successfully. And if only two of mask_rho, mask_u, mask_v comply with up characters, the river also cannot be added.

But these case is applied for the river flow from west to east. If I have a river flow east to west, there will no point can have characters like:
mask_rho(lat,lon-1)==1 & mask_rho(lat,lon)==0 & mask_rho(lat,lon+1)==0
mask_u(lat,lon-1)==1 & mask_u(lat,lon)==0 & mask_u(lat,lon+1)==0
mask_v(lat,lon-1)==1 & mask_v(lat,lon)==0 & mask_v(lat,lon+1)==0

Because the mask_u are calculated like this:
[Mp,Lp]=size(mask_rho);
M=Mp-1;
L=Lp-1;
mask_u=mask_rho(:,1:L).*mask_rho(:,2:Lp);
so if mask_rho(lat,lon-1)==1 & mask_rho(lat,lon)==0, mask_u(lat,lon-1)must be 0 instead of 1.

I also tried lots of points, I cannot add river which flows from east to west.
And there are two rivers in model domain which flow from north to south, I also canon add them in model.

So, is there anyone who can add a river flow from east to west and the river did not disappear and the salt value is also ok ?? Or the river can be added into model without the characters mentioned above?

Any suggestion will be appreciated!

[mathieu]

The river inflow in ROMS is quite complicated. In the netcdf file of rivers you have the following arrays:

Code: Select all

        double river_Xposition(river) ;
                river_Xposition:long_name = "river ETA-position at RHO-points" ;
                river_Xposition:units = "nondimensional" ;
                river_Xposition:field = "river_Xposition, scalar" ;
        double river_Eposition(river) ;
                river_Eposition:long_name = "river XI-position at RHO-points" ;
                river_Eposition:units = "nondimensional" ;
                river_Eposition:field = "river_Xposition, scalar" ;
        double river_direction(river) ;
                river_direction:long_name = "river runoff direction" ;
                river_direction:units = "nondimensional" ;
                river_direction:field = "river_direction, scalar" ;
        double river_transport(river_time, river) ;
                river_transport:long_name = "river runoff vertically integrated mass transport" ;
                river_transport:units = "meter3 second-1" ;
                river_transport:field = "river_transport, scalar, series" ;

The river_Xposition and river_Eposition gives the position of the river. The river_direction gives the direction of the river, i.e. east-west or south-north. The point is that whether it is going north or south, west or east is determined by the river_transport itself.
This makes for complicated formulas. I advise you to use map_rivers.m which allows you to draw the river from the grid file and the river file.
Now for the formula you have 4 cases.
Let us write eEta, eXi, eDir and eTrans for the Eta point, Xi point, direction and transport in the river file. Let us call iEtaSea, iXiSea and iEtaLand, iXiLand for the points of the grid where the river arrives and where it is coming from (we set 1<=iEtaSea<=eta_rho and similarly for other variables)
Then:

• If eDir=0 and eTrans>0 then iEtaSea=eEta+1, iXiSea=eXi+1, iEtaLand=eEta+1, iXiLand=eXi (river goes east)
  If eDir=1 and eTrans>0 then iEtaSea=eEta+1, iXiSea=eXi+1, iEtaLand=eEta, iXiLand=eXi+1 (river goes north)
  If eDir=0 and eTrans<0 then iEtaSea=eEta+1, iXiSea=eXi, iEtaLand=eEta+1, iXiLand=eXi+1 (river goes west)
  If eDir=1 and eTrans<0 then iEtaSea=eEta, iXiSea=eXi+1, iEtaLand=eEta+1, iXiLand=eXi+1 (river goes south)

Naturally enough, all this are interpretation. There is no such thing in ROMS as testing if the river goes north or south, only inflow/outflow matters and this is set by the sign. But if your river has a constant direction then the above formula should help you.

[CBian]

Thank you very much, mathieu, this is really helpful!

I never know the river_transport could be negative. When I change the sign of the river_transport, I can add all the rivers in model domain.

I cannot find map_rivers.m in the web, could you attach this code?

[austinctodd]

You can easily see if you have the correct cell by taking all rivers that you wish to flow east-west and plot the river_Xposition and river_Eposition on a map of the land mask at u-points, and similarly on the v mask for all north-south points.

An eastward flowing river should be a flux through the land mask, just to the left of the land/sea boundary.

[winterbellgirl]

The river inflow in ROMS is quite complicated. In the netcdf file of rivers you have the following arrays:

Code: Select all

        double river_Xposition(river) ;
                river_Xposition:long_name = "river ETA-position at RHO-points" ;
                river_Xposition:units = "nondimensional" ;
                river_Xposition:field = "river_Xposition, scalar" ;
        double river_Eposition(river) ;
                river_Eposition:long_name = "river XI-position at RHO-points" ;
                river_Eposition:units = "nondimensional" ;
                river_Eposition:field = "river_Xposition, scalar" ;
        double river_direction(river) ;
                river_direction:long_name = "river runoff direction" ;
                river_direction:units = "nondimensional" ;
                river_direction:field = "river_direction, scalar" ;
        double river_transport(river_time, river) ;
                river_transport:long_name = "river runoff vertically integrated mass transport" ;
                river_transport:units = "meter3 second-1" ;
                river_transport:field = "river_transport, scalar, series" ;

The river_Xposition and river_Eposition gives the position of the river. The river_direction gives the direction of the river, i.e. east-west or south-north. The point is that whether it is going north or south, west or east is determined by the river_transport itself.
This makes for complicated formulas. I advise you to use map_rivers.m which allows you to draw the river from the grid file and the river file.
Now for the formula you have 4 cases.
Let us write eEta, eXi, eDir and eTrans for the Eta point, Xi point, direction and transport in the river file. Let us call iEtaSea, iXiSea and iEtaLand, iXiLand for the points of the grid where the river arrives and where it is coming from (we set 1<=iEtaSea<=eta_rho and similarly for other variables)
Then:

• If eDir=0 and eTrans>0 then iEtaSea=eEta+1, iXiSea=eXi+1, iEtaLand=eEta+1, iXiLand=eXi (river goes east)
  If eDir=1 and eTrans>0 then iEtaSea=eEta+1, iXiSea=eXi+1, iEtaLand=eEta, iXiLand=eXi+1 (river goes north)
  If eDir=0 and eTrans<0 then iEtaSea=eEta+1, iXiSea=eXi, iEtaLand=eEta+1, iXiLand=eXi+1 (river goes west)
  If eDir=1 and eTrans<0 then iEtaSea=eEta, iXiSea=eXi+1, iEtaLand=eEta+1, iXiLand=eXi+1 (river goes south)

Naturally enough, all this are interpretation. There is no such thing in ROMS as testing if the river goes north or south, only inflow/outflow matters and this is set by the sign. But if your river has a constant direction then the above formula should help you.

Hi, mathieu, could you please give some explanations on the formula for river going in the 4 directions? I don't quite understand why it sets the location of the river point in that way. In your case, are all the coordinates the ones for the rho points? Why don't we put the river point at the rho point of the land/sea grid? Thanks!

[mathieu]

I have not derived the formula myself. I followed for that the map_rivers.m file which I put in attachment. This is not the original file but a hacked version of it. I do not know the original author. I think it is mandatory to use this program before any ROMS run using rivers.

The formulas are that way because the rivers are not put as from this rho point to this rho point but at this u- or v-point. If you look at step3d_uv.F then you see the following:

Code: Select all

u(i,j,k,nnew)=Qsrc(is,k)*cff1

and the same for v. So the flux is applied at one u- or v-point and this means that if the corresponding rho points are both wet then you move water from one cell to the next. If only one of the rho point is wet then you actually do river inflow (water goes to the sea) or outflow (water goes in the river, in tides for example). If both rho points are land then nothing happens. Please notice how land and sea are exchanged when the sign of eTrans changes.

The formula I gave are translations from the u- and v-points to rho-point setting but ROMS is an Arakawa C-grid model and it shows.

[winterbellgirl]

Thanks, mathiu, and I'm going to doublecheck it.

[tony1230]

Hi mathieu
Thank you for your .m script. But one question : what's the function for gmean

Code: Select all

lat=gmean(LAT_rho(:));

on the line 218 .
when i try to run the script , it collapsed.

tks anybody who reply me.

[tony1230]

ok i know . it's the geometric mean . Under my MATLAB version, the function is

geomean

!

[mathieu]

Nothing so complicated. it is simply the mean with the NaN values removed.

[tony1230]

Hi mathieu
I can run the .m script and the location-pic shown me where the river-sources are. But, something i still can not understand.
1. What's the two type of hint-points mean created by

Code: Select all

plot(LonDot,LatDot,'ro');   % arrow "tail"
plot(LonStar,LatStar,'*');    % arrow "head"

i see, if my point loc is (x,y), then the '*' point is(x+1,y+1) and 'o'point is(x,y+1).
so my question is which one of the 3 points is the actually one that i would put inflow on?
2.As you see hereinafter, 6 rivers in my domain.In my application the first 3 rivers are inflow from west and toward to east,the fourth one
is inflow south to north, the fifth and sixth inflow north and toward south.But after the run, nbRight=3 nbUp=3 nbLeft=0 nbDown=0.As i comprehend , nbRight=3 and nbUp=3 tells that 3 rivers run toward east and north, respectively.Am i right? If so , i can
not understand why i got the abnormal result!

nbRiver = 6
iRiver = 1 eRivDir = 0 eTransAvg = 640
LonStar=119.1726 LatStar=37.7466
iRiver = 2 eRivDir = 0 eTransAvg = 3.2
LonStar=118.9327 LatStar=37.3292
iRiver = 3 eRivDir = 0 eTransAvg = 28.5
LonStar=117.7776 LatStar=38.962
iRiver = 4 eRivDir = 1 eTransAvg = 6.4
LonStar=119.3275 LatStar=39.4383
iRiver = 5 eRivDir = 1 eTransAvg = 158.5
LonStar=121.8113 LatStar=40.9494
iRiver = 6 eRivDir = 1 eTransAvg = 57.1
LonStar=121.5663 LatStar=40.8809
nbRight=3 nbUp=3 nbLeft=0 nbDown=0

Any suggestion or reply will be greatly appreciated.

w.shou

[mathieu]

Some general remarks without looking at the details.
First, the points that are put in the river file are u- and v-points.
Second the coordinates are put with respect to the internal coordinate system of ROMS and there is an index shift documented by Hernan recently in the forum.
Lastly, the way rivers are going is indicated by the river_flux, it can change sign and so the river downstream
or upstream.

[xiaozhu557]

ok i know . it's the geometric mean . Under my MATLAB version, the function is

geomean

!

Hi Tony and all,
When I run the script map_rivers.m and changed the code 'gmean' to 'geomean', I have got an error that

Error using geomean(line 17)
X may not contain negative values.

ones the array LAT_rho(:) have negative values for Southern Hemisphere.

Anyone know how to fix it?

Thanks.

[rduran]

how bout sign(y_1)*geomean(abs(y_i))? (need y_i < 0 for all i)

[xiaozhu557]

Hi rduran,
Thanks for your reply. I have tried your suggestion, that line
lat=sign(LAT_rho(:)).*geomean(abs(LAT_rho(:)));
can work, but the other line
xfac=cos(lat*pi/180);
set(gca,'DataAspectRatio',[1 xfac 1]);
can not work, and I have got the error message with

Error using horzcat
CAT argument dimensions are not consistent.

Error in map_rivers (Line 220)
set(gca,'DataAspectRatio',[1 xfac 1]);

Then I try to use

Code: Select all

lat=sign(LAT_rho(:,:)).*geomean(abs(LAT_rho(:,:)));

or

Code: Select all

lat=sign(LAT_rho).*geomean(abs(LAT_rho));

or

Code: Select all

lat=sign(LAT_rho(:))*geomean(abs(LAT_rho));

Then I have got the error message with
Error using *(.*)
Inner matrix dimensions must agree.

I don't know the reason. But I would appreciate you very much still.
Anyway, I have tried to comment those three lines, and got the same result. So please don't take too much time on that problem any more.

Thanks.

[rduran]

Reminds me when I was starting out with matlab: you need to pay attention to the error messages. In this case the dimensions are not working out.

type

whos

to see the size of your variables.

type

help horzcat to see what it means.

Another option is copy-paste the error message into google, youd be surprised how much Dr. Google knows.

whatever you do you need to get your hands dirty and experiment with matlab until you learn what is going on and how to fix it.