Horizontal surface pressure gradient forcing in 2D and 3D

[m.hadfield]

I have created this post in order to discuss the inclusion of forcing by horizontal surface (atmospheric) pressure gradients in ROMS. As I understand it, having read this thread...

viewtopic.php?t=324

..., run a test simulation and perused the code, there is a preprocessor option called ATM_PRESS, described as "use to impose atmospheric pressure onto sea surface". If this is activated and if SOLVE3D is also activated, then the model requires atmospheric pressure to be specified, either analytically or from the forcing file. However as far as I can see--and as my test simulation indicates--the air pressure has no dynamic effect.

On the thread, the code required to feed the surface pressure into the model is discussed: in 3D, apply it as an upper boundary condition on hydrostatic pressure; in 2D, add an atmospheric-pressure gradient term to the surface-height gradient term.

I propose that the code be modified to support the ATM_PRESS option in 2D as well as 3D, and to include the air pressure term in the pressure-gradient in 2D and 3D. I can do most of the coding, I think, but may need some support from Hernan on the adjoint.

Mark

[m.hadfield]

Attached is a ZIP file containing source code modifications (relative to the current trunk) to implement atmospheric pressure gradient forcing in 2D and 3D.

Variable Pair is read in (or calculated by the analytic routine) if ATM_PRESS is defined, whether or not SOLVE3D is defined. This has been done for the Adjoint, Nonlinear, Representer and Tangent models.

For the Nonlinear code only (so far) I have added pressure or pressure-gradient terms in the momentum equation for 2D and for 3D (all baroclinic options). The multiplying factor for the pressure in these calculations is either 100/rho0 or 100/g, where the 100 converts millibars to Pascals and the 1/rho0 or 1/g depends on the dimensions for pressure in the equations, as determined by my reading of the code plus a bit of trial and error.

I set up a simple test: water in a basin is exposed to a sinusoidal perturbation in surface pressure, with amplitude 20 mb; no other forcing. I tested the 2D option and all the 3D baroclinic options. Results are almost identical: the 20 mb pressure perturbation generates a 0.198 m zeta perturbation at equilibrium. However with the PJ_GRADP, PJ_GRADPQ2 and PJ_GRADPQ4 options, the zeta perturbation is 0.195 m, i.e. 1.5% smaller.

I will have a look at the Tangent, Representer and Adjoint models shortly. Or maybe not.

[m.hadfield]

Thanks, Hernan, for merging my changes into the trunk:

https://www.myroms.org/projects/src/ticket/303

Note that the modifications to the 2D code are essentially those posted by rongzr in this thread:

viewtopic.php?t=324

and the 3D mods come from Hernan originally (I think), via Paul Budgell and Kate Hedstrom.

If anyone can suggest why the 3D PJ_GRAD* variants give a slightly smaller equilibrium response than the 2D and other 3D variants, I'd love to hear it.

[jcwarner]

Did anyone test these !!!!!!????

prsgrd31.h:
#ifdef ATM_PRESS
phie(i)=phix(i)+fac*(Pair(i,j)-Pair(i,j-1))
#endif
Probably a bug !! looks like it should be phie = phie + ....

Maybe we should test these ???
This should get corrected, and a new roms v3.2.1 tar file made.
IT WOULD BE VERY NICE TO SPEND SOME TIME TO TEST THE CODE, MAKE CORRECTIONS, AND THEN RELEASE A STABLE VERSION. PERHAPS WE COULD START HERE. LETS STOP ADDING NEW STUFF FOR A WHILE, TEST WHAT IS HERE, CORRECT IT, AND THEN RELEASE A STABLE VERSION.
PERHAPS YOU ALL COULD DISCUSS THIS DOWN UNDER ??
(sorry i wont be there)

[arango]

Yes, good catch Typos like this can happen. I only tested the prsgrd32.h, the default algorithm. It was my understanding that Mark tested all of them. Fortunately, prsgrd31.h is deprecated algorithm and rarely used nowadays. It is not the first time that thing like this can happen and sometimes idealized test cases do not expose the problem.

[m.hadfield]

...It was my understanding that Mark tested all of them. ... It is not the first time that thing like this can happen and sometimes idealized test cases do not expose the problem.

Especially test cases in which surface pressure varies only in the xi direction!

[m.hadfield]

If anyone can suggest why the 3D PJ_GRAD* variants give a slightly smaller equilibrium response than the 2D and other 3D variants, I'd love to hear it.

This difference is not confined to forcing by atmospheric pressure; it also applies to forcing by surface stress. For a test, I set up a closed box, everything uniform, with walls at the eastern and western boundaries, no Coriolis force, eastward surface stress at 1.E-4 m2/s2 specified in ana_smflux.h. The equilibrium response has a linear gradient in zeta, with the pressure-gradient due to the slope balancing the stress. That slope is 2.5% smaller in 3D runs with PJ_GRAD* than it is for 2D and other 3D variants. The same difference (to within 0.1%) is seen in my test simulations of the equilibrium zeta response to surface pressure variations, as I mentioned above.

Hmm, 2.5% sounds suspiciously close to the diffence between typical oceanic and freshwater densities (1025 vs 1000 kg/m3).

PS: I had RHO0 set to 1025 in the input files.

[jcwarner]

why is "OneAtm" removed from the Pair in prsgrd32??
I think this would not matter.

[kate]

I too would like to put in a vote for a testing framework. Georgina just got bitten by a change to a "magic" number Natt (number of attributes) which was recently changed, but not in my routines. And why isn't Natt a global if it matters globally what its value is?

[m.hadfield]

why is "OneAtm" removed from the Pair in prsgrd32?? I think this would not matter.

You'd think so, wouldn't you, but with one of the pressure-gradient options (I forget which) an instability developed that went away when the OneAtm term was introduced. Something to do with truncation error??

[m.hadfield]

You'd think so, wouldn't you, but with one of the pressure-gradient options (I forget which) an instability developed that went away when the OneAtm term was introduced.

All three PJ_GRADP* options develop this instability. This plus their slightly different response to surface forcing (pressure and stress) suggests they should be used with caution until the reasons are understood.

Does anyone use them?

[m.hadfield]

If you use ATM_PRESS with M2FLATHER boundary conditions, you need to add a term to allow for the surface atmosphere pressure in calculating the difference between the model SSH and the boundary SSH. Something like this:

Code: Select all

#  ifdef ATM_PRESS
      OneAtm=1013.25_r8                  ! 1 atm = 1013.25 mb
      pfac=100.0_r8/(g*rho0)
#  endif

Code: Select all

#  ifdef ATM_PRESS
          ubar(Istr,j,kout)=bry_val-                                    &
     &                      Cx*(0.5_r8*(zeta(Istr-1,j,know)+            &
     &                                  zeta(Istr  ,j,know)+		&
     &                                pfac*(FORCES(ng)%Pair(Istr-1,j)+  &
     &                                      FORCES(ng)%Pair(Istr  ,j)-  &
     &                                      2.0_r8*OneAtm))-            &
     &                          BOUNDARY(ng)%zeta_west(j))
#  else
          ubar(Istr,j,kout)=bry_val-                                    &
     &                      Cx*(0.5_r8*(zeta(Istr-1,j,know)+            &
     &                                  zeta(Istr  ,j,know))-           &
     &                          BOUNDARY(ng)%zeta_west(j))
#  endif

This term will be needed unless the boundary SSH was calculated by a model in which surface pressure forcing was included.

Neglect of this term (or inclusion when it's not appropriate) can lead to surprisingly strong spurious flows and this can result in erroneous statements being made at ROMS workshops.

[arango]

Okay, I know that this is an old post, but I have been looking at the 3D and 2D pressure gradient terms and the various algorithms available. I implemented the astronomical tide generation forces (TGF) in terms of the equilibrium tide (ζeq), defined as the shape that the sea surface (m) would assume if it were motionless and in equilibrium with the TGF on a fluid planet (see trac ticket).

The first thing that I noticed is a typo/bug in prsgrd32.h:

Code: Select all

          P(i,j,N(ng))=GRho0*z_w(i,j,N(ng))+                            &
#ifdef ATM_PRESS
     &                 fac*(Pair(i,j)-OneAtm)+                          &
#endif
     &                 GRho*(rho(i,j,N(ng))+cff2)*                      &
     &                 (z_w(i,j,N(ng))-z_r(i,j,N(ng)))

We are computing here kinematic pressure (P/rho0), so P has units of m2/s2. This algorithm is based on the density Jacobian, and we compute kinematic pressure, vertically integrate, and evaluate its gradient. The issue here is the first term that includes the surface elevation, z_w(i,j,N(ng)), needs to be scaled by g instead of Grho0 (g/rho0). That term is only affecting the value of g. It is around a 0.97 smaller. This bug has not yielded significant effects on my testing solutions. Recall that the value of P is not as important as its horizontal grandient.

However, I have lots of issues when I implemented the TGF term in the pressure gradient. Recall that the TGF term is defined as: -g [∂(ζ - ζeq)]/∂x . Therefore, the correct code is:

Code: Select all

#ifdef ATM_PRESS
      OneAtm=1013.25_r8                  ! 1 atm = 1013.25 mb
      fac=100.0_r8/rho0
#endif
      ...
 
          P(i,j,N(ng))=g*z_w(i,j,N(ng))+                                &
#ifdef ATM_PRESS
     &                 fac*(Pair(i,j)-OneAtm)+                          &
#endif
     &                 GRho*(rho(i,j,N(ng))+cff2)*                      &
     &                 (z_w(i,j,N(ng))-z_r(i,j,N(ng)))
#ifdef TIDE_GENERATING_FORCES
          P(i,j,N(ng))=P(i,j,N(ng))-g*eq_tide(i,j)
#endif
        END DO

Below are two figures that show the sea surface elevation for the center of the Gulf of Maine (top) and at the entrance of the Bay of Fundi (bottom). ROMS is configured with boundary tidal forcing, tide generating force, and inverted barometer effect in the pressure gradient.




As expected, the TGF makes very small corrections to the boundary tidal forcing. Users will need to activate both tidal forcings (SSH_TIDES, UV_TIDES, TIDE_GENERATING_FORCES). The ATM_PRESS is activated for all solutions.

Now, I am wondering if the stability issues that Mark was talking about when implementing the inverted barometer (IB) effect term needed (Pair-OneAtm). See the ATM_PRESS option. Maybe the issue is due to the bug of the scaling of the sea surface term.

Recall that the inverted barometer effect is usually defined as: -g [∂(ζ - ζib)]/∂x where ζib = P'/(rho0*g). See Wunsch and Stammer (1997) paper: Atmospheric Loading and the Oceanic Inverted Barometer Effect.

Therefore, I think that there is also a bug in the barotropic inverted barometer effect. In step2d_LF_AM3.h we have:

Code: Select all

#if !defined SOLVE3D && defined ATM_PRESS
      fac3=0.5_r8*100.0_r8/rho0
#endif
      ...

#if defined ATM_PRESS && !defined SOLVE3D
          rhs_ubar(i,j)=rhs_ubar(i,j)+                                  &
     &                  fac3*on_u(i,j)*                                 &
     &                  (h(i-1,j)+h(i,j)+                               &
     &                   gzeta(i-1,j)+gzeta(i,j))*                      &
     &                  (Pair(i-1,j)-Pair(i,j))
#endif

It seems to me that if we follow Wunsch and Stammer (1997), we need to have instead:

Code: Select all

#if defined ATM_PRESS && !defined SOLVE3D
          rhs_ubar(i,j)=rhs_ubar(i,j)-                                  &
     &                  fac3*on_u(i,j)*                                 &
     &                  (h(i-1,j)+h(i,j)+                               &
     &                   gzeta(i-1,j)+gzeta(i,j))*                      &
     &                  (Pair(i,j)-Pair(i-1,j))
#endif

We need the minus sign and flip the order of the spatial derivative of Pair. Similarly, we need to do the same for the other component rhs_vbar (m4/s2). We haven't noticed this because this part of the code is only activated for shallow-water applications (SOLVE3D is off). We need to design an application for testing this. Say a shallow-water storm surge application. Any volunteers? I am swamped a the moment, and it has very low priority on my long list.

[jcwarner]

I think the original code was correct. P is only an internal variable for this routine. Maybe you want to re-define it for the tides, but the way it was originally seemed correct. It was:

Code: Select all

         P(i,j,N(ng))=GRho0*z_w(i,j,N(ng))+                         &
     &                 GRho*(rho(i,j,N(ng))+cff2)*                      &
     &                 (z_w(i,j,N(ng))-z_r(i,j,N(ng)))

dimensionally this is:

Code: Select all

P = GRho0*z_w   +  GRho*( rho+cff2 )* (z_w - z_r)
  = 1000* g /rho0 * m   + g/rho0 (  rho + rho) * m

(Recall that the '1000' is really a rho)

  = g*m  +  g*m
  = g*m

so this is grav * m  = m^s / s^2

The units need to be this way to compute ru and rv below:

Code: Select all

            ru(i,j,k,nrhs)=0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k))*on_u(i,j)*
     &                     (P(i-1,j,k)-P(i,j,k)-

.....

This is cumbersome, but just look at the first 2 lines we see

Code: Select all

ru = Hz* on_u * (P ...
   = m *  m *  m^2/s^2   = m^4 / s^2

and that is what we need for ru.
-j

[arango]

Nope, I am taking about the first term: GRho0*z_w(i,j,N(ng)). It needs to be g*z_w(i,j,N(ng)), which also gives m2/s2 for kinematic pressure. I am sure that this term needs to be multiplied by gravity only and not scaled by density like it is done the step2d_LF_AM3.h for the barotropic pressure gradient using the non-dimensional vertically-averaged density (rhoA) and density perturbation (rhoS). Recall, that gravity wave phenomena is resolved in step2d , and we always recommend to activate VAR_RHO_2D. Anyway, the error is small.

[jcwarner]

GRho0 * z_w =
1000.0_r8 * GRho * z_w =
1000.0_r8* g/rho0 * z_w

so this is g*z_w * 1000/rho0
and you want it to be
g*z_w

well, the difference is 1000/1025 = 0.97, as you said.
But it does not have the wrong units. The units are correct.
The 1000/rho0 is density over density and those units cancel.

I think the reason for the 1000/rho0 scaling is from Sasha's paper:
A Method for Computing Horizontal Pressure-Gradient Force in an
Oceanic Model with a Non-Aligned Vertical Coordinate
equation 1.7 it has

[arango]

My understanding is that the nondimensional density scale (see rhoA and rhoS in rhos_eos.F) enters only in the barotropic pressure gradient in step2d with VAR_RHO_2D. Notice that the paper's equation shows the density evaluated at the free surface, which is only obtained by vertically integrating that equation. In a 3D model application, we don't have state values at the sea surface. It is similar to the extrapolation problem that we have with SST and SSS. This density scale depends on space and time, so it affects the horizontal gradient. Then, that term is accounted for during the coupling of the 3D- and 2D-momentum equations. A constant (1000/rho0) everywhere doesn't have a gradient. I think that Sasha's usage of rhoA and rhoS are very solid numerically. Maybe he will reply to these postings. The only thing that is doing is to reduce the value of gravity, which is assumed to be constant in ROMS. However, on Earth, the value of gravity is a function of latitude. Anyway, it is all a game of tiny numbers. As I said, the impact is very negligible in all my tests.

Yes, I mistyped about the wrong unit in the above post.

[jpcrms]

Dear arango

I try to run 2D mode in ROMS version 4.1.
In step2d_LF_AM3.h, we have:

Code: Select all

!-----------------------------------------------------------------------
!  Compute pressure gradient terms.
!-----------------------------------------------------------------------
!
      cff1=0.5_r8*g
      cff2=1.0_r8/3.0_r8
#if !defined SOLVE3D && defined ATM_PRESS
      fac3=0.5_r8*100.0_r8/rho0
#endif
...
#if defined ATM_PRESS && !defined SOLVE3D
          rhs_ubar(i,j)=rhs_ubar(i,j)+                                  &
     &                  fac3*on_u(i,j)*                                 &
     &                  (h(i-1,j)+h(i,j)+                               &
     &                   gzeta(i-1,j)+gzeta(i,j))*                      &
     &                  (Pair(i,j)-Pair(i-1,j))
!#endif

The modified code is:

Code: Select all

#if defined ATM_PRESS && !defined SOLVE3D
          rhs_ubar(i,j)=rhs_ubar(i,j)-                                  &
     &                  fac3*on_u(i,j)*                                 &
     &                  (h(i-1,j)+h(i,j)+                               &
     &                   gzeta(i-1,j)+gzeta(i,j))*                      &
     &                  (Pair(i,j)-Pair(i-1,j))
#endif

The same changes are made to rhs_vbar.

This case is storm surge caused by a typhoon. As a result, the original code shows the water level dropping at the center of the typhoon. On the other hands, the modified code shows the water level rising at the center of the typhoon due to inverted barometer effect.
I hope this helps.

Original code:

Modified code:

[arango]

I don't have a shallow water (2D kernel) application to test this term. But, yes, a cyclone rotates counterclockwise around a center of low atmospheric pressure in the northern hemisphere. Consequently, the inverted barometer effect causes the water level to rise to form a dome around the low-pressure center. The opposite would happen over a cell of high air pressure causing the water level to drop.

Good catch, thank you I agree that we will need a minus sign in the barotropic pressure gradient above. However, we didn't discover this bug because we rarely run the shallow water setup, much less with atmospheric air pressure forcing.

Checking on my notes, I did notice the plus sign when I was coding the TIDE_GENERATING_FORCES adjacent term in step2d for a similar shallow water setup. I was supposed to set a test case with air pressure forcing, but it was low on my long priority list. So, I forgot. Thus, your application is a good one to test the ATM_PRESS term. Thanks.