Analytical Tide defined in O.B.C. in COAWSTv3.1

[CeeCee.X]

Hi,All,

I am now using the COAWST v3.1. I tried to input an analytical M2 tide through the east open boundary conditions for an ideal test case. My domain is rectangular. What I did is add the following sentence into ana_fsobc.h and ana_m2obc.h:

# elif defined IDEAL
fac=1.0_r8 ! zeta0
omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8) ! M2 Tide period
# ifdef EAST_FSOBC
IF (EASTERN_EDGE) THEN
val=1.0_r8
phase=(0.0_r8)*deg2rad
DO j=JstrR,JendR
BOUNDARY(ng)%zeta_east(j)=fac*val*COS(omega-phase)
END DO
END IF
# endif

Meanwhile, I defined my boundary conditions in ocean.in as following:
LBC(isFsur) == Clo Per Cla Per ! free-surface
LBC(isUbar) == Clo Per Red Per ! 2D U-momentum
LBC(isVbar) == Clo Per Red Per ! 2D V-momentum
LBC(isUvel) == Clo Per Red Per ! 3D U-momentum
LBC(isVvel) == Clo Per Red Per ! 3D V-momentum
LBC(isMtke) == Clo Gra Gra Gra ! mixing TKE

Therefore, in my cpp option I didn't include ADD_M2OBC and ADD_FSOBC.

However, the results shows there is no oscillation in the east boundary throughout my simulation (I ran for 5 days with time step equal to 60 second). I checked the analytical.f90 and find the following, which seems incorrectly to me:
!-----------------------------------------------------------------------
! Free-surface open boundary conditions.
!-----------------------------------------------------------------------
!
fac=1.0_r8 ! zeta0
omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8) ! M2 Tide period
RETURN
END SUBROUTINE ana_fsobc_tile
SUBROUTINE ana_initial (ng, tile, model)

Now I have to manually modify the analytical.f90 to add the open boundary condition. I was wondering since the new version of ROMS have moved into input file, how shall we define in the cpp option for it to recognize the open boundary condition. Are those old cpp boundary condition still useful? Or shall I just simply remove the if sentences: # ifdef EAST_FSOBC; IF (EASTERN_EDGE) THEN to make it run?

Thanks

[jcwarner]

suggest you use Inlet_test as an example, although that is for a northern 'tide'. If you look in Projects/Inlet_test/Coupled/ana_fsobc you will see a slightly different logic to impose the free surface
#if defined INLET_TEST
IF (LBC(inorth,isFsur,ng)%acquire.and. &
& DOMAIN(ng)%Northern_Edge(tile)) THEN
cff=-1.0_r8*sin(2.0_r8*pi*time(ng)/(12.0_r8*3600.0_r8))
DO i=IstrR,IendR
BOUNDARY(ng)%zeta_north(i)=cff
END DO
END IF
#else ...

so you want to use something like
IF (LBC(ieast,isFsur,ng)%acquire.and. &
& DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=JstrT,JendT
BOUNDARY(ng)%zeta_east(j)= .....
END DO
END IF

give that a try.
-john

[LiuZHQ]

Wondering what version of ROMS you're using. You defined the OBC for three-dimensional velocity as Red, which supposedly represented a reduced-physics OBC. I checked the ROMS code (Version 645 2013-01-22 23:21:54Z), and didn't find that option.
If the actual OBC is a closed one, due to, possibly, wrong typing, and through the coupling process between the external and internal mode, the tidal forcing can not be well implemented, then.
For your reference~~

[kate]

The new COAWST is using the new boundary formulation with LBC. You therefore should get rid of all cppdefs like:

# ifdef EAST_FSOBC

[CeeCee.X]

Thank you, Dr.John.Now it works.

I have two additional question related to this ideal simulation:
I found a header file named as ana_m2obc.h, which looks to me that it also deals with O.B.C. Since in my test case, the M2 component is imposed on the EAST boundary, I thought it is necessary to call this header file as well. But it turned out in the inlet test case, it is not used. I was wondering under what condition does this m2obc header file to be called?

[CeeCee.X]

Thanks, Dr.Kate.

[kate]

You would need to provide an m2obc (either via a file or ana_m2obc) if you use boundary conditions that need a value for ubar/vbar on the boundary. You have picked reduced, which doesn't need ubar/vbar. We have had good luck providing tidal flow as part of the boundary condition and therefore use Shc for LBC for ubar/bvar.

[CeeCee.X]

Hi, Dear All,

During my simulation, I encountered a new problem when trying to add constant winds through analytical function in ROMS. I ran my simulations for five days with 1 min time step, however, after I define svstr=10 in the enana_smflux.h, I encountered such problem:

NL ROMS/TOMS: started time-stepping: (Grid: 01 TimeSteps: 00000001 - 00007200)


STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME
C => (i,j,k) Cu Cv Cw Max Speed

0 0 00:00:00 0.000000E+00 4.861661E+03 4.861661E+03 6.219716E+13
(000,000,0) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
DEF_HIS - creating history file, Grid 01: ocean_his.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
1 0 00:01:00 3.302214E+01 4.861662E+03 4.894684E+03 6.219715E+13
(001,001,8) 3.593298E-03 2.696580E+01 0.000000E+00 8.988599E+02
DIAG speed trouble 898.859924658695
DIAG speed ijk 898.859908697992 1 1 8

Blowing-up: Saving latest model state into RESTART file

Why the ROMS cannot receive constant wind stress?

[kate]

A more typical surface stress is of the order 0.1 Newton/m^2.

[CeeCee.X]

Thank you, Kate. I adjusted svstr(i,j)=0.1_r8, however, it only lasted 2 steps more.

STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME
C => (i,j,k) Cu Cv Cw Max Speed

0 0 00:00:00 0.000000E+00 4.861661E+03 4.861661E+03 6.219716E+13
(000,000,0) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
DEF_HIS - creating history file, Grid 01: ocean_his.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
1 0 00:01:00 3.301685E-03 4.861661E+03 4.861665E+03 6.219716E+13
(001,001,8) 2.244533E-21 2.689694E-01 0.000000E+00 8.965645E+00
2 0 00:02:00 1.334396E-02 4.861661E+03 4.861674E+03 6.219715E+13
(001,001,8) 4.383571E-21 5.354244E-01 1.706769E-21 1.784748E+01
3 0 00:03:00 1.793436E-02 4.861659E+03 4.861677E+03 6.219713E+13
(001,001,8) 5.978067E-21 7.455880E-01 2.498850E-23 2.485293E+01
DIAG speed trouble 24.8529330651598
DIAG speed ijk 24.8529330651598 1 1 8

Blowing-up: Saving latest model state into RESTART file

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

Elapsed CPU time (seconds):

Node # 0 CPU: 2.967
Total: 2.967

[kate]

I wrote on the wiki about this stuff. At least you know where the trouble is: point 1,1,8.

[kate]

Also, if k=8 is at the surface, how deep is your domain? What is your vertical stretching like? How well resolved is the surface boundary layer? Applying a surface stress with N=8 for 5000 m of water is asking for trouble. The BODYFORCE option was created for this case, applying the wind forcing over the top few levels as a body force.

[CeeCee.X]

Hi, Kate,

Here I attached my ideal study domain. The point(1,1,8) is located at boundary of the shallow water. But I was wondering if grid cell (1,1) is ghost cell or not? Basically when I create the domain, I didn't mask any cells dry because it is a simple test for an assumption in our real case study.

I set number of vertical layers to be 8. Here are my basic grid dimension parameters:

Lm == 174 ! Number of I-direction INTERIOR RHO-points
Mm == 124 ! Number of J-direction INTERIOR RHO-points
N == 8 ! Number of vertical levels,feng-try-2D first

Can you provide me more details about "The BODYFORCE option" and "applying the wind forcing over the top few levels as a body force"?

Thanks.

[kate]

Nothing beats actually looking at the source code for some of these questions:

Code: Select all

              u2v2=u(i  ,j,k,nstp)*u(i  ,j,k,nstp)+                     &   
     &             u(i+1,j,k,nstp)*u(i+1,j,k,nstp)+                     &   
     &             v(i,j  ,k,nstp)*v(i,j  ,k,nstp)+                     &   
     &             v(i,j+1,k,nstp)*v(i,j+1,k,nstp)

For rho point (1,1), it is getting input from both interior and boundary points. Because the west is closed and the south is periodic, I wouldn't expect trouble to come from the boundary points. You haven't driven 5 m of water away in 4 timesteps, have you? What is zeta there?

Depending on your vertical mixing, you might be resolving the surface boundary layer in your shallowest water, but probably not in the deep water. Your jump from 50 m to 1000 m depth is scary and will lead to spurious currents along that wall.

The BODYFORCE option applies the forcing over some uniform number of k levels, surface forcing from k=N down to k=levsfrc, bottom forcing from k=1 up to levbfrc. Both LEVSFRC and LEVBFRC are set in the ocean.in.

[CeeCee.X]

zeta along the west boundary is no more than 1.6e-04 m. What happened is the V-component is too high (over 35 m/s) along the west boundary. If Coriolis effect is not considered. Everything went fine. However, if f is not zero, the cell gets blow up.

[CeeCee.X]

And now I am trying BODYFORCE. In order to do that, I was wondering if the following steps are correct:
1) define BODYFORCE in cpp option
2) Set LEVSFRC = 8 in ocean.in

Besides these two steps, I was wondering if I need to generate a new netcdf file for the wind? Also, I was confused about the index of vertical layers. As we have discussed before, N =8 is the surface layer. Therefore, what does it mean by saying "from k=N down to k=levsfrc" ?

[kate]

Yes, that's right for steps 1 and 2.

No need to redo the winds.

For the loop, the code tells the story:

Code: Select all

      DO k=N(ng),levsfrc(ng),-1

Another loop is the other way around:

Code: Select all

      DO k=levsfrc(ng),N(ng)

If you grep the code for BODYFORCE, you can find this stuff in rhs3d.F.

[CeeCee.X]

I got you. Thank you, Kate. Let me play with it and see how thing's going.

[CeeCee.X]

Hi,
I have tried with BODYFORCE, but there is no much improvement. The model still blows up either at western boundary or at the junction between shallow and deep water, due to the speed trouble. Now I am testing svstr = 0.5N/m^2. And the current velocity along y_rho direction often gets into a very large value( > 40 m/s) either at bottom or at surface.

I could eliminate the surface stress near western boundary, but that will be against my initial goal. Since I would like to see how the free surface elevation in the shallow water behaves driven by the wind with a scale of hurricane-status. Also, for the steep slope, if I smooth it, it will be bias from the reality which I am simplifying.

So I was wondering if the above issue is due to sigma grid. If so, can 2D simulation avoid the above issue? Also, can ROMS do hybrid grid or z gird? If so, how shall I implement it?

I also appreciated any ideals to improve my simulation.

Thanks.

[kate]

If you do a 2D simulation, you will avoid the pressure gradient errors associated with the jump in bathymetry. I'm not sure you get rid of all the errors associated with it, though. I know the real world has some pretty steep bathymetry - that drives some of our desire for increased horizontal resolution. You can replace the jump with say a tanh function and increase the resolution there - but you want to change resolution smoothly too. I know that's not the problem you want to be solving, but it's the problem the model was designed to handle.

The BODYFORCE option was designed for domains which are deep everywhere and with insufficient vertical resolution. Our more recent runs have just gone for more vertical resolution (N=50) and I haven't used BODYFORCE in over a decade.

I'm not a hurricane expert, but is that really the surface stress from them? That's a very large surface sress. I've heard that in nature, increasing the winds doesn't lead to ever increasing surface stress and that the default ROMS bulk flux routine overestimates stress at high winds.

[kate]

OK, checking my Bering Sea storm case, the wind stresses get up to 1 N/m^2. It is running fine. I've got hmin=5, N=50, tides, WET_DRY, and enough bathymetry smoothing to bring the rx0 down to 0.3 (search on Beckmann in the ROMS stdout).

[CeeCee.X]

Good to know all of these. Thanks a lot, Kate.
I masked some cells that blowing up to be dry and now it works fine. I am also trying to make a new grid just as you said making a transition from the shallow to deep. But what I have found that is the cell during this transition portion is very weird, they will have 2km in eta_rho but 180m in x_rho direction in order to have a smooth slope of each cell grid (delta_depth = 0.9m). I am not sure that if this is a good idea. Also, I was wondering since I make the cells in the transition portion to be finer, do I need to apply the same fine resolution for my shallow water area?
Thanks.

[kate]

you want to change resolution smoothly too.

That rx0 is a measure of how much your depth changes from one grid to the next. One could also create a metric for how much your dx changes from one grid to the next, plus a metric for grid aspect ratio. I don't know of any careful studies of either of those.