River discharge issues

[gmaze]

Hey all,

I am having some problems getting river discharge in my model. I have been looking at the threads but I have not found anything that has fixed my problem.
It is compiling without issue and running without blowing up.

See below...

Any help would be appreciated....


Grace

I have defined to following in my options:


#define T_PASSIVE
#define TS_PSOURCE
#define UV_PSOURCE


in ana_psource I have the following:
!
!-----------------------------------------------------------------------
! Set tracer and/or mass point sources and/or sink.
!-----------------------------------------------------------------------
!
IF (iic(ng).eq.ntstart(ng)) THEN
!
! Set-up point Sources/Sink number (Nsrc), direction (Dsrc), I- and
! J-grid locations (Isrc,Jsrc). Currently, the direction can be along
! XI-direction (Dsrc = 0) or along ETA-direction (Dsrc > 0). The
! mass sources are located at U- or V-points so the grid locations
! should range from 1 =< Isrc =< L and 1 =< Jsrc =< M.

#elif defined MINE
Nsrc=2
Dsrc(1)=0.0_r8
Isrc(1)=75
Jsrc(1)=64
Dsrc(2)=0.0_r8
Isrc(2)=94
Jsrc(2)=69

!
! Set-up vertically integrated mass transport (m3/s) of point
! Sources/Sinks (positive in the positive U- or V-direction and
! viceversa).
# if defined MINE
IF ((tdays(ng)-dstart).lt.0.5_r8) THEN
fac=1.0_r8+TANH((time(ng)-43200.0_r8)/43200.0_r8)
ELSE
fac=-1.0_r8
END IF
DO is=1,Nsrc
Qbar(is)=fac*1500.0_r8
END DO
!
! Set-up tracer (tracer units) point Sources/Sinks.
# if defined MINE
IF (DOMAIN(ng)%NorthEast_Test(tile)) THEN
DO k=1,N(ng)
DO is=1,Nsrc
Tsrc(is,k,itemp)=0.0_r8
Tsrc(is,k,isalt)=0.0_r8
Tsrc(is,k,idye)=0.0_r8
END DO
END DO
END IF





./oceanG<ocean_mine.in

Model Input Parameters: ROMS/TOMS version 3.6
Tuesday - October 23, 2012 - 1:40:39 PM
-----------------------------------------------------------------------------

Wind-Driven Upwelling/Downwelling over a Periodic Channel Testing KB

Operating system : Darwin
CPU/hardware : x86_64
Compiler system : gfortran
Compiler command : /Applications/bin/gfortran
Compiler flags : -frepack-arrays -g -fbounds-check -ffree-form -ffree-line-length-none -ffree-form -ffree-line-length-none -ffree-form -ffree-line-length-none

SVN Root URL : https://www.myroms.org/svn/src/trunk
SVN Revision : exported

Local Root : /Users/gracemaze/ROMS/trunk
Header Dir : /Users/gracemaze/ROMS/Projects/Mine
Header file : mine.h
Analytical Dir: /Users/gracemaze/ROMS/Projects/Mine

Resolution, Grid 01: 0098x0098x020, Parallel Threads: 1, Tiling: 001x001


Physical Parameters, Grid: 01
=============================

2880 ntimes Number of timesteps for 3-D equations.
300.000 dt Timestep size (s) for 3-D equations.
30 ndtfast Number of timesteps for 2-D equations between
each 3D timestep.
1 ERstr Starting ensemble/perturbation run number.
1 ERend Ending ensemble/perturbation run number.
0 nrrec Number of restart records to read from disk.
T LcycleRST Switch to recycle time-records in restart file.
288 nRST Number of timesteps between the writing of data
into restart fields.
1 ninfo Number of timesteps between print of information
to standard output.
T ldefout Switch to create a new output NetCDF file(s).
72 nHIS Number of timesteps between the writing fields
into history file.
1 ntsAVG Starting timestep for the accumulation of output
time-averaged data.
72 nAVG Number of timesteps between the writing of
time-averaged data into averages file.
1 ntsDIA Starting timestep for the accumulation of output
time-averaged diagnostics data.
72 nDIA Number of timesteps between the writing of
time-averaged data into diagnostics file.
0.0000E+00 nl_tnu2(01) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 01: temp
0.0000E+00 nl_tnu2(02) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 02: salt
0.0000E+00 nl_tnu2(03) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 03: dye_01
5.0000E+00 nl_visc2 NLM Horizontal, harmonic mixing coefficient
(m2/s) for momentum.
1.0000E-06 Akt_bak(01) Background vertical mixing coefficient (m2/s)
for tracer 01: temp
1.0000E-06 Akt_bak(02) Background vertical mixing coefficient (m2/s)
for tracer 02: salt
1.0000E-06 Akt_bak(03) Background vertical mixing coefficient (m2/s)
for tracer 03: dye_01
1.0000E-05 Akv_bak Background vertical mixing coefficient (m2/s)
for momentum.
3.0000E-04 rdrg Linear bottom drag coefficient (m/s).
3.0000E-03 rdrg2 Quadratic bottom drag coefficient.
2.0000E-02 Zob Bottom roughness (m).
1.0000E-01 Dcrit Minimum depth for wetting and drying (m).
2 Vtransform S-coordinate transformation equation.
4 Vstretching S-coordinate stretching function.
0.0000E+00 theta_s S-coordinate surface control parameter.
0.0000E+00 theta_b S-coordinate bottom control parameter.
1.0000E+16 Tcline S-coordinate surface/bottom layer width (m) used
in vertical coordinate stretching.
1025.000 rho0 Mean density (kg/m3) for Boussinesq approximation.
0.000 dstart Time-stamp assigned to model initialization (days).
0.00 time_ref Reference time for units attribute (yyyymmdd.dd)
0.0000E+00 Tnudg(01) Nudging/relaxation time scale (days)
for tracer 01: temp
0.0000E+00 Tnudg(02) Nudging/relaxation time scale (days)
for tracer 02: salt
0.0000E+00 Tnudg(03) Nudging/relaxation time scale (days)
for tracer 03: dye_01
0.0000E+00 Znudg Nudging/relaxation time scale (days)
for free-surface.
0.0000E+00 M2nudg Nudging/relaxation time scale (days)
for 2D momentum.
0.0000E+00 M3nudg Nudging/relaxation time scale (days)
for 3D momentum.
0.0000E+00 obcfac Factor between passive and active
open boundary conditions.
T VolCons(1) NLM western edge boundary volume conservation.
T VolCons(2) NLM southern edge boundary volume conservation.
F VolCons(3) NLM eastern edge boundary volume conservation.
T VolCons(4) NLM northern edge boundary volume conservation.
14.000 T0 Background potential temperature (C) constant.
35.000 S0 Background salinity (PSU) constant.
1027.000 R0 Background density (kg/m3) used in linear Equation
of State.
1.7000E-04 Tcoef Thermal expansion coefficient (1/Celsius).
0.0000E+00 Scoef Saline contraction coefficient (1/PSU).
1.000 gamma2 Slipperiness variable: free-slip (1.0) or
no-slip (-1.0).
T LtracerSrc(01) Processing point sources/Sink on tracer 01: temp
T LtracerSrc(02) Processing point sources/Sink on tracer 02: salt
T LtracerSrc(03) Processing point sources/Sink on tracer 03: dye_01
T Hout(idFsur) Write out free-surface.
T Hout(idUbar) Write out 2D U-momentum component.
T Hout(idVbar) Write out 2D V-momentum component.
T Hout(idUvel) Write out 3D U-momentum component.
T Hout(idVvel) Write out 3D V-momentum component.
T Hout(idWvel) Write out W-momentum component.
T Hout(idOvel) Write out omega vertical velocity.
T Hout(idTvar) Write out tracer 01: temp
T Hout(idTvar) Write out tracer 02: salt
T Hout(idUsms) Write out surface U-momentum stress.
T Hout(idVsms) Write out surface V-momentum stress.
T Hout(inert) Write out inert passive tracer 01: dye_01

T Aout(idFsur) Write out averaged free-surface.
T Aout(idUbar) Write out averaged 2D U-momentum component.
T Aout(idVbar) Write out averaged 2D V-momentum component.
T Aout(idUvel) Write out averaged 3D U-momentum component.
T Aout(idVvel) Write out averaged 3D V-momentum component.
T Aout(idWvel) Write out averaged W-momentum component.
T Aout(idOvel) Write out averaged omega vertical velocity.
T Aout(idTvar) Write out averaged tracer 01: temp
T Aout(idTvar) Write out averaged tracer 02: salt
T Aout(idUsms) Write out averaged surface U-momentum stress.
T Aout(idVsms) Write out averaged surface V-momentum stress.
T Aout(inert) Write out averaged inert passive tracer 01: dye_01

T Dout(M2rate) Write out 2D momentun acceleration.
T Dout(M2pgrd) Write out 2D momentum pressure gradient.
T Dout(M2fcor) Write out 2D momentum Coriolis force.
T Dout(M2hadv) Write out 2D momentum horizontal advection.
T Dout(M2xadv) Write out 2D momentum horizontal X-advection.
T Dout(M2yadv) Write out 2D momentum horizontal Y-advection.
T Dout(M2hvis) Write out 2D momentum horizontal viscosity.
T Dout(M2xvis) Write out 2D momentum horizontal X-viscosity.
T Dout(M2yvis) Write out 2D momentum horizontal Y-viscosity.
T Dout(M2sstr) Write out 2D momentum surface stress.
T Dout(M2bstr) Write out 2D momentum bottom stress.

T Dout(M3rate) Write out 3D momentun acceleration.
T Dout(M3pgrd) Write out 3D momentum pressure gradient.
T Dout(M3fcor) Write out 3D momentum Coriolis force.
T Dout(M3hadv) Write out 3D momentum horizontal advection.
T Dout(M3xadv) Write out 3D momentum horizontal X-advection.
T Dout(M3yadv) Write out 3D momentum horizontal Y-advection.
T Dout(M3vadv) Write out 3D momentum vertical advection.
T Dout(M3hvis) Write out 3D momentum horizontal viscosity.
T Dout(M3xvis) Write out 3D momentum horizontal X-viscosity.
T Dout(M3yvis) Write out 3D momentum horizontal Y-viscosity.
T Dout(M3vvis) Write out 3D momentum vertical viscosity.

T Dout(iTrate) Write out rate of change of tracer 01: temp
T Dout(iTrate) Write out rate of change of tracer 02: salt
T Dout(iTrate) Write out rate of change of tracer 03: dye_01
T Dout(iThadv) Write out horizontal advection, tracer 01: temp
T Dout(iThadv) Write out horizontal advection, tracer 02: salt
T Dout(iThadv) Write out horizontal advection, tracer 03: dye_01
T Dout(iTxadv) Write out horizontal X-advection, tracer 01: temp
T Dout(iTxadv) Write out horizontal X-advection, tracer 02: salt
T Dout(iTxadv) Write out horizontal X-advection, tracer 03: dye_01
T Dout(iTyadv) Write out horizontal Y-advection, tracer 01: temp
T Dout(iTyadv) Write out horizontal Y-advection, tracer 02: salt
T Dout(iTyadv) Write out horizontal Y-advection, tracer 03: dye_01
T Dout(iTvadv) Write out vertical advection, tracer 01: temp
T Dout(iTvadv) Write out vertical advection, tracer 02: salt
T Dout(iTvadv) Write out vertical advection, tracer 03: dye_01
T Dout(iThdif) Write out horizontal diffusion, tracer 01: temp
T Dout(iThdif) Write out horizontal diffusion, tracer 02: salt
T Dout(iThdif) Write out horizontal diffusion, tracer 03: dye_01
T Dout(iTxdif) Write out horizontal X-diffusion, tracer 01: temp
T Dout(iTxdif) Write out horizontal X-diffusion, tracer 02: salt
T Dout(iTxdif) Write out horizontal X-diffusion, tracer 03: dye_01
T Dout(iTydif) Write out horizontal Y-diffusion , tracer 01: temp
T Dout(iTydif) Write out horizontal Y-diffusion , tracer 02: salt
T Dout(iTydif) Write out horizontal Y-diffusion, tracer 03: dye_01
T Dout(iTvdif) Write out vertical diffusion, tracer 01: temp
T Dout(iTvdif) Write out vertical diffusion, tracer 02: salt
T Dout(iTvdif) Write out vertical diffusion, tracer 03: dye_01

Output/Input Files:

Output Restart File: ocean_rst.nc
Output History File: ocean_his.nc
Output Averages File: ocean_avg.nc
Output Diagnostics File: ocean_dia.nc
Input Grid File: wfsbathl2.nc

Tile partition information for Grid 01: 0098x0098x0020 tiling: 001x001

tile Istr Iend Jstr Jend Npts

0 1 98 1 98 192080

Tile minimum and maximum fractional grid coordinates:
(interior points only)

tile Xmin Xmax Ymin Ymax grid

0 0.50 99.50 0.50 99.50 RHO-points

0 0.00 99.00 0.50 99.50 U-points

0 0.50 99.50 0.00 99.00 V-points

Lateral Boundary Conditions: NLM
============================

Variable Grid West Edge South Edge East Edge North Edge
--------- ---- ---------- ---------- ---------- ----------

zeta 1 Radiation Radiation Closed Radiation

ubar 1 Radiation Radiation Closed Radiation

vbar 1 Radiation Radiation Closed Radiation

u 1 Radiation Radiation Closed Radiation

v 1 Radiation Radiation Closed Radiation

temp 1 Radiation Radiation Closed Radiation

salt 1 Radiation Radiation Closed Radiation

dye_01 1 Radiation Radiation Closed Radiation

Activated C-preprocessing Options:

MINE Wind-Driven Upwelling/Downwelling over a Periodic Channel Testing KB
ANA_BPFLUX Analytical bottom passive tracers fluxes.
ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ANA_INITIAL Analytical initial conditions.
ANA_PASSIVE Analytical initial conditions for inert tracers.
ANA_PSOURCE Analytical point sources and sinks.
ANA_SMFLUX Analytical kinematic surface momentum flux.
ANA_SPFLUX Analytical surface passive tracer fluxes.
ANA_SSFLUX Analytical kinematic surface salinity flux.
ANA_STFLUX Analytical kinematic surface temperature flux.
ANA_VMIX Analytical vertical mixing coefficients.
ASSUMED_SHAPE Using assumed-shape arrays.
AVERAGES Writing out time-averaged nonlinear model fields.
DIAGNOSTICS_TS Computing and writing tracer diagnostic terms.
DIAGNOSTICS_UV Computing and writing momentum diagnostic terms.
DJ_GRADPS Parabolic Splines density Jacobian (Shchepetkin, 2002).
DOUBLE_PRECISION Double precision arithmetic.
MASKING Land/Sea masking.
MIX_S_TS Mixing of tracers along constant S-surfaces.
MIX_S_UV Mixing of momentum along constant S-surfaces.
NONLINEAR Nonlinear Model.
!NONLIN_EOS Linear Equation of State for seawater.
POWER_LAW Power-law shape time-averaging barotropic filter.
PROFILE Time profiling activated .
!RST_SINGLE Double precision fields in restart NetCDF file.
SALINITY Using salinity.
SOLVE3D Solving 3D Primitive Equations.
SPLINES Conservative parabolic spline reconstruction.
T_PASSIVE Advecting and diffusing inert passive tracer.
TS_U3HADVECTION Third-order upstream horizontal advection of tracers.
TS_C4VADVECTION Fourth-order centered vertical advection of tracers.
TS_DIF2 Harmonic mixing of tracers.
TS_PSOURCE Tracers point sources and sinks.
UV_ADV Advection of momentum.
UV_COR Coriolis term.
UV_U3HADVECTION Third-order upstream horizontal advection of 3D momentum.
UV_C4VADVECTION Fourth-order centered vertical advection of momentum.
UV_LDRAG Linear bottom stress.
UV_PSOURCE Mass point sources and sinks.
UV_VIS2 Harmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.
WET_DRY Wetting and drying activated.

Process Information:

Thread # 0 (pid= 78192) is active.

INITIAL: Configuring and initializing forward nonlinear model ...


Vertical S-coordinate System:

level S-coord Cs-curve Z at hmin at hc half way at hmax

20 0.0000000 0.0000000 0.000 0.000 0.000 0.000
19 -0.0500000 -0.0025000 -0.000 -2.714 -1.357 -2.714
18 -0.1000000 -0.0100000 -0.000 -5.428 -2.714 -5.428
17 -0.1500000 -0.0225000 -0.000 -8.142 -4.071 -8.142
16 -0.2000000 -0.0400000 -0.000 -10.856 -5.428 -10.856
15 -0.2500000 -0.0625000 -0.000 -13.570 -6.785 -13.570
14 -0.3000000 -0.0900000 -0.000 -16.284 -8.142 -16.284
13 -0.3500000 -0.1225000 -0.000 -18.998 -9.499 -18.998
12 -0.4000000 -0.1600000 -0.000 -21.713 -10.856 -21.713
11 -0.4500000 -0.2025000 -0.000 -24.427 -12.213 -24.427
10 -0.5000000 -0.2500000 -0.000 -27.141 -13.570 -27.141
9 -0.5500000 -0.3025000 -0.000 -29.855 -14.927 -29.855
8 -0.6000000 -0.3600000 -0.000 -32.569 -16.284 -32.569
7 -0.6500000 -0.4225000 -0.000 -35.283 -17.641 -35.283
6 -0.7000000 -0.4900000 -0.000 -37.997 -18.998 -37.997
5 -0.7500000 -0.5625000 -0.000 -40.711 -20.356 -40.711
4 -0.8000000 -0.6400000 -0.000 -43.425 -21.713 -43.425
3 -0.8500000 -0.7225000 -0.000 -46.139 -23.070 -46.139
2 -0.9000000 -0.8100000 -0.000 -48.853 -24.427 -48.853
1 -0.9500000 -0.9025000 -0.000 -51.567 -25.784 -51.567
0 -1.0000000 -1.0000000 -0.000 -54.281 -27.141 -54.281

Time Splitting Weights: ndtfast = 30 nfast = 42

Primary Secondary Accumulated to Current Step

1-0.0008094437383769 0.0333333333333333-0.0008094437383769 0.0333333333333333
2-0.0014053566728197 0.0333603147912792-0.0022148004111966 0.0666936481246126
3-0.0017877524645903 0.0334071600137066-0.0040025528757869 0.1001008081383191
4-0.0019566842408176 0.0334667517625262-0.0059592371166046 0.1335675599008453
5-0.0019122901320372 0.0335319745705535-0.0078715272486418 0.1670995344713988
6-0.0016548570247459 0.0335957175749547-0.0095263842733877 0.2006952520463536
7-0.0011849025289723 0.0336508794757796-0.0107112868023600 0.2343461315221331
8-0.0005032751608631 0.0336903762267453-0.0112145619632232 0.2680365077488784
9 0.0003887272597151 0.0337071520654408-0.0108258347035081 0.3017436598143192
10 0.0014892209965583 0.0336941944901169-0.0093366137069498 0.3354378543044362
11 0.0027955815694920 0.0336445537902317-0.0065410321374578 0.3690824080946679
12 0.0043042707117221 0.0335513677379153-0.0022367614257356 0.4026337758325831
13 0.0060106451121704 0.0334078920475245 0.0037738836864348 0.4360416678801076
14 0.0079087469427945 0.0332075372104522 0.0116826306292293 0.4692492050905598
15 0.0099910761708920 0.0329439123123590 0.0216737068001213 0.5021931174029188
16 0.0122483446563884 0.0326108764399960 0.0339220514565097 0.5348039938429148
17 0.0146692120341107 0.0322025982847830 0.0485912634906204 0.5670065921276978
18 0.0172400033810439 0.0317136245503127 0.0658312668716643 0.5987202166780105
19 0.0199444086685725 0.0311389577709445 0.0857756755402368 0.6298591744489550
20 0.0227631639997064 0.0304741441486588 0.1085388395399432 0.6603333185976138
21 0.0256737146312911 0.0297153720153352 0.1342125541712342 0.6900486906129490
22 0.0286498597812016 0.0288595815276255 0.1628624139524359 0.7189082721405746
23 0.0316613792205220 0.0279045862015855 0.1945237931729578 0.7468128583421600
24 0.0346736416507075 0.0268492068942347 0.2291974348236653 0.7736620652363948
25 0.0376471948657328 0.0256934188392112 0.2668446296893981 0.7993554840756060
26 0.0405373376992232 0.0244385123436867 0.3073819673886213 0.8237939964192927
27 0.0432936737565710 0.0230872677537126 0.3506756411451924 0.8468812641730054
28 0.0458596469320356 0.0216441452951603 0.3965352880772280 0.8685254094681656
29 0.0481720587108284 0.0201154903974257 0.4447073467880565 0.8886408998655914
30 0.0501605672561820 0.0185097551070648 0.4948679140442384 0.9071506549726561
31 0.0517471682814030 0.0168377361985254 0.5466150823256415 0.9239883911711815
32 0.0528456577069106 0.0151128305891453 0.5994607400325521 0.9391012217603267
33 0.0533610761022577 0.0133513086655816 0.6528218161348098 0.9524525304259084
34 0.0531891349131379 0.0115726061288397 0.7060109510479478 0.9640251365547481
35 0.0522156244733761 0.0097996349650684 0.7582265755213239 0.9738247715198165
36 0.0503158038019030 0.0080591141492892 0.8085423793232269 0.9818838856691057
37 0.0473537721847153 0.0063819206892258 0.8558961515079423 0.9882658063583315
38 0.0431818225418188 0.0048034616164019 0.8990779740497611 0.9930692679747334
39 0.0376397765791564 0.0033640675316746 0.9367177506289175 0.9964333355064080
40 0.0305543017255206 0.0021094083123694 0.9672720523544381 0.9985427438187774
41 0.0217382098544504 0.0010909315881854 0.9890102622088885 0.9996336754069628
42 0.0109897377911118 0.0003663245930371 1.0000000000000004 0.9999999999999999

ndtfast, nfast = 30 42 nfast/ndtfast = 1.40000

Centers of gravity and integrals (values must be 1, 1, approx 1/2, 1, 1):

1.000000000000 1.047601458608 0.523800729304 1.000000000000 1.000000000000

Power filter parameters, Fgamma, gamma = 0.28400 0.18933

Minimum X-grid spacing, DXmin = 1.36680654E+00 km
Maximum X-grid spacing, DXmax = 2.25429594E+00 km
Minimum Y-grid spacing, DYmin = 1.12125285E+00 km
Maximum Y-grid spacing, DYmax = 1.74711789E+00 km
Minimum Z-grid spacing, DZmin = 5.00000000E-16 m
Maximum Z-grid spacing, DZmax = 2.71406742E+00 m

Minimum barotropic Courant Number = 2.48902798E-09
Maximum barotropic Courant Number = 2.25749765E-01
Maximum Coriolis Courant Number = 2.05965957E-02


Maximum grid stiffness ratios: rx0 = 1.000000E+00 (Beckmann and Haidvogel)
rx1 = 3.900000E+01 (Haney)


Initial basin volumes: TotVolume = 5.4829748446E+11 m3
MinVolume = 1.4393857369E-09 m3
MaxVolume = 8.6887313418E+06 m3
Max/Min = 6.0364161734E+15

[nacholibre]

You must be defining them on a masked cell if you don't see the discharge comming out.
use mask_u or mask_v in Matlab
the wet cell at the boundary is the source location e.g. mask_u(439,16) for a river in u direction
in ROMS the indices will be one less and swapped e.g. river_Xposition=15, river_Eposition=438

[gmaze]

nacholibre,

I have checked a few times and the river positions are not on masked land. I am still not seeing the discharge no matter where I place the rivers. Any other suggestions?

[gmaze]

Also, I would like to use a netCDF forcing file instead of the ana_psource file. I have read the posts here but honestly I have no idea how to make it. Is there a step by step guide?

[kate]

If your confusion is about NetCDF proper, I put together a talk about it some time ago which you can find here.

There is an example of the cdl format for the sources file under Data/ROMS/CDL in the ROMS source.

I'm sure someone has an Matlab script for creating such a file. My favorite was written in Perl Data Language (PDL), but that's a bit obscure.

[hetland]

Here is a (slightly old) python script that will give you some guidance in making river forcing files.

Also, if you don't see the rivers at all, you need to back up some steps. Make sure you can get the RIVERPLUME test case to work, and then work from there. Make sure you are compiling the code you think you are in the f90 files.

Code: Select all

#!/usr/bin/env python
# encoding: utf-8
"""
untitled.py
Created on 2008-04-24

Copyright (c) 2008 Robert D. Hetland
"""
__docformat__ = "restructuredtext en"

import os
import sys
from datetime import datetime

import numpy as np
import matplotlib.pyplot as plt
import netCDF4
import netcdftime

import octant

time_units = 'days since 1970-01-01 00:00:00'
utime = netcdftime.utime(time_units)

river_data = np.load('tarbert_landing.npy')
river_jd = utime.date2num(river_data['date'])

river = np.arange(1, 56)

river_Xposition = \
    np.array([  86.,   86.,   87.,   87.,   87.,   87.,   87.,   86.,   86.,
             87.,   88.,   89.,   90.,   91.,   92.,   93.,   94.,   95.,
             96.,   97.,   98.,  100.,  100.,  100.,  100.,  100.,   98.,
             98.,   98.,   98.,   98.,   98.,   99.,   99.,   99.,   99.,
             99.,   97.,   97.,   97.,   82.,   82.,   82.,   31.,   32.,
             33.,   34.,   35.,   36.,   37.,   38.,   39.,   40.,   41.,   42.])

river_Eposition = \
    np.array([ 28.,  27.,  26.,  25.,  24.,  23.,  22.,  21.,  19.,  19.,  19.,
            19.,  18.,  13.,  16.,  16.,  16.,  16.,  15.,  16.,  16.,  17.,
            18.,  19.,  20.,  21.,  22.,  23.,  24.,  25.,  26.,  27.,  28.,
            29.,  30.,  31.,  32.,  33.,  34.,  35.,  18.,  19.,  20.,  63.,
            63.,  63.,  63.,  63.,  63.,  63.,  63.,  63.,  63.,  63.,  63.])

river_direction = \
    np.array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,
            1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,
            0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
            0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
            1.,  1.,  1.])

# multiplication factor for each point source.  Multiply by Mississippi River
# discharge to get the right point source flux.
river_factor = \
    np.array([-0.01624998, -0.01624998, -0.01624998, -0.01624998, -0.01624998,
           -0.01624998, -0.01624998, -0.01624998, -0.01624998, -0.01624998,
           -0.01624998, -0.01624998, -0.01624998, -0.01624998, -0.01624998,
           -0.01624998, -0.01624998, -0.01624998, -0.01624998, -0.01624998,
           -0.01624998,  0.01624998,  0.01624998,  0.01624998,  0.01624998,
            0.01624998,  0.01624998,  0.01624998,  0.01624998,  0.01624998,
            0.01624998,  0.01624998,  0.01624998,  0.01624998,  0.01624998,
            0.01624998,  0.01624998,  0.01624998,  0.01624998,  0.01624998,
           -0.11666688, -0.11666688, -0.11666688, -0.03449878, -0.03812705,
           -0.04130254, -0.04385654, -0.04564636, -0.04656848, -0.04656848,
           -0.04564636, -0.04385654, -0.04130254, -0.03812705, -0.03449878])

river_flag = 5

# Set river temp equal to local air temp
nc = netCDF4.Dataset('../dasilva_data/ss_air_temp.cdf')
air_temp = nc.variables['sat']
air_temp.set_auto_maskandscale(True)
air_temp = air_temp[:,118,270]
nc.close()

sat = []
for yr in range(1980, 2010):
    for mo in range(1, 13):
        jd = utime.date2num(datetime(yr, mo, 15))
        sat.append((jd, air_temp[mo-1]))

sat = np.asarray(sat, dtype=[('jd', 'd'), ('air_temp', 'd')])
river_temp = np.interp(river_jd, sat['jd'], sat['air_temp'])
river_temp = river_temp[:,np.newaxis, np.newaxis] * np.ones((20, len(river)))

river_salt = 0.0

river_O2 = octant.ocean.o2_sat(river_temp, river_salt)

river_Miss = np.zeros_like(river_temp)
river_Miss[:,:,river_Xposition > 60] = 1.0
river_Atch = np.zeros_like(river_temp)
river_Atch[:,:,river_Xposition < 60] = 1.0

vshape_default = \
    np.array([ 0.0047619 ,  0.00952381,  0.01428571,  0.01904762,  0.02380952,
            0.02857143,  0.03333333,  0.03809524,  0.04285714,  0.04761905,
            0.05238095,  0.05714286,  0.06190476,  0.06666667,  0.07142857,
            0.07619048,  0.08095238,  0.08571429,  0.09047619,  0.0952381 ,
            ])[:,np.newaxis]

vshape_swpass = \
    np.array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
            0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
            0.        ,  0.        ,  0.02777778,  0.05555556,  0.08333333,
            0.11111111,  0.13888889,  0.16666667,  0.19444444,  0.22222222,
            ])[:,np.newaxis]


river_Vshape = vshape_default * np.ones( (1, len(river)), 'd')
river_Vshape[:,40:43] = vshape_swpass

river_transport = river_data['q'][:, np.newaxis] * river_factor[np.newaxis, :]


###################################
# Crate netCDF file

# dimensions:
#   xi_rho = 128 ;
#   eta_rho = 64 ;
#   s_rho = 20 ;
#   river = 55 ;
#   river_time = 3652 ;
# variables:
#   double river(river) ;
#   double river_Xposition(river) ;
#   double river_Eposition(river) ;
#   double river_direction(river) ;
#   double river_Vshape(s_rho, river) ;
#   double river_time(river_time) ;
#   double river_transport(river_time, river) ;
#   double river_temp(river_time, s_rho, river) ;
#   double river_salt(river_time, s_rho, river) ;
#   double river_flag(river) ;
#   double river_O2(river_time, s_rho, river) ;
#   double river_Miss(river_time, s_rho, river) ;
#   double river_Atch(river_time, s_rho, river) ;
# 
# // global attributes:
#       :title = "MRP river forcing file" ;
#       :type = "FORCING file" ;
#       :author = "Rob Hetland" ;
#       :date = "20-Nov-2003 18:35:37" ;
# }

nc = netCDF4.Dataset('mch_river_frc.nc', 'w', format='NETCDF3_CLASSIC')

nc.createDimension('x_rho', 128)
nc.createDimension('eta_rho', 64)
nc.createDimension('s_rho', 20)
nc.createDimension('river', len(river))
nc.createDimension('river_time', len(river_jd))

def write_nc_var(name, dimensions, var, units=None):
    nc.createVariable(name, 'f8', dimensions)
    if units is not None:
        nc.variables[name].units = units
    nc.variables[name][:] = var

write_nc_var('river', ('river', ), river, 'River index')
write_nc_var('river_Xposition', ('river', ), river_Xposition, 'River xi index')
write_nc_var('river_Eposition', ('river', ), river_Eposition, 'River eta index')
write_nc_var('river_direction', ('river', ), river_direction, 'River direction')
write_nc_var('river_flag', ('river', ), river_flag, 'River flag')

write_nc_var('river_Vshape', ('s_rho', 'river'), river_Vshape, 'River vertical shape')
write_nc_var('river_time', ('river_time', ), river_jd, time_units)
write_nc_var('river_transport', ('river_time', 'river'), river_transport, 
             'River transport [m3/s]')
write_nc_var('river_temp', ('river_time', 's_rho', 'river'), 
              river_temp,
             'River temperature [deg C]')
write_nc_var('river_salt', ('river_time', 's_rho', 'river'), river_salt,
             'River salinity [psu]')

# Tracers
write_nc_var('river_O2', ('river_time', 's_rho', 'river'), river_O2,
            'River oxygen [mmol/l]')
write_nc_var('river_Miss', ('river_time', 's_rho', 'river'), river_Miss,
            'Mississippi river tag [none]')
write_nc_var('river_Atch', ('river_time', 's_rho', 'river'), river_Atch,
            'Atchafalaya river tag [none]')

nc.close()

[gmaze]

Kate and Hetland,

Thank you for your help. I have gone back and run the river plume case and it is working properly. I also looked at the information you gave Kate and thank you it has helped. I do still have some questions.

I am not sure how to specify the data in the forcing file. I have the following but I am not sure if it is correct, do I need to have river_Vshape? And should I have a transport value for each time and river?

Code: Select all

netcdf frc_rivers {

dimensions:
	xi_rho = 100 ;
	xi_u = 99 ;
	xi_v = 100 ;
	eta_rho = 100 ;
	eta_u = 100;
	eta_v = 99 ;
	s_rho = 20 ;
	river = 2 ;
	time = UNLIMITED ; // (0 currently)

variables:
	double river(river) ;
		river:long_name = "river runoff identification number" ;
		river:units = "nondimensional" ;
		river:field = "river, scalar" ;
	double river_Xposition(river) ;
		river_Xposition:long_name = "river XI-position at RHO-points" ;
		river_Xposition:units = "nondimensional" ;
		river_Xposition:valid_min = 1. ;
		river_Xposition:valid_max = 99. ;
		river_Xposition:field = "river_Xposition, scalar" ;
	double river_Eposition(river) ;
		river_Eposition:long_name = "river ETA-position at RHO-points" ;
		river_Eposition:units = "nondimensional" ;
		river_Eposition:valid_min = 1. ;
		river_Eposition:valid_max = 99. ;
		river_Eposition:field = "river_Eposition, scalar" ;
	double river_direction(river) ;
		river_direction:long_name = "river runoff direction" ;
		river_direction:units = "nondimensional" ;
		river_direction:field = "river_direction, scalar" ;
	double river_flag(river) ;
		river_flag:long_name = "river runoff tracer flag" ;
		river_flag:option_0 = "all tracers are off" ;
		river_flag:option_1 = "only temperature is on" ;
		river_flag:option_2 = "only salinity is on" ;
		river_flag:option_3 = "both temperature and salinity are on" ;
		river_flag:units = "nondimensional" ;
		river_flag:field = "river_flag, scalar" ;
	double river_Vshape(s_rho, river) ;
		river_Vshape:long_name = "river runoff mass transport vertical profile" ;
		river_Vshape:units = "nondimensional" ;
		river_Vshape:field = "river_Vshape, scalar" ;
	double river_time(time) ;
		river_time:long_name = "river runoff time" ;
		river_time:units = "days since 1992-01-01 00:00:00" ;
		river_time:add_offset = 2448623. ;
		river_time:field = "river_time, scalar, series" ;
	double river_transport(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" ;
		river_transport:time = "river_time" ;
	double river_temp(time, s_rho, river) ;
		river_temp:long_name = "river runoff potential temperature" ;
		river_temp:units = "Celsius" ;
		river_temp:field = "river_temp, scalar, series" ;
		river_temp:time = "river_time" ;
	double river_salt(time, s_rho, river) ;
		river_salt:long_name = "river runoff salinity" ;
		river_salt:units = "PSU" ;
		river_salt:field = "river_salt, scalar, series" ;
		river_salt:time = "river_time" ;
	double river_dye1(time, s_rho, river) ;
		river_dye1:long_name = "river runoff dye" ;
		river_dye1:units = "PSU" ;
		river_dye1:field = "river_dye1, scalar, series" ;
		river_dye1:time = "river_time" ;

// global attributes:
		:type = "ROMS FORCING file" ;
		:title = "River Forcing" ;
		:grd_file = "wfsbathl2.nc" ;
		:rivers = "(1)Peace River (2) Caloosahatchee River" ;

data:
	river=1, 2;
	river_Xposition=64, 69;
	river_Eposition=75, 94;
	river_direction=1, 0;
	river_flag=3, 3;
	river_transport=50 , -100;
	river_temp=0, 0;
	river_salt=0, 0;
	river_dye1 =0, 0;
	
}

[hetland]

I am not sure how to specify the data in the forcing file. I have the following but I am not sure if it is correct, do I need to have river_Vshape? And should I have a transport value for each time and river?
[/code]

Yes, and yes. I am not sure if there is a 'default' Vshape, but I think it is best to specify explicitly what you are doing even if there is such a thing. And you do need a transport value for each flux point (sometimes rivers have more than one grid cell source) and at each time you have data -- say once per day. ROMS will lineally interpolate to the model time, just as with other forcing files.

[gmaze]

Hi all,
I have the rivers working now, however the model will not finish running, it giver the following error:


Blowing-up: Saving latest model state into RESTART file

WRT_RST - wrote re-start fields (Index=1,2) into time record = 0000003

I have checked the restart file and the locations that I have my rivers give a large salinity value (>50) even though I am specifying a discharge of 20 salinity. The temperature and other tracer I am using are not having an issue. I am not sure how to fix this problem.

[tony1230]

Apparently,model told that you have placed the rivers at points that not should be them. Refer to viewtopic.php?f=14&t=2200&hilit=river. Hope to make your work done.

- Shou

[fagundesmo]

Dear Hetland,

I am trying to use your python script to make my river. Therefore, I would like to know how your tarbert_landing.npy looks like. Because I have stuck in the line where it asks for the river_jd.

I have pasted below a piece of how my data. May you please give me some help?


Thank you in advance.

# Year month P1 P2 temp
1980 01 079.756 132.386 26.1
1980 02 129.282 199.633 25.7
1980 03 295.348 952.241 25.6
1980 04 152.666 408.861 25.8
1980 05 080.884 234.951 25.9
1980 06 063.296 133.935 25.9
1980 07 056.141 061.283 25.6
1980 08 051.546 033.302 25.9
1980 09 049.332 026.841 26.3
1980 10 094.267 118.854 26.6
1980 11 056.810 048.177 26.9
1980 12 055.793 038.445 26.7
1981 01 069.653 072.073 26.1
1981 02 108.795 171.042 25.7
1981 03 079.136 104.761 25.6
1981 04 170.014 156.083 25.8
1981 05 065.415 073.774 25.9
1981 06 052.819 037.025 25.9
1981 07 047.817 022.991 25.6
1981 08 044.433 020.674 25.9
1981 09 068.814 021.423 26.3
1981 10 044.024 017.311 26.6
1981 11 047.681 036.347 26.9
1981 12 054.634 018.101 26.7