Begginer Questions about Modifying ROMS: Forcing & Domain

[robjenkinsiii]

Hi All,

I'm still relatively new to using ROMS and I have some thoughts on taking little strides to get familiar with the model by modifying existing cases.
I've been reading a lot of the literature and opening up the include and ana_XXX files etc, to get an idea of how things work but I'm still really stuck and don't know where to begin.

Specifically what I want to do is to take a test case, Specifically the UPWELLING test case and change:
- Domain
- Initial Conditions
- Boundary Conditions
- Forcing
- Grid Size

I think I've identified in 'ocean_upwelling.in' where to set a forcing file, and I see in this file how to change the grid but not the domain.
Furthermore, while I think I see where to set a forcing file, I've no idea how to set one.
What I want to do is set the wind stress analytically (The equation for stress is based on something like a hurricane, where stress depends on the radius from the center of the domain.

Thanks in advance for any helpful words or tips to get going.
It may seem obscure, but I'm doing this specifically to compare the ROMS results to another model.

[kate]

I suggest you start with the FAQ. The switches between analytic formulae and files are C preprocessor options starting with ANA_. The files to change are all of the form ana_xxx.h, ana_grid.h for example.

[robjenkinsiii]

Thanks for responding, Kate.

Your direction to the FAQ's is helpful and it reaffirms everything I've read elsewhere.
I would like to ask another question just to make sure I'm understanding this process correctly.

Setting the cpp flags/options only tells the compiler where to look for some certaininformation, correct?
You say the files to actually change are of the ana_~~~.h family.

I open up ana_smflux.h where the FAQ state momentum fluxes can be set analytically, but do I not also need to create a forcing file?
I see where in that script, the u & v are supposed to be set analytically but I'm not 100% about how so. Meaning

You point to ana_grid. I see in that file that it controls the domain, however I'm confused as to how to change any of the domain sets.
They're written like so:

Code: Select all

#if defined MY_APPLICATION
      Xsize=???
      Esize=???
      depth=???
      f0=???
      beta=???
#else
      ana_grid: no values provided for Xsize, Esize, depth, f0, and beta.
#endif

With ??? as part of the actual code and I don't understand how or where those variables are set for, for instance, the upwelling/downwelling test case.

The analytical have a similar "???" syntax which I can't make sense of.

[kate]

If you tried to compile that code with the ???, you would get compiler errors. This is a signal to you to do something different there. If you want to start with the UPWELLING case, search all the files for the string UPWELLING. The UPWELLING case is described in the ROMS manual. It also covers the C preprocessor, as does the ROMS wiki.

[robjenkinsiii]

Okay that's wonderful!
I feared the "???" was something far more cryptic and horrible.

What I'll try to do now is to define "ANA_GRID" in upwelling.h, then set some different domain for the upwelling case in ana_grid.h. I'll post the results here.

Thanks again!

[robjenkinsiii]

So I simply changed

Code: Select all

#define ANA_GRID

to

Code: Select all

define ANA_GRID

and set some values in ana_grid.h not too far off from the original ones for the test case.

I typed ./build.bash into the command line and it complains:

Code: Select all

$ ./build.bash
rm -f -r core *.ipo /home/Robert/software/roms/Projects/Upwelling2V/Build /home/Robert/make_macros.mk
/home/Robert/make_macros.mk:1: 
*** missing `endef', unterminated `define'.  Stop.

is that because of the change I made to ana_grid?

EDIT:: I undefined ANA_GRID and the compilation ran until a later complaint about my ana_grid file. So I suppose that answers my original question but raises new ones.

[rduran]

You should have noticed that all cpp options that are defined have

Code: Select all

#define XXX

while all options that are undefined have

Code: Select all

#undef YYY

.

# is not used to comment a line.

[robjenkinsiii]

Thanks, I did notice this fact and moved ahead accordingly.

[robjenkinsiii]

I've been playing around now with definitions and analytical files for a while and I still feel like I'm missing something fundamental.

If I simply build and run the test case, I can see what analytical header files are used. Next, I tried two seemingly equivalent arrangements but saw two different results.

1) Copy & modify only ana_grid.h in my working directory.
I want to begin my changing the domain, so in ana_grid.h I define and undef "MY_APPLICATION" and set the parameters I want. changing from the original 80km to 100km. I also change the depth from 150m to 200m and give it some coriolis parameter so that it does not complain about the "???"

The result runs but I lose the bathymetry & see no variation in temperature with time or space. It is simply 14 degrees everywhere in the domain. (I see that 14* is set in ana_initial. h which I did NOT modify.)

2) Copy all used analytical header files to working directory. Modify only ana_grid.h

This time my initial time point shows the temperature distribution as before. According to the output however, it blows up and stops running. Again I lose the bathymetry.

What I don't understand is how, if at all, my two situations are different. In each all I did was to change ana_grid.h file like so.

Code: Select all

#define MY_APPLICATION
#if defined MY_APPLICATION
      Xsize=100000
      Esize=100000
      depth=200
      f0=0.0001
      beta=0.0001
#else
      ana_grid: no values provided for Xsize, Esize, depth, f0, and beta.
#endif
#undef MY_APPLICATION

One major misunderstanding is how the original test case defines the domain if MY_APPLICATION is never defined originally. Is there a seperate input file not listed among the analytical?

[kate]

The string MY_APPLICATION is simply a place holder for the name of your current domain (like the ??? thing). If you have run UPWELLING, that is the name of a specific setup. The string UPWELLING is not set in ana_grid.h, but higher up, more globally. I set my string to be say ARCTIC in a makefile.Arctic while others set it in the build script. In any case, the system then looks for an arctic.h file to include with all its C preprocessor definitions.

You can see what's defined by watching the compile happen:

Code: Select all

/usr/bin/cpp -P -traditional -DMPI -DLINUX -DX86_64 -DPGI -D'ROOT_DIR="/import/home/u1/uaf/kate/feedme"' -DNEP6 -D'HEADER="nep6.h"' -IApps/NEP -D'ROMS_HEADER="nep6.h"' -DNestedGrids= -D'ANALYTICAL_DIR="Apps/NEP"' -D'MY_ANALYTICAL="on"' -D'SVN_REV="exported"' -IROMS/Include -IApps/NEP -IROMS/Nonlinear -IROMS/Nonlinear/Biology -IROMS/Nonlinear/Sediment -IROMS/Utility -IROMS/Drivers -IROMS/Functionals -IApps/NEP -IMaster -ICompilers -D'HEADER_DIR="Apps/NEP"'  ROMS/Nonlinear/ice_frazil.F > /center/w/kate/Build_NEPG/ice_frazil.f90
ROMS/Bin/cpp_clean /center/w/kate/Build_NEPG/ice_frazil.f90
cd /center/w/kate/Build_NEPG; /usr/local/pkg/openmpi/openmpi-1.4.3.pgi/bin/mpif90 -c  -g -C -O0 ice_frazil.f90

Here you can see that MY_APPLICATION is actually set to NEP6, with header file nep6.h. It is looking in the directory Apps/NEP for the ana_xxx.h files because that's what I told it to do.

Note that this is actually three commands, the first invoking cpp, the second cleaning up the mess (blank lines) left by cpp, the third invoking the compiler.

[robjenkinsiii]

Hi Again.

I think I see what you mean about MY_APPLICATION being defined "higher up" and so I've tried something new.

This time I copied ana_grid.h and made an ana_grid2.h where in I change the domain just as described above and place this file under the ROMS/functionals directory.

I also add #define ANA_GRID2 to upwelling.h and undefine ANA_GRID.

The build runs smoothly but upon execution I get this error only a moment or two after the run begins.

Code: Select all

 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

 READ_PHYPAR - could not find input file:  ocean_grd.nc

 Elapsed CPU time (seconds):

 ROMS/TOMS - Output NetCDF summary for Grid 01:

 ROMS/TOMS - I/O error ............... exit_flag:   4

 ERROR: I/O related problem.

Shouldn't the execution create ocean_grd.nc? Any tips on sources for this kind of error?

Further, if I'm having I/O errors after changing the domain should I assume that the problem is that it is because I'm adjusting the domain while leaving everything else unchaged (i.e. there are grid disagreements between my domain, and forcing etc.)?

[kate]

If you look at analytical.F, you will see:

# ifdef ANA_GRID
# include <ana_grid.h>
# endif

Because you have undefined ANA_GRID, it will now look for an input netcdf grid file. Instead, you want to keep ANA_GRID defined and tell it how to find the appropriate ana_grid.h for your domain. Or you could change analytical.F to have:

# ifdef ANA_GRID2
# include <ana_grid2.h>
# endif

but that's now how anyone else has done it (that I know of). Also, this would require changing the code to see if it needs to read a grid file in initial.F. Better to learn to use ANA_GRID.

[kate]

Also, instead of editing upwelling.h, I suggest you copy it to say hurricane.h. Then set your application name to HURRICANE. Then in ana_grid.h, have the #ifdef HURRICANE code. Later, you can search for your modifications by searching on the string HURRICANE.

[robjenkinsiii]

I've done some playing around and now I have a new sort of problem.
ROMS compiled fine but the execution blew up after one time step and stops.

A plot of that time step shows stratification, no bathymetry (okay for my purposes) and my new domain.

Is it likely that this blowing up is because of the parameters I set in my ana_grid?
Is it possibly due to some grid disagreement because of the changed domain?

I don't think it simply has to do with the way I changed application names.

Thanks. I aplogize for this being so long. Below is exactly what I did:

• I changed Analytic.f to call my new ana_grid2.h as you suggested. There was one complaint about undefined variables in ana_vmix.h because to no application being defined. To correct it I changed "UPWELLING" to my new application name "HURRICANE" (used that name suggestion for simplicity) the name of which was defined in my buildscript.
  
  Now I compiled, ran and I got a solution with a constant temperature & no stratification & no change in time.
  
  The domain however was successfully changed.
  
  I guessed this was because I made no changes to the other ana_files involved.
  
  Using what you just told me I tried to be clever.
  I saved originals of all of the ana_ files called by this case as given by the output.
  
  Then I added If statements for "HURRICANE" that were identical to segments for "UPWELLING" in any files which had a place to do so. (Of the 7 called, I found places to make this change in 4.)

My hope was that by doing so I could approach the original test case, with maybe only the domain changed.

Compilation was smooth, no complaints. Execution blew up after first time step.

[arango]

I think that you need to read more about the instructions that we provided in WikiROMS. In particular, you need to understand the build.sh or build.bash scripts.

There is NOT a reason to rename any of ROMS header files (*.h) and/or analytical files (ana_*.h). Just copy all what you need to a working directory including the build script. Then, set MY_ANALYTICAL_DIR in the build script to your working path Your header files and analytical files in the working directory take priority during compilation and not the ones in your source code directory containing the official version of ROMS. That's all. Several of us have mentioned this several times in this forum.

If you download the ROMS test repository, you will see what I am talking about. This information is also in WikiROMS:

svn checkout https://www.myroms.org/svn/src/test MyTest

See the upwelling sub-directory in that repository. I don't understand why some users keep ignoring this repository. All the ROMS test are there. They are separated from the source code. There is not a reason to modify the version of ROMS downloaded from the repositories, except for new algorithm developing or modifications to the numerical kernel.

All the user of ROMS, novice and experts, need to read slowly all the provided instructions and released notes (new and all) describing a lot of information and tricks. You need to do your homework and learn the basic stuff. User's need to pay attention to the forum messages and read other people posts and answers so we don't have to repeat this over and over...

[robjenkinsiii]

Thanks for the corrections and advice.I was able to see in the output the priority given to the analytical files located in my working directory once MY_ANALYTICAL_DIR was defined accordingly.Since then I've been successful in playing with the test case and running the model with sensible outputs.

The last major question I have is in regards to forcing. I realize this question has been asked many times before, but I think my question is a variation.

In ana_smflux.h I've already played some with the stresses and I understand that I could use some kind of function here rather than the constant winds prescribed.

I have a wind stress field that is radially defined ( I'd be happy to show the form if you'd like to see) that I want to apply.

By my understanding the MATLAB scripts for making forcing in repository are more for real-world forcing data and don't seem suitable for what I want to do. Is this correct thinking? I'm trying to avoid using these.

I would like to ask another question about defining this wind in ana_smflux.h or another ana_*.h file of my own making but I'm not sure how to phrase it. I have this radially defined stress. I would need, I think, to build it first and decompose it into U & V stresses.

I can't even begin if I don't have real information from the grid. I see a lot of intent(in) variables here but I'm not sure that I can discern if any among them are the actual meter values of grid points which is exactly what I need to proceed. I'm sorry I know I'm not being perfectly clear.

My question is about how to call actual values from the grid in ana_smflux.h or a similar file if they are not already being called by default.

[kate]

If you look at ROMS/Functionals/ana_smflux.h, it has:

Code: Select all

#ifdef SPHERICAL
     &                      GRID(ng) % lonr,                            &
     &                      GRID(ng) % latr,                            &
#else
     &                      GRID(ng) % xr,                              &
     &                      GRID(ng) % yr,                              &
#endif

The GRID structure also contains either lat,lon or x,y on the other parts of the staggered grid (those above are at rho points). If SPHERICAL is not defined, you might already have all the x,y information you need. If you do want x,y on the u,v points, simply pass them in the same way, then add them to the ana_smflux_tile routine likewise. You've also got the angler variable, containing the grid's deviation from due east,north at each point.

[robjenkinsiii]

Thanks for those tips. I see what you mean about passing variables to ana_smflux_tile subroutine, which is shown in the code just below the part you highlighted.

From reading through ana_smflux, and looking at the grid design according to the ROMS wiki, I'm wondering also if it matters weather the sustr and svstr are defined at RHO points or at U & V points?

Many of the stresses seem to be defined using 'xr' and 'yr' and are calculated over loops defined by JstrR to JendR and Istr to IendR. So I am right to think it's okay to calculate stresses as if they are at RHO points?

Also, would you mind clarifying 'angler?' You said it contains "the grid's deviation from due east,north at each point."

Did you mean East or North?

[kate]

Many fields are on the rho points, including winds. However, the model wants sustr and svstr to be on u and v points, respectively. If it's easier to do on rho points, compute them there, then average to where they are needed after any rotation (if needed).

angler is the angle between the i-direction and east, also between the j-direction and north.

[robjenkinsiii]

That is more information from angler than I would have expected. 2 value from one point. Your description of that variable suggests to me that angler for the upwelling channel would be zero everywhere as the i direction and j direction are aligned to east and north for this case throughout three domain. I might be still misunderstanding.

[kate]

Yes, you are right, many idealized domains have angler=0 everywhere. Ditto grids aligned with latitude, longitude lines. ROMS grids ought to be orthogonal, so there's only one angler to describe the local rotation.

[robjenkinsiii]

Great! That information has prompted me to find another, simpler way to describe my stress analytically.

I'm going to attempt to pass along dx & dy from ana_grid.h to ana_smflux.h as we discussed was possible for other arrays above and then I'll finally have been able to describe the stress in question.

Beyond that, I think I'm pretty close to my goal here.

I'm a little wary of adding a new calculation variable to ana_smflux. Not sure if it should be defined in both routines in that file if it's only being used to calculate stresses.

Finally I think I'll have to change the boundary conditions to 'Rad' in ocean_upwelling.in if I'm to see the physics I'm expecting for. Any trip ups there? I see it's probably best to set the VolCons terms to True for mass conservation, as well as the ad_VolCons terms.

[robjenkinsiii]

Minor update:

I managed to get the build to go through with no errors.

Execution blew up on the first time step.

I'm going to read the forum for others with this problem but wanted to post that update here.

No values jump out at me as being to large or too small in anything I've changed.

Could this be because I am starting the model with forcing in the first time-step with no "ramp-up?"


EDIT:

BLOWUP.rtf

The output from this attempt to run

(23.73 KiB) Downloaded 383 times

[kate]

I'm a little wary of adding a new calculation variable to ana_smflux. Not sure if it should be defined in both routines in that file if it's only being used to calculate stresses.

You need to add an argument to a subroutine call both to the place a routine is called and where the routine is written. The number and type of arguments have to match.

Finally I think I'll have to change the boundary conditions to 'Rad' in ocean_upwelling.in if I'm to see the physics I'm expecting for. Any trip ups there? I see it's probably best to set the VolCons terms to True for mass conservation, as well as the ad_VolCons terms.

You don't need ad_VolCons unless you are using the adjoint. I have never used VolCons either. I don't have recent experience with 'Rad' alone - I use 'RadNud' for the 3D fields, Chapman/Flather for the barotropic mode. Setting analytic boundary conditions of zero to nudge to might work for you, depending.

[kate]

Execution blew up on the first time step.

Getting 'exit_flag of 1' is coming from diag.F. I have made my diag.F more verbose in the case of blowing up, telling where and what value. Actually, it did tell you that:

Code: Select all

           (20,38,16)  9.097913E-01  1.731204E+00  0.000000E+00  1.077320E+02

I suggest you try a shorter timestep as the first step.

[robjenkinsiii]

In your speculation or experience, do poorly defined boundary conditions often cause blowing up errors?

Also, I made a slight edit, allowing the stress to increase gradually from zero to full stress over a two hour time period. Now it blows up at the 4th time step( about 20 minutes in to the full 12 hours). Not necessarily progress. Might not have defined this cleanly.

Also, doesn't rule out possible errors from other sources.

Thanks again for the tip, I'll try adjusting B.C.'s as I move ahead.

[kate]

Yes, but not in the first timestep. Is this point on the boundary? (20,38,16)

[robjenkinsiii]

Sorry, didn't see your second post.

Time step was another concern that I had. Is it possible to set a dynamic time step only to the beginning? I thought time step was fixed in ocean_upwelling.in

[robjenkinsiii]

(20, 38, 16) would be a point on the vertical boundary only 16, is the number of vertical levels. This however wouldn't be a boundary that the BC are concerned with.

[robjenkinsiii]

That is possibly insightful....

Adjusting the time step down get's me an extra step or two before blowing up. I can try adjusting it lower than 200s.

However My grid is (41,80,16) (arbitrary; same as upwelling case) so the point at which it is blowing up (20,38,16) is near the center of the grid.

My stress is radially defined and should be zero at the center, and very small near the center. I'll have to review also the form of my stress equation.

[kate]

Just yesterday I set up an idealized stress that might be similar enough for you to find it interesting. See attached, but note that I want it to be wider in the i-direction than in the j-direction. Making it symmetric is an exercise for the reader.

Yes, problems at the surface boundary can lead to trouble in the first step.

Edit: I just realized my reasons for adding the assymetry are invalid. Onward.