Ocean Modeling Discussion

wilkin

Parker,

For the purposes of debugging your concerns about reproducibility, you might find it easier to repurpose a few passive tracers rather than endeavor to add to NBT by creating new variables with new entries in varinfo.yaml that require you to define a new index_code, metadata and I/O.

You ...

wilkin

You appear to be imposing boundary zeta and ubar independently as in-phase functions of sin(omega*t) with arbitrary amplitude.

These variables are dynamically linked in amplitude (through sqrt(gh) factor), and with their relative phase depending on whether the external tide is a progressive wave ...

wilkin

For LuvSrc the volume flux through the cell faces modifies FX and FE in pre_step3d.F. The tracer enters when, subsequently, in step3d_t.F we take the divergence of those fluxes. For LwSrc the divergence of tracer simply appears directly in step3d_t.F. Look in step3d_t.F for "LwSrc" to see where this ...

wilkin

I suggest you try setting h = 10 where mask_rho = 0 and #undef WET_DRY. You have very small (10^-4) values of h in the land that might cause trouble.

I notice that your land/sea mask is such that the Malaysia Peninsula connects to the southern boundary. The Singapore Strait is closed.

wilkin

I agree with Hernan that there might be a problem with your grid. The TOTAL VOLUME is reported as NaN at initialization, so there must be a cell where ∆x*∆y*(h+zeta) is NaN. You have WET_DRY defined, so this should have taken care of any h=0, and min(∆z) on water points is > 0.

But do you really ...

wilkin

Estimated maximum CFL values U∆t/∆x for u,v,w velocity are reported in the log as Cu,Cv,Cw. At the start of your run,
Cv = NaN so you have v velocity values that are NaN (in your initial conditions) which via continuity equation make w also NaN. You
cannot have NaN values in initial (or boundary ...

wilkin

I wrote that code a long time ago to process the rather peculiar way in which ECMWF were serving those data. Unpacking the time accumulated variables (precip and heat fluxes) was tricky, and is coded for only the specific download conditions of 3 hourly interval data 0 3 6 9 ... as noted in the ...

wilkin

Previously, taur due to rain was added as distinct vector components to taux,tauy near the end of the calculation. But in the new formulation we add the magnitude of the momentum of the rain to the magnitude of the wind stress, and then break it out into vector components in one step. This is ...

wilkin

We have observed some sensitivity to the compiler version using PIO.
On our cluster, ifort 19.1.1 works with ROMS PIO options ...

INP_LIB = 1
OUT_LIB = 2

! PIO library methods for reading/writing NetCDF files:
...
! [2] serial read and write of NetCDF3 (64-bit offset)
...
PIO_METHOD = 2 ...

wilkin

Free surface variability is driven by the divergence of the depth integrated velocity in the 2D engine. Attempting to nudge to a zeta field will be overcome by this unless a kinematically consistent adjustment is also make to ubar,vbar.

wilkin

Some NCEI Ocean NOMADS product (https://nomads.ncep.noaa.gov/) services are going to be decommissioned on May 15, 2025. The archived data will remain available, but no new data will be added.

Description of change:
New updates of Navy operational ocean model fields will be discontinued. This ...

wilkin

Check your f90 files in Build to see if you actually got what you intended.

wilkin

If the values you are taking from NEMO are, say, hourly resolution and show tidal variability, then if you use those as boundary conditions you not add tides.

However, if the NEMO outputs are averages over a day or more (so effectively "de-tided" by the time filter) then it makes sense to add ...

wilkin

If you look in set_tides.F you'll see how the 2-D harmonics are used in conjunction with the Flather open boundary condition ... If you have set the harmonic amplitudes to zero inside the domain, the averaging to u,v points will halve the values.

! If appropriate, load tidal forcing into boundary ...

wilkin

I can confirm that SSH_TIDES and UV_TIDES add harmonic tidal variability strictly at the model perimeter. That the information is provided via a 2-D field on the full model grid is a matter of legacy going way back. The values in the interior are unused for forcing. However, if the DETIDE option is ...

wilkin

I don't see why your inverse transform can't be an interpolation by table lookup. You have z(k) and k, so just do a linear interpolation on k and z (both are guaranteed monotonic). This will be robust for all ROMS output files because the ROMS ocean_s_coordinate_g1 etc. have CF-Convention ...

wilkin

Hello Erik,

I don't quite follow what you want the inverse relation for. Maybe you are rewriting the vertical advection for parcels in the native ROMS s-coordinate for better accuracy, as we did with the ROMSpath* offline particle tracking. But if that's the case you need ROMS Hz (layer thicknesses ...

wilkin

It is documented, with references, in ./Nonlinear/rho_eos.F

wilkin

There is xroms (Python) which will do this for you https://github.com/xoceanmodel/xroms

Also, the functions stretching.m and set_depth.m in the myroms matlab repo
https://github.com/myroms/roms_matlab

And ROMS itself will write z-coordinates to the output netcdf files if you activate the Hout ...

wilkin

Not interpolation method but extrapolation method (assuming you are using Matlab interp … you don't say whose code you are using)

wilkin

I suspect that in the vertical interpolation either the new surface or bottom level is outside the range of the old levels and extrapolation creates NANs. You might need to create a rule to catch this, or explicitly set the extrapolation method to 'nearest'

wilkin

GLS is fine. But I suggest you use the same advection scheme choice for both tracers. You have A4 for temperature and MPDATA for salinity which will introduce small differences that might be your noise.

wilkin

Mixing schemes for temp and salt might be ill-advised.

wilkin

So, in your roms.in you have ?

! The user has the option of specifying the full Keyword or the first
! two letters, regardless if using uppercase or lowercase. If nested
! grids, specify values for each grid (see glossary below).

Hadvection == U3 \ ! temperature
U3 ! salinity

Vadvection == C4 ...

wilkin

Re-reading this thread, I notice you don't say what horizontal and vertical advection schemes you are using, and the values of any explicit horizontal mixing coefficients for comparisons to your grid resolution. It may be that these are contributing to the noise and the vertical mixing is not the ...

wilkin

Many NCEI services hosted in Asheville, NC, are down until further notice due to the impacts of hurricane Helene.

Hurricane Helene made landfall Friday September 27th and proceeded to
move inland toward the Blue Ridge Mountain area and more specifically
Asheville, NC. Asheville, NC is home to NOAA ...

wilkin

At https://www.myroms.org/wiki/Numerical_Solution_Technique#Horizontal_Discretization scroll down to the small table that shows

Variable Interior Range Full Range

The valid index range for u-type points is i=1:L, and j=0:M.

You have specified u-face point sources (river_direction = 0) and i ...

wilkin

To debug this I recommend you keep your working test case UPWELLING more or less as is, and then successively add features (#define ...) from stormsurge.h to work toward the dynamics configuration you seek. Then you can zero in on the specific feature that is triggering the problem.

wilkin

If you literally have the CPP defs you posted earlier ...

/*************** Surface Forcing */
# define SOLAR_SOURCE
# define QCORRECTION
# define DIURNAL_SRFLUX
#LONGWAVE_OUT
then you do not have LONGWAVE_OUT defined, because there is no #define

Be careful with # define DIURNAL_SRFLUX. If the ...

wilkin

The only way to do this is to stop the model, add the passive dye to the last history file, and restart with the modified file as initial conditions.

wilkin

When using radiation data from ECMWF (or any meteorological model analysis) pay attention to the units.

If they are not Watts m -2 but Joules m -2 then they are accumulated over a time interval
https://codes.ecmwf.int/grib/param-db/147
which might be 3, 6 or even 9 up to 12 hours from the start of ...

wilkin

#define LONGWAVE_OUT is an option to compute outgoing (upward) longwave radition based on model SST. One would never specify upward longwave separately.
This post might be helpful: https://www.myroms.org/forum/viewtopic. ... OUT#p25453

wilkin

Note that there was a code update 6 weeks ago to introduce COARE3.5 bulk fluxes:
https://www.myroms.org/projects/src/ticket/970
The changes from 3.0 to 3.5 make the algorithm more appropriate at higher wind speeds.

If your issue is arising in shallow water you might need to limit the bottom stress ...

wilkin

The advective term in the xi-direction Huon_temp is H z times u over n times temp
where Hz is the layer thickness (so, dz) and over n is 1/n so dy.
Hence, the elemental cell face area dz*dy is already in the flux term. To vertically integrate over a region you just add up the terms over k.

In ...

wilkin

I thought the most correct way to get it right was to divide my calculated horizontal diffusion by Hz

Yes, that is correct. 1/m,1/n and Hz are the dx, dy and dz of the cell volume. Where you might lose some precision in your calculation is in the fact that Hz is time varying because it depends on ...

wilkin

Is there a way to double-check if my calculation is correct? Can it be compared somehow with the diagnostic diffusion term temp_hdiff?

Equation (10) at the link you gave to the Wiki https://www.myroms.org/wiki/Horizontal_Mixing#Horizontal_Diffusion shows that the net diffusion is the divergence ...

wilkin

ROMS should be able to find what it needs in the system defaults, so your sysadmin build of MPI might be problematic.

If you need to customize the paths, it's best to do it in ./Compilers/my_build_paths.csh which will be used when (in build_roms.sh) you activate ...

export USE_MY_LIBS=yes # use ...

wilkin

The effective drag coefficient for air-sea momentum exchange depends on many factors associated with marine boundary layers conditions. There is no "default value".

Read the papers listed in the opening lines of the subroutine.

! Fairall, C.W., E.F. Bradley, D.P. Rogers, J.B. Edson and G.S ...

wilkin

Sensible and latent heat are not input variables. They are only outputs that are calculated and optionally saved.

In the varinfo.yaml file look to the right of the variable name and you see an entry that is either #Input, #Output or #Input/Output .

sensible and latent are outputs only ...

wilkin

There's you answer:

Update your metadata file: /scratch/cdschai/ROMS_Working/nested_Bob_Exp/varinfo.dat

Your varinfo.dat -- really you should be updated to using varinfo.yaml -- is not consistent with your code.

wilkin

This is a COAWST specific problem, so you should pose your questions to that group.

That said, it looks like you are running an old version (you still have varinfo.dat, not varinfo.yaml) so it might be time to update in case you are encountering a solved problem.

The myroms.org repository has an ...

wilkin

I can restart with rst, but I need to undef PERFECT_RESTART and remove TWO dimensions,Or I use his file to restart.

No, removing dimensions from an existing netcdf file is tricky. Don't you have a HIS or AVG file you can use instead?

wilkin

The additional dimension in a RST file is there to enable a "perfect restart" as if the model never stopped. The long_name of the ru variable documents that it is the RHS of total u-momentum (and likewise for rv, rzeta), which is also required for a perfect restart of the time stepping.

You are ...

wilkin

Are you merging data from different files to create your restart? It looks like you have mismatched dimensions.

Restart RST files have an extra dimension to carry the n-1 and n time levels to restart the time stepping. But if you add in new biology variables that are not dimensioned the same way it ...

wilkin

Your nudging time scale is 365 days, so it won't do much on the 30-day time frame you are concerned about. Try making that time scale smaller, like 30 days.

You are nudging to an analytical climatology, so be sure that is correct.

wilkin

If you just want SST, use the "quicksave" file outputs and in roms.in set to TRUE the temp_sur option ...

Qout(idsurT) == T F ! temp_sur, salt_sur surface temperature and salinity

which is explained in the GLOSSARY at the end of every roms.in ...

wilkin

It doesn't make any physical/logical sense to have some variables periodic and others not. Consequently, the logical flags EWperiodic and NSperiodic are set to .TRUE. if any boundary is periodic. See Utility/inp_decode.F

If you look in Nonlinear/t3dbc_im.F you'll see there are entires there for ...

wilkin

If you have a linear salt profile then you should have constant d(rho)/dz and then the vertical diffusion term (d/dz*AKT*constant) would be identically zero. So, I doubt that explicit vertical diffusion is the the issue. Besides, AKT_BAK is just a background minimum to the chosen vertical turbulence ...

wilkin

Restart files include TWO time levels - the model output time step n and the proceeding time step n+1 . This is necessary to exactly restart the time stepping algorithm as if the model had never stopped.

You can add the DYE variable to a history file and start from there, in which case ROMS takes ...

wilkin

Judging from your plots the geometry of your grids is simple east-west rectangular coordinates. But the log reports a range of values for angler from -90 degrees to 270 degrees:

Grid 01:
GET_GRID_NF90 - angle between XI-axis and EAST: angler
(Grid = 01, File: /home/roms/input/NH_grid.nc)
(Min ...

wilkin

but I still have a question about what operations swrad will undergo after roms defines DIURNAL_SRFLUX. What are the characteristics of the calculated temp (what is the difference from not defining #DIURNAL_SRFLUX)

This is a good example of where going to the source code and reading the comments ...

wilkin

We are advertising for a Postdoctoral Associate in Physical Oceanography & Ocean Modeling at the home of myroms.org

The project is to work on elements of the design, development and evaluation of a prototype East Coast Community Ocean Forecast System (ECCOFS) that implements data-assimilation of ...

wilkin

No, the inputs and outputs will not be exactly equal, and depending on your options could be quite different.

In the history files, ROMS outputs the instantaneous radiation at the model time step.

If you provide input swrad that is a long term (daily, monthly ...) average, and use #define ...

wilkin

Something to check is whether your circumpolar integral of Hvom in the long term mean is zero (volume conservation). Then for a closed volume Tref doesn't matter because < Hvom > Tref = 0. It's a sanity check on your methods and a possible indicator of accumulated errors associated with the ...

wilkin

The Chapman condition does not "fix" zeta on the boundary. It is designed to radiate out barotropic waves generated inside the domain.

Read ... https://www.myroms.org/wiki/Boundary_Conditions#Chapman_boundary_condition

It is designed to be used in conjunction with the Flather condition. Read ...

wilkin

I use time varying air pressure to generate oscillating flow in a "1-dimensional" ( actually 3 x 3 horizontal grid) doubly periodic test case for teaching ... similar to the BIO_TOY example.

For this to work, I have to fill the entire grid and disable the periodic exchange.

So, in ana_pair.h note ...

wilkin

Isn't is just ...

 = - = - Tref because Tref is constant.

You have Huon_temp and Huon, so calculate Huon_temp - Tref*Huon.

wilkin

Your choice of reference time is not recommended. See https://www.myroms.org/forum/viewtopic.php?p=119#p119

Since you are working with MATLAB, you can easily convert a chosen reference time to datenum convention

In ROMS, e.g.:
TIME_REF = 20240101.00

In MATLAB:
>> dnum = datenum(2024,1,1) + ocean ...

wilkin

In bowl_safe_new2_rivers_u.nc the first river_Xposition index (501) is invalid FORTRAN indexing of rho goes from 0 to L (green points in this figure: https://www.myroms.org/wiki/Numerical_Solution_Technique) which for you is 0 to 500. The indexing of u-points (in red) is therefore 1 to 500, so 501 ...

wilkin

You say ...

defined rivers both in lateral boundary and mask areas (mask = 0), but I still have the same problem with LuvSrc.
I tried LwSrc (mask = 1),

... but do you mean u,v or rho masks?

LwSrc points can be anywhere in the water mask_rho=1 but with i,j set by river_Xposition,river_Eposition ...

wilkin

Is your initial sea level (zeta) less than -0.1 at the start. If so, Fz0 = -0.1 could be out of physical range.

If sea level is OK, I suggest experimenting with release locations in grid units instead of lon/lat.

wilkin

I believe Hsbl and Hbbl are available for output only when using LMD_MIXING. You have MY25_MIXING.

Running with both BODYFORCE and a vertical turbulent mixing closure scheme might be inconsistent. Using MY25_MIXING (or GLS_MIXING or LMD_MIXING) assumes the stress is applied as the surface boundary ...

wilkin

The BODYFORCE option is designed to directly impart momentum into the cell centers, by-passing the role of vertical turbulence to transfer surface stress into the interior by mixing.

The momentum flux is thus distributed across the chosen set of layers as a stress divergence (Tau_k+1 minus Tau_k ...

wilkin

In bulk_flux.F, once the stress is calculated according to the Fairall et al. COARE algorithm the scalar stress needs to be broken out into its u- and v-direction components.

This occurs in code ...

...
Cd=Wstar(i)*Wstar(i)/(Wspeed*Wspeed+eps)
...
! Compute wind stress components (N/m2), Tau ...

wilkin

(1) The Dsrc = 2 option is not really a vertical flux. It is a flux divergence imposed at the rho-points cell center. So, as such it is neither upward nor downward, it just appears. This is what you want if your vent is adding buoyancy but not momentum.

(2) The fractional distribution of the mass ...

wilkin

Clamping and nudging won't help you if you have an advection dominated outflow because the exterior conditions know nothing of the interior solution that is propagating out along characteristics of the current. They will act as a wall.

You might be better off simply using gradient conditions ...

wilkin

A radiation condition for 3-D variables is potentially sensible if you have 3-D wave motions emanating from the model interior, such as internal gravity waves, that might reflect at the boundary. This was the general idea behind the original Orlanski (1976) paper that inspired this approach.
https ...

wilkin

No, that is incorrect. For flow through the v face river_direction = 1 but if flow is in the negative y direction then river_transport is negative.

wilkin

If you are using analytical wind stresses, like in the UPWELLING test case, those are set in Functionals/ana_smflux.F (for surface momentum flux).

Open Functionals/ana_smflux.F and search for UPWELLING. You will see that sustr = 0 and svstr ramps up over 2 days as a SINE function. You can modify ...

wilkin

The river_direction variable specifies whether flow is crossing a u-face (river_direction = 0), a v-face (river_direction = 1) or enters as a flux divergence in the center of the cell at a rho-point (river_direction = 2). That is spelled out here: https://www.myroms.org/wiki/River_Runoff#river ...

wilkin

You have CPP flags in there I don't recognize from myroms.org .... OBC_NUDGING, CORE_FORCING, ALBEDO_CLOUD ... etc.

If this is COAWST or CROCO I suggest you post your query on those forums.

wilkin

You've got a bit of a mess of CPP directives there.

You have BODYFORCE to add the wind stress directly into the surface layers, but you also have GLS_MIXING and BULK_FLUXES to add the stress at a surface bojundary condition to vertical mixing. Without digging into the code I don't know which you ...

wilkin

you need to set NDEFHIS and NDEFAVG not equal to zero.
roms*.in file has an extensive description about these switches

That is incorrect. NDEFHIS = 0 is valid. It means create only one history file and write all history output to it.

The GLOSSARY comments at the end of every example roms*.in ...

wilkin

1. When I prepare the wind forcing, should I rotate the wind direction from geometrical eastward/northward to grid xi/eta in my forcing .nc file ?

If you are providing Uwind,Vwind in a netcdf file and having ROMS interpolate from the Uwind,Vwind lon,lat coordinates then ROMS will apply the angle ...

wilkin

It did read them. It's just telling you that somewhere in the file you have NaN values. They aren't on the boundary though otherwise ROMS would not have run.

wilkin

There is a diagram of the ROMS computational grid here:
https://www.myroms.org/wiki/Numerical_Solution_Technique

The rho cells number from 0 to L (in green), so there are L+1 values. This is the netcdf file dimension xi_rho.

In roms.in you indicate parameter Lm, which stands for L minus 1. So, xi ...

wilkin

Adding solar radiation will heat the ocean. An ice covered lake is insulated from solar heating.

wilkin

You say you want to keep the ocean surface temperature at 0 degrees ... so it freezes?

Please post a more detailed question and do so in a new thread because this no longer has anything to do with River Inputs.

wilkin

c_contact.m is for creating the contact file for two nested grids once the two grid files are already made, so I don't follow why you would be running that routine when you say you are trying to create a grid anew.

wilkin

You want the LwSrc option. This introduces tracer mass without momentum by augmenting the tracer flux divergence at the cell center. It does not impact u,v.

For both LuvSrc and LwSrc sources the rate at which tracer is added (in kg/s) is the volume flux (river_transport, m^3/s) multiplied by the ...

wilkin

Your log shows ,,,

Maximum barotropic Courant Number = 9.00832155E+00

which seems much too large, but I can't where that comes from on the basis of initial conditions.

wilkin

You might need ...

** ANA_SRFLUX if analytical surface shortwave radiation flux **

wilkin

% bbox - A 4-element vector in the format of Matlab axis defining the
% lon/lat bounding box region to subset with the OPeNDAP query
% e.g. bbox = [-110 -30 0 55] or [250 330 0 55] for West Atlantic.
% NOTE: The longitude coordinate in ERA5 breaks at the prime meridian.
% This function detects ...

wilkin

When I read these papers some time ago I concluded the quantity most representative of ROMS surface later (k = N) temperature would be to compute sea_surface_temperature minus sses_bias . That's what we do in operational Doppio.

wilkin

The AMSR file format includes a sses_bias correction term that a user can apply to obtain the best estimate of foundation temperature. See Petrenko et al. https://doi.org/10.1175/JTECH-D-15-0166.1

wilkin

Correct. I am not updating the svn version. In the contents.m ...

% This collection of tools was originally created March 2008 (some are much
% older) for release to the ROMS User community. Prior to December 2020
% the code set was managed by svn and accessible with login credentials for
% the ...

wilkin

In the real-time Doppio analysis system we run for MARACOOS/IOOS, we currently access AMSR2 L2 from dataset ID AMSR2-REMSS-L2P-v8.2 at PO.DAAC.

I'd have to check, but I believe the granules that intersect with the model domain bounding box are identified using the PO.DAAC subsetter.

wilkin

Yes, there were some changes at NCAR. I updated my code some weeks ago and it seems to be working robustly. I just uploaded the changes to my GitHub https://github.com/johnwilkin/roms_wilkin so please refresh and see if that solves the problem.

wilkin

This is explained in WikiROMS:

For uv source you must be mindful of the placement of land and water to the "left" or "right" of the cell interface:
https://www.myroms.org/wiki/River_Runoff#river_Xposition_and_river_Eposition

You need to change the sign of the river entering the ocean in the ...

wilkin

You can #undef SOLVE3D and see what happens
https://www.myroms.org/wiki/SOLVE3D

wilkin

(There was another post about this a couple of years ago ... https://www.myroms.org/forum/viewtopic.php?p=23078&hilit=thin#p23078)

Exercise some caution in simply setting mask_u(i)=0 to impose a barrier because the calculation at u(i+1) may use values at u(i-1) in addition to u(i) to complete the ...

wilkin

Do you have ...

Code: Select all

  #define UV_LOGDRAG

wilkin

This isn't a nesting issue, it's simply that you are asking for tile partitions that are too small.

See this link in WikiROMS https://www.myroms.org/wiki/File:communications.png to visualize the tiles and their associated halo regions that communicate the data from adjacent tiles to complete the ...

wilkin

The coastal ocean modeling and dynamics community has lost an influential and insightful colleague and mentor. John S. Allen authored over 100 papers https://oregonstate.box.com/s/ykejvhu8t3y935o2kvq4v6dhri98xtda on coastal ocean and shelf sea currents, instabilities, coastal trapped waves ...

wilkin

Your open boundary conditions don't make sense.

You have periodic conditions for u,v but a mix of closed and Flather conditions for ubar, vbar. If u,v is periodic then ubar,vbar must be also. Without knowing more about your grid it's hard to say what to do, but since you have MASKING defined I ...

wilkin

That looks like an issue with your libraries or compiler. If it were just ROMS the blow-up message would give more guidance as to where in the code the issue arose. This is an outright crash. Make sure you are linking to MPI libraries built with the same compiler you use for ROMS.

That said, the ...

wilkin

The trac ticket you note #943 explains ....

The ROMS source code was consolidated due to changes for new developments:

... including ...

Porting the Waves Effect on Currents (WEC) based on the Vortex Force formulation of Uchiyama et al. (2010) and Kumar et al. (2012) based on COAWST ...

wilkin

There is no specific example in the myroms repository.

When you think of how many features ROMS has it's overwhelming to construct examples for everything. FLOAT_OYSTER is pretty straightforward - I'm sure you can figure it out.

But please feel free to contribute a test case - perhaps an option ...

wilkin

Do you want to do this with Lagrangian particles, or Eulerian tracers?

With tracers:

You can do this with the sediment model. Configure it with just a single sediment class and set the sinking velocity you want. Turn off all the options for resuspension etc. (unless you ultimately those). You can ...

wilkin

Please try running this ... modified ... code again. The NCAR helpdesk claim the issue is resolved.

wilkin

This error message ...

ERROR: Cannot open file '/home/roms_project/sobhan/varinfo.yaml': No such file or directory

indicates that file varinfo.yaml is not present in that location, or its permissions make it unreadable.

Most users keep a copy of varinfo.yaml in the directory where they run ...

wilkin

If you got to that place in the code ...

Error in roms_get_era5_NCARds633_bulkflux_update (line 478)
itime = double(ncread(url,'forecast_initial_time'))/24 + epoch;

... then it must have read the first set of data files, right? So the urlbase is not wrong. The File Not Found error in this case ...

wilkin

The code at line 149 in get_data.f 90 is
& BOUNDARY(ng) % zetaG_west)
Is a problem with the boundary conditions?

It does seem that's a likely source of the problem.

I can't help here. Your boundary file is in a very different format to what I'm used to, with 1-dimensional vectors for the ...

wilkin

#6 0x564becef30db in get_data_
at /home/felipejbarragan11/roms_project/test/MyTest/ROMSSantander3D/Build_romsG/get_data.f90:149

What is the code at line 149 in get_data.f90?

wilkin

Read the Backtrace for this error.

This shows the ultimate failure is associated with trying to read from a netcdf file. But which file? The particular netcdf read was initiated at

Backtrace for this error:
...
#3 0x561b93338719 in __mod_netcdf_MOD_netcdf_get_fvar_1d
at /home/felipejbarragan11 ...

wilkin

It looks like you want to prescribe temperature and salinity surface boundary conditions with analytical functionals (because you have #define ANA_STFLUX). In that case you also need #define ANA_SSFLUX, otherwise the default action is to look for the data amongst the forcing netcdf file list.

wilkin

It appears NCAR changed the hostname for the THREDDS server. It is now https://thredds.rda.ucar.edu not just https://rda.ucar.edu.

If you update the line in the code that sets urlbase to read ...
urlbase = 'thredds.rda.ucar.edu/thredds/dodsC/files/g/ds633.0';
... it should work.

It also ...

wilkin

Thank you for alerting us to this. Something has broken very recently following recent updates. We are looking into it.

wilkin

I showed a way to do this with MATLAB in this post:
viewtopic.php?p=23769#p23769

wilkin

Well, yes, air-sea flux options and inputs is just another one of the many things we need to add to WikiROMS. Perhaps in my retirement ... !

Scanning those entries in varinfo.yaml I do spot an error in the #Input tag for dQdSST:

- variable: dQdSST # Input
standard_name: derivative_of_surface ...

wilkin

The comments in varinfo.yaml tell you the units ROMS expects in a netCDF file, # Input: [Watt/m2]

- variable: shflux # Input/Output
standard_name: surface_downward_heat_flux_in_sea_water
long_name: surface net heat flux
units: watt meter-2 # Input: [Watt/m2]
field: surface heat flux ...

wilkin

I don't follow the question.

If ANA_PSOURCE is not defined but logical flags in roms.in indicate sources are to be used, then the netcdf file is interrogated.

If ANA_PSOURCE is defined, an arbitrary allocation of 200 is accommodated (there should probably be a warning in ana_psourceF about this).

wilkin

Looks like you missed this part of the guidance at https://www.myroms.org/wiki/River_Runoff ...

If passive tracers are being used, by having set #define T_PASSIVE and the associated options in roms.in, then for those tracers ...

In your tracer_1r.h file you need ...

#define T_PASSIVE

Also, I ...

wilkin

Look at the figure here:
https://www.myroms.org/wiki/River_Runoff#river_Xposition_and_river_Eposition

It is for you, as user, to decide whether your source points falls on a coastline with land on the left or right. If it is Case A (in the figure), then river_direction = 0, river_Xposition = i and ...

wilkin

You may be off by one in your river_Xposition etc if you have not adjusted for FORTRAN indexing as opposed to Matlab. This is explained in WikiROMS https://www.myroms.org/wiki/River_Runoff#river_Xposition_and_river_Eposition.

roms_plot_river_source_locations.m from my ROMS tools at https://github ...

wilkin

(3)
and OBCFAC == 0.0d0 Tnudge=0. I think this means that I didn't use the nudging.

If you set the nudging time scale to zero, then yes, the "Nud" part of "RadNud" does nothing. The temp and salt data in your BRYNAME list is ignored. Moreover, even if you had a finite Tnudg your OBCFAC=0 setting ...

wilkin

I'm struggling to envision what geophysical problem it is you are trying to model, but if you want to add mass (volume) at some rate (cubic meters per second) distributed all across the surface of your lake, then you have to create all the point sources. It sounds like this is an idealized problem ...

wilkin

This is not a use case that anyone has run, to my knowledge anyway.

In principle it should work, because the volume flux is introduced via Nonlinear/step2d_LF_AM3.h (search for LwSrc to see where it is done). There are no loops over k, of course, and so the additional information like river_Vshape ...

wilkin

I have written code (roms_isosurface.m) to use Matlab's isosurface function to create 3D visualizations. This is part of my roms_wilkin toolbox at https://github.com/johnwilkin/roms_wilkin

Some examples are shown below generated by roms_isosurface_demo.m in the GitHub repo.The demo calls an OPeNDAP ...

wilkin

There is some additional, admittedly very succinct, documentation at https://www.myroms.org/wiki/floats.in

As an example, the code below will initialize 15 floats at i = 3 (Fx0 =3), j = 2 (Fy0 = 2) and k = 1:2:31 (Fz0 = 1, Fdz = 2):

! C Initial horizontal coordinate type (0: grid units, 1 ...

wilkin

A paper describing larval transport results using ROMSpath in a nested grid configuration is here:
Garwood, J.C., Fuchs, H.L., Gerbi, G.P., Hunter, E.J., Chant, R.J. and Wilkin, J.L., 2022. Estuarine retention of larvae: Contrasting effects of behavioral responses to turbulence and waves. Limnology ...

wilkin

Another approach you might consider is offline particle tracking using ROMSPath https://gmd.copernicus.org/articles/15/4297/2022/

ROMSpath uses almost exactly the same floats algorithm internal to ROMS. Importantly, the velocity is expressed in ROMS native xi,eta,s coordinates using u,v,omega, so ...

wilkin

In my experience LtracerSrc = F is quite stable when used with Dsrc = 2, i.e. the updated LwSrc cell-centered divergence method for sources.

wilkin

romsS is a S for "serial" meaning you intend to run on a single processor.

romsM is a M for "multiple" or "MPI" meaning you will run on multiple cores by setting nNtileI and NtileJ in roms.in greater than 1. Their product is the number of concurrent (parallel) cores to run on with MPI.

Whether or ...

wilkin

If you want to set the concentration of biology tracers in the rivers, i.e. use LtracerSrc = T for some variables, then those concentrations are determined by:

Default is to read the file specified by SSFNAME in roms .in
! Input Sources/Sinks forcing (like river runoff) file name.

SSFNAME ...

wilkin

ROMS should create the file defined by HARNAME. When we first introduced this option, ROMS would add to the forcing harmonics file but that was not an option once nesting was introduced. See ticket https://www.myroms.org/projects/src/ticket/736

To get to that line in the def_tides.F code where you ...

wilkin

Posted on behalf of Glenn Flierl, MIT

Dale and I talked during a time when we were both visiting NCAR; it was great seeing someone bringing more advanced numerical methods to ocean models, albeit idealized, that began to represent the eddies. This expertise in both the physics and the numerics ...

wilkin

The Boulder "Daily Camera" newspaper published this obituary online today ...
https://www.legacy.com/us/obituaries/da ... d=51548008

wilkin

You have ...

SOLVE3D Solving 3D Primitive Equations

so you need enough vertical layers for ROMS to be able to compute vertical gradients for such things as vertical mixing and hence to impose surface and bottom stresses, etc.,

If you want a strictly 2-D solution where you are solving the ...

wilkin

This code already exists ... #define SOLAR_SOURCE. The work is done in Nonlinear/lmd_swfrac.F using two length scales according to the method of Paulson and Simpson (1977) https://doi.org/10.1175/1520-0485(1977)007%3C0952:IMITUO%3E2.0.CO;2

The user chooses the length scales with the water type ...

wilkin

My suggestion was to #define LONGWAVE_OUT and #undef LONGWAVE. Don't set both - that makes no sense. Look in bulk_fluxes.F to see what these options do.

LONGWAVE_OUT will cause ROMS to read lwrad_down (or whatever you renamed it to be in varinfo.dat) but compute outward longwave from model SST ...

wilkin

I see you are using net longwave radiation forcing.

Does IMDAA have a separate analysis for downward longwave only (i.e. what ROMS expects as lwrad_down)?

If it does, you might find #define LONGWAVE_OUT helps because your growing SST might be reined in by increasing upward outgoing longwave ...

wilkin

On behalf of Ichiro Fukumori ...

I first met Dale in the early ‘90s. Friendly and kind, and generous with his time, Dale introduced me to ocean modeling for exploring ways to conduct ocean data assimilation. The model was SPEM, which I found to be elegantly written in a manner one could readily ...

wilkin

The oceanography community has lost a respected pioneer of regional and coastal ocean modeling. Dale Haidvogel passed away unexpectedly on March 14 at his home in Buffalo, NY. Though recently retired from Rutgers University, Dale was still active as an emeritus instructor in graduate GFD courses.

A ...

wilkin

Is that the salt initial condition you want? I doubt it. I think you need to sort out the basic solution before you worry about dyes.

And why all the little lakes in the land?

wilkin

You could simplify your test case by choosing Rad instead of RadNud for the dye open boundary condition ... then ROMS should not need use (or even read) the boundary dye values.

You have OBCFAC = 0.0, which I think disables the nudging on inflow. Maybe start with OBCFAC = 1.0 so you have boundary ...

wilkin

It is hard to tell from those diffs, but I suspect you might have confused ROMS in the restart process.

It looks like you are asking ROMS to write 10 records to a history before refreshing a new file, but you are restarting only part way through that sequence. ROMS might be looking for file 0002 to ...

wilkin

Parker and Aurora,

I don't think this is your issue, but could you just check whether the point source giving you trouble happens to fall on a parallel tile boundary? If so, or even if not, can you test whether the blow-up you get also occurs at the same time in the run for a different choice of ...

wilkin

If there was a bug here, and I haven't yet tunneled into the code to double check, it would affect whether the source was identified in the active i,j tile, or not, for all time. I don't see how it could be the direct cause of an instability that appears well into a run.

wilkin

I think you are confusing dimensions and coordinate variables.

It does not matter whether a variable is dimensioned salt_time or mud_time. That conveys no information to ROMS about the time coordinate . Which time coordinate to use to interpret the index into the time dimension can be specified ...

wilkin

Yes, as I noted previously
Since the dz is already in there, the vertical integral of the tracer flux, weighted by cell width, is simply the summation from k=1,N. To also integrate over "x", just take the sum over "i" because the dx is already in there.
so you should sum over layers to get a ...

wilkin

What you say is, in principle, correct. Generally speaking, in method B you need to subtract the mass of tracer at the start, but that is zero for this specific set up, right?

How are you calculating dz in the volume integral? It is easy to get that wrong if you are coding from scratch.

Replying ...

wilkin

We are advertising for a Postdoctoral Associate in Physical Oceanography & Ocean Modeling

The project is to work on elements of the design, development and evaluation of a prototype East Coast Community Ocean Forecast System (ECCOFS) that implements data-assimilation of real-time observations in a ...

wilkin

You want ...

v.dx.dz
v.T.dx.dz

The diagnostics outputs Hvom tracer are exactly these terms. In "Hvom" the "H" does not refer to water depth, h, but to the "Hz" variable that is ROMS time-varying layer thickness, i.e., dz. (The sum k=1 to N of Hz equals h+zeta.) The "o" refers to "over", and "m ...

wilkin

My own guess is that there is some bug in the zeta calculation in the presence of a vertical source that leads to an instability. I can't see how a baroclinic error would show up as a huge SSH jump.

Hi Parker,

Once you get an instability in the velocity the extreme values drive extreme ...

wilkin

I think the sudden jump in sea level and surface velocities are symptoms of the instability, and not its origin.

Looking at the density profile and AKt vertical mixing eddy diffusivity immediately before and after the vertical homogenization of the water column can you convince yourself this is ...

wilkin

I see your point.

The individual files, e.g., https://rda.ucar.edu/thredds/dodsC/files/g/ds083.3/2021/202101/gdas1.fnl0p25.2021010100.f09.grib2.html (one every 3 hours)
have variable Downward_Short-Wave_Radiation_Flux_surface_3_Hour_Average

Perhaps that is what should have been aggregated in ...

wilkin

Brainstorming here ...

If your Vshape causes the source to flow entirely into the bottom cell, then a transport Q is going to drive vertical flow up out of that cell at velocity w = Q/(dx*dy). If the thickness of that cell, dz, is small, such that w*dt/dz > 1, then you might have a CFL instability ...

wilkin

I think we might be speaking at cross purposes here.

You are plotting values extracted from individual files, which have 3- and 6-hour average values in each cycle. I was plotting "3-hour" data from the THREDDS FMRC aggregations. I may be wrong, but I assumed whoever configured the FMRC correctly ...

wilkin

Looks to me like your magenta dots are simply half way between the orange dots, and if you plot them together you get a nice diurnal cycle resolved at 3-hour intervals? I think the files you are reading have correctly parsed data from the 3 and 6 hour fields in each cycle. No further tinkering ...

wilkin

As I understand it, the reftime tells when that forecast cycle started. There is a new forecast every 6 hours, and data from that cycle is saved at 3-hour intervals. Hence, reftime repeats in pairs 0, 0, 6, 6, whereas the time coordinate advances by 3 hours for each record.

This is a Forecast ...

wilkin

I have added codes roms_{get,write}_gfs_NCARds083_bulkflux.m for accessing the NCAR archive of GFS data to my tools https://github.com/johnwilkin/roms_wilkin. They run similarly to the ERA5 downloads, though there is presently an issue with updating the GFS THREDDS aggregation that the NCAR Helpdesk ...

wilkin

The MATLAB codes I wrote (roms_get_era5_NCARds633_bulkflux.m and roms_write_era5_NCARds633_frcfile.m) at https://github.com/johnwilkin/roms_wilkin extract hourly data for all fields from the ERA5 dataset in the Research Data Archive at NCAR (rda.ucar.edu).

wilkin

The error says "Permission denied" which suggests you are trying to write to a non-existent directory or you lack the permissions to create files there.

Matlab has great debugging tools. Try:

>> dbstop in write_roms_otps_ncfile_v2

and when Matlab pauses there, check what directory you are in ...

wilkin

The 4th order predictor corrector scheme for ROMS online floats interpolates the gridded velocity field to each float's position within a cell and integrates over the model time step DT in such a way that as a float approaches a land boundary the velocity normal to the boundary tends zero (by virtue ...

wilkin

Is the bathymetry continuous at the nest perimeter? If not, it's not possible to match velocity ubar,vbar with transport ubar*h,vbar*h and volume will not be conserved. Two-way nesting compensates for this by establishing a pressure gradient that will drive flow into the nest to rebalance the sea ...

wilkin

Graduate Assistantship in a NOPP Internal Wave Project (Physical Oceanography)

The Division of Marine Science of the University of Southern Mississippi (USM) invites applicants for an M.S. or Ph.D. study in Physical Oceanography starting in summer/fall 2023. We seek a skilled and motivated ...

wilkin

I think you need to take some time to carefully quality control your initial conditions.

Quickly looking at one of your log files 2018/3600/..out2 I see initial condition values of:

- potential temperature
(Min = -7.18577153E+00 Max = 3.30190636E+01)
- salinity
(Min = -4.23978987E+00 Max = 4 ...

wilkin

(Min = -1.09291988E-05 Max = -4.45711885E-06) regrid = F

Those seem like very small numbers for dQdT.

varinfo.yaml indicates that the assumed units are watt meter-2 Celsius-1

ROMS works in kinematic units so this factor is divided by rho0*Cp to convert it to m/s. Then in set_vbc.F the ...

wilkin

Are you using #define ANA_DQDSST to set dqdt with the analytical functional file? If so, check analytical.f90 in the Build directory to see what code is actually there. Maybe dqdt is set to zero. If you are reading it from netcdf, the logfile should report the vales read. Check for that.

To revert ...

wilkin

If you are using the myroms.org code then the file varinfo.yaml sets the default names for time coordinate, sms_time, which is not in the file you made with croco tools.

The myroms code also allows over-riding the default name with a "time" attribute for each variable, but you don't have that ...

wilkin

The Mathworks have several online tutorials for getting started with using Matlab
https://www.mathworks.com/learn/tutoria ... nramp.html

wilkin

Stefan,

I see what you are saying. The reduced physics option to estimate ubar(Istr) (the physical boundary u-point) for use in Flather (when ubar_west is not provided) is solving d(ubar)/dt = -g*d(zeta)/dx and there appears to be a factor of 1/2 multiplying the pressure gradient RHS term.

This ...

wilkin

Those CPP options don't look like myroms.org or COAWST. Are you running ROMS Agrif or CROCO? You might get more help from that user community.

wilkin

I believe the 1/2 is there because, while the staggered grid formally places zeta_west 1/2 cell to the "left" of the u-face point (this is the code block where UV_TIDES is undefined so we are substituting a dynamical approximation for it), when zeta_west and ubar_west are created by users they are ...

wilkin

That does look like a bug.

There aren't a lot of of users running 2D ROMS, so it's possible some bug crept in with file handling.

wilkin

From those plots it does appear that the English Channel is where blow-up is likely to occur, but the error flags grid point (i,j,k) = (0,0,0), not a local blow-up point.

Your initial and boundary values appear to be benign, so it is odd that the blow-up occurs on the very first step and is ...

wilkin

All that said, I'm re-running the model now as it was to replicate the error and it's not crashing, which is in some ways scarier than the error I was getting previously.

The fact that your earlier run terminated with no error report from ROMS does suggest it was a crash of one of the processors ...

wilkin

One perplexing thing is the way in which your job exited. There are quite extensive error trapping routines in ROMS so that you get some information on where the code failed before it exits. Yours just stopped, judging by the log you posted. It did not indicate a BLOW UP, or give a line number for ...

wilkin

Does the log really end abruptly like that, with no error reporting from ROMS?

wilkin

You are correct. Something has changed with Matlab netcdf handling. I'm using R2022b (the latest) and get the same error.

Your fix to add netcdf.reDef at line 100 works

% svn diff nc_drename.m
Index: nc_drename.m
===================================================================
--- nc_drename ...

wilkin

That looks like a write permissions error.

ROMS is not reading the restart file ... because you have ANA_INITIAL and it creates the initial conditions internally.

The log shows the error occurs in DEF_RST_NF90, i.e., when ROMS is trying to define/create the restart file.

You may not have ...

wilkin

Yes, that's a bug I just discovered myself. There are some other similar bugs but not for "spherical" grids using lon/lat coordinates.

wilkin

The ROMS sign convention is that positive heat flux warms the ocean. This is so for swrad, lwrad and lwrad_down.

In your region you would typically have 400 W/m^2 of downwelling longwave radiation https://link.springer.com/article/10.1007/s42452-019-1172-2/figures/4 so the fact you have negative ...

wilkin

If you really want to do this, the logic would be to perform the advection time step update for all tracers except 6 with ...

IF (itrc.ne.6) THEN
t(i,j,k,3,itrc)=Hz(i,j,k)*(cff1*t(i,j,k,nstp,itrc)+ &
& cff2*t(i,j,k,nnew,itrc))- &
& cff*pm(i,j)*pn(i,j)* & ....

When the code gets to itrc = 6 it ...

wilkin

Found Error: 02 Line: 1006 Source: ROMS/Nonlinear/get_data.F

Look at line 1006 in your get_data and see what variable it is trying to read.

In mine (not COAWST) it's reading boundary conditions so yours must be quite different. Presumably, it's related to your modified setup.

wilkin

Often, compiling with debug option leads to more informative error messages. Try that first.

wilkin

You definitely want to remove the perimeter wall.

Then you need some sensible open boundary conditions. Presently you have east-east periodic, which would only be admissible in this situation if you contrived to match the land/sea mask at the left and right boundaries, or ensured both were ...

wilkin

You haven't #define MASKING

wilkin

Thank you for your reply. Is it possible if I replace nan with 0 in the file? And will roms do the screening operation based on mask?

Yes, zero values (any finite value) in the land should solve your initialization problem. Then the masks will be applied - zero times anything finite is zero ...

wilkin

You can't have NaNa in any input fields, not even in the area masked as land. ROMS computes on all points and then applies the 0/1 land mask afterward. Zero times NaN is NaN, so they won't go away.

Your log reports that both forcing and boundary data were read for initial time 2019-08-01 00:00:00 ...

wilkin

Your air pressure data are in Pascals ranging from 98,040 to 102,590. By default, ROMS expects Pair in millibars unless you have modified varinfo.dat to do the conversion via a scale factor of 0.01. It is likely you have crazy air pressures in your run.

Quickly plotting your Pair data I find that ...

wilkin

Read the documentation.

The | joins files with the same variables that are broken into time segments.

You want \ which goes to the next file - here you have different variables in each file in a list. You have 7 files.

wilkin

The "|" character is used to effectively concatenate files with sequential time records, e.g. one year at a time. When this is used it counts as only one file.

You should use the "\" character to end one file and go to the next, with a different variable, in a list.

What you have given is ...

wilkin

Your file output.txt indicates that the ERA5 forcing data were read successfully.

The next thing ROMS is doing when the error occurs is reading the open boundary conditions file, and the first variable it reads is zeta, so that's why you are getting an error about zeta_time.

The log says your ...

wilkin

The netcdf files created by roms_write_era5_NCARds633_frcfile.m have a time coordinate variable called "time" (not zeta_time) and the "time" attribute of all variables will point ROMS to read "time" in preference the default in varinfo.dat.

But, I can't imagine what meteorology forcing file would ...

wilkin

I have updated the svn repository, but at some time in the medium term future we will shut that down and keep everything just in GitHub.

wilkin

You are correct that the svn version is out of date. I'll try to get that updated soon, but the latest version will always be at GitHub for the foreseeable future.

wilkin

Do you have a current version of this function from https://github.com/johnwilkin/roms_wilkin ?

I updated it 4 months ago in response to a similar post from another user, and I just tested that ...

>> E = roms_get_era5_NCARds633_bulkflux(2020,1,[99 140 0 26])

works just fine with the present ...

wilkin

Is the error repeatable?

It looks like a failure on the server side which might be a temporary problem at NCAR.

wilkin

I've been making simple rectangular equi-spaced grids in Matlab using an arbitrary Lambert projection which you can center on any point and rotate. There is a post here: https://www.myroms.org/forum/viewtopic.php?p=23769#p23769. I have not tried, but I don't see any reason it would fail spanning the ...

wilkin

I thought I had expunged all the old calls to SNC tools functions, but missed this one.

The first line just needs to be:

Code: Select all

cpp = ncreadatt(file,'/','CPP_options');

I've updated the repo on GitHub.

wilkin

The University of Arizona, Princeton University, and NOAA’s Geophysical Fluid Dynamics Laboratory (GFDL) are seeking applicants for a position to analyze simulations of the Southern Ocean from GFDL’s CMIP6-generation ocean-only and coupled climate models using observational constraints. The position ...

wilkin

Where did you get your source code?

My up-to-date checkout from myroms.org has flt_test.h at revision 1099, and the #define FLOATS is definitely there.

/*
** svn $Id: flt_test.h 1099 2022-01-06 21:01:01Z arango ...

wilkin

#define FLOATS

then ROMS reads the "floats.in" file name set by option FPOSNAME in roms.in. There is an example floats.in file in External that is documented further on WikiROMS here
https://www.myroms.org/wiki/floats.in

wilkin

You chose 4-th order centered advection for tracers. Is that what you used previously with the stable model?

Your barotropic Courant number is 0.6 which is getting a bit large, though still stable. You might need to increase NDTFAST.

wilkin

Just code it yourself using the scatteredInterpolant class:

% some ROMS output
f = 'http://tds.marine.rutgers.edu/thredds/dodsC/roms/doppio/2017_da/his/History_Best';

% get some zeta output
t0 = 1000;
tcount = 100;
zeta = ncread(f,'zeta',[1 1 t0],[Inf Inf tcount]);

% get the lon/lat coordinates ...

wilkin

I notice you have salinities in the initial conditions (and boundary conditions) of minus 11

NLM: GET_STATE - Reading state initial conditions, 2005-08-27 00:00:00.00
(Grid 01, t = 53609.0000, File: Katrina_ini.nc, Rec=0001, Index=1)
...
- salinity
(Min = -1.15371107E+01 Max = 3.68318607E+01 ...

wilkin

I get that error recomputing the contacts points file for my DOPPIO and PIONEER application. I had previously computed the file some time ago with an earlier version of contact.m.

So, yes, it looks like a bug has crept in.

wilkin

I don't get that error. At that place in the code u_west and v_west are the same dimension because they have been interpolated to the receiver rho points (for rotation). On exit from obc_roms2roms they have been placed on their respective grids so have different dimensions.

But, my code has a few ...

Search found 937 matches