Problem with nearshore_mellor05.h

[mathieu]

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: Select all

   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: Select all

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: Select all

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: Select all

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

replaced by

Code: Select all

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

On line 914 we do

Code: Select all

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

replaced by

Code: Select all

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

and on line 920 we simply remove

Code: Select all

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: Select all

  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: Select all

  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: Select all

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: Select all

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: Select all

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: Select all

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: Select all

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

by

Code: Select all

/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: Select all

fac2=1.0_r8+SCALARS(ng)%Cs_r(k)
COSH(kD(i,j)*fac2)

instead of the thermocline preserving

Code: Select all

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.

[jcwarner]

Mathieu-

Sorry for the looooong delay in response. This Mellor method has been problematic, but i wanted to provide some feedback so at least i could try to bring some closure to this post. We initially coded in the method as described in Mellor 03 paper, but then realized there were some errors. He published an update in the 05 paper. That is the version coded here. Since then, there was a '08 update, which we coded (not distributed on Rutgers). We made some modifications, are have been using that for a while. Now there is a 2011 submission, suggesting another update. We looked at that but it does not seem to completely resolve all the issues.

So for now, i have made changes mostly as you suggested, to modify this Mellor 05 method. Instead of going thru the gory details, the summary of the changes include:
- change kdmax = 300.
- replace Cs_r and CS_w with appropriate z_w and z_r computations.
- Use vertical w points and avg them in computations of FCC, FCS, and FSS, and stokes vels.
- did not remove the multiplication for pm pn, etc. As this is required to get the correct units for output purposes. need to have rustr2d 3d etc as m2/s2, because for output they are *rho0 and this makes them in Pa.
- modify the UFe computation to use Sxyl_psi info.
- added curvgrid for 2D.
- correct 2d computation in 3D simulation so it is the same as a 2D only computation.

Please try this attached file. If this works for you we can push it for distribution on the Rutgers trunk.

-john