# Ocean Modeling Discussion

ROMS/TOMS

Search for:
 It is currently Mon Feb 18, 2019 11:00 am

 Page 1 of 1 [ 2 posts ]
 All times are UTC
Author Message
 Post subject: Problem with nearshore_mellor05.hPosted: Fri Nov 05, 2010 12:27 pm

Joined: Fri Sep 17, 2004 2:22 pm
Posts: 74
Location: Institut Rudjer Boskovic
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)
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.

Top

 Post subject: Re: Problem with nearshore_mellor05.hPosted: Fri Mar 25, 2011 7:19 pm

Joined: Wed Dec 31, 2003 6:16 pm
Posts: 758
Location: USGS, USA
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.
- 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

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 2 posts ]

 All times are UTC

#### Who is online

Users browsing this forum: No registered users and 1 guest

 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forumYou cannot post attachments in this forum

Search for:
 Jump to:  Select a forum ------------------ News, Events & Job Opportunities    Meetings/Workshops    Job Opportunities    Ocean News ROMS/TOMS    ROMS Adjoint    ROMS Benchmarks    ROMS Bugs    ROMS Discussion    ROMS Documentation    ROMS Ecosystem    ROMS FAQ    ROMS Ice    ROMS Information    ROMS Installation    ROMS Known Problems    ROMS Messages    ROMS Problems    ROMS Releases    ROMS Results    ROMS Sediment    ROMS Source    ROMS Tips    ROMS Tools and Techniques    ROMS Trivia    ROMS Usage    ROMS Webinar    ROMS Wish List ROMS/TOMS Applications    User Applications    Adriatic Sea Ocean Modeling    Ocean Modeling FAQ