I found a number of problems with the implementation of the nearshore formulation of the radiation stress.

The problem can be detected very easily since there is both a formulation of the stress for barotropic and baroclinic in

nearshore_mellor05.h. Hence if one removes the

#ifdef SOLVE3D ... #else ... #endif and rename the variables then one sees that the

*2d variables are completely differently. This should not be the case since the formula of Mellor 2003 when vertically integrated gives the Longuet-Higgins Stewart 1962 formulation.

Putting everything together takes time. Here is the result of my investigation:

Issue 1: The fact that Mellor 2003, when vertically integrated gives the classical solution relies on the following integral identities:

**Code:**

int_0^D FCS FCC dz = 1/(2k) + D/sh(2kD)

int_0^D FSS FCS dz = 1/(2k)

int_0^D FCS FCC - FSS FCS dz = D/sh(kD)

Hence if the code reproduces the same result it means that it computes those integrals correctly, which is certainly hard. The problem is that the effective computation of

SINH(kD) and friends requires a limiter on the value of

kD. Currently, the limitation is at 5. The problem is that in the vertical integration

Hz is not limited. This can result in discrepancy in the numerically integrated quantity.

This is not a small problem. The situation below:

**Code:**

ThetaS=7 ThetaB=0.1 TCLINE=0 Vstretching=2 Vtransform=2

dep=1160m Lwave=10 Hwave=1.01

Sxx(numerical integral)=18 Sxx(2D formula)=0.12

which is realistic shows a 144 times larger value of

Sxx.

Solution to the problem: First of all the value of

kDmax=5 is really extremely conservative. ROMS works in double precision, and

SINH(710) is still finite while

SINH(711) is infinite. Hence setting

kDmax=300 sounds reasonable to me as far as I know.

But setting this does not solve completely the problem. The amplification factor is still 2.5 for

kDmax=300. The reason is that the integrals of

FCS FCC and others are extremely singular. Thus I propose to replace the value of

FCS FCC by their average values over the interval

[z_w(i,j,k-1) z_w(i,j,k)]. This is possible since those functions have explicit expressions that can be explicitly integrated easily.

Issue 2: The term

rustr3d is multiplied by the quantity

**Code:**

cff1=0.25_r8*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))

at line 846 of the code. Then when vertically averaging is done

rustr2d is again multiplied by this quantity on line 920. For the barotropic code there is only one multiplication by this quantity; it is on line 1033.

Proposed solution: remove the multiplication by this quantity.

Issue 3: after that the vertical integral still diverge. What I think is the problem is that

Sxx is multiplied by

Hz two times. Once in line 506 and then in line 909. Right now I cannot think why this should be.

Proposed solution: remove the multiplication by cff5 and subsequent averaging. If one does that in line 909 we do

**Code:**

rustr2d(i,j)=cff5*rustr3d(i,j,1)

replaced by

**Code:**

rustr2d(i,j)=rustr3d(i,j,1)

On line 914 we do

**Code:**

rustr2d(i,j)=rustr2d(i,j)+cff5*rustr3d(i,j,k)

replaced by

**Code:**

rustr2d(i,j)=rustr2d(i,j)+rustr3d(i,j,k)

and on line 920 we simply remove

**Code:**

rustr2d(i,j)=rustr2d(i,j)*cff3*cff4

(remember the unnecessary multiplication in issue 2)

Issue 4: After that the term

rustr2d still diverge. One reason is that the integral is that the multiplication by

Hz are not done in the same way in barotropic and baroclinic for

UFe.

Proposed solution: replace

**Code:**

UFe(i,j)=0.25_r8*(Hz(i,j ,k)+Hz(i-1,j,k)+ &

& Hz(i,j-1,k)+Hz(i-1,j-1,k))* &

& Sxyl_psi(i,j)

by

**Code:**

UFe(i,j)=0.25_r8*(Hz(i,j,k)*Sxyl(i,j)+ &

& Hz(i-1,j,k)*Sxyl(i-1,j)+ &

& Hz(i-1,j-1,k)*Sxyl(i-1,j-1)+ &

& Hz(i,j-1,k)*Sxyl(i,j-1))

Issue 5: after that one last step is needed to get the same result: Add the

#ifdef CURVGRID ... #endif corrections to

rustr2d since we want to run in curvilinear grid right? That correction was obvious. Another is that the variable

Sxy_psi is used in barotropic mode while not declared. The solution is simply to replace it by

Sxyl_psi. This is a compile error which shows that maybe some regression tests should be made on ROMS. After that

rustr2d give the same value. The reason is that the vertical correction terms are always found to be very small at least in the case that I considered.

Issue 6: Now it remains to correct for ubar_stokes. The key formula is

**Code:**

1/D int_0^D cosh(2kz)/sinh(2kD)dz = 1/(2kD)

This is simpler than for

Sxx, but the formula as written is problematic:

**Code:**

u_stokes(i,j,k)=cff2* &

& (wavenx(i-1,j)+wavenx(i,j))/ &

& (wavec (i-1,j)+wavec (i,j))* &

& COSH(cff3*fac2)*0.5_r8* &

& (o2sinh(i-1,j)+o2sinh(i,j))

The problem is that the cosine hyperbolic is taken for the middle point, while the sinh is taken at the points and then averaged which is not consistent. The problem is that we always have the inequality

**Code:**

1/sinh((a+b)/2) <= 1/2sinh(a) + 1/2sinh(b)

and that the difference can be very significant. See below an example (which is still in the case

kDmax=5, i.e. with the conservative estimate)

**Code:**

1/sinh(7.5) = 0.0011 and (1/sinh(5)+1/sinh(10))/2 = 0.0068

i.e. the error is of a factor 6. If one replaces

**Code:**

*0.5_r8*(o2sinh(i-1,j)+o2sinh(i,j))

by

**Code:**

/SINH(cff3)

then the error is smaller but can still be of a factor of 10 always in overestimating magnitudes of barotropic velocities. The solution then is to replace u_stokes by the mean value over the interval

[z_w(i,j,k-1), z_w(i,j,k)] just as we did for

Sxx.

Miscellaneous issues: In the computation of depth, the thermocline is not used. It is as if TCLINE=0. I.e. in the existing code we have

**Code:**

fac2=1.0_r8+SCALARS(ng)%Cs_r(k)

COSH(kD(i,j)*fac2)

instead of the thermocline preserving

**Code:**

COSH(waven(i,j)*z_r(i,j,k))

But I do agree that using

z_r and still having the limiter will considerably complexify the code.

There is a code section

#ifdef GEO_ROTATION ... #endif which is actually never triggered (top of file is

#undef GEO_ROTATION). What was the idea or the problem there?

Modified file is in attachment. All comments welcomed.