## How to force a 1D vertical case?

Moderators: arango, robertson

Message
Author
csherwood
Posts: 36
Joined: Fri Apr 02, 2004 4:46 pm
Location: USGS, Woods Hole, USA

### How to force a 1D vertical case?

I'd like to be able to simulate steady, uniform open-channel flow efficiently, and get the same results that 1DV models get for turbulence, stress, and velocity profiles when forced with a sloping sea surface. Small periodic domains (sed_toy, bio_toy) with body force or surface stresses don't produce the same vertical structure. Nudging is ok, but only if you know the right answer for the turbulence model, and does not allow experiments with the turbulence model to look at effects of, for example, wave-current interactions at the sea bed or dynamic tracers (e.g., sediment-induced stratification). Is there a way to do this correctly?
Chris Sherwood, USGS
1 508 457 2269

csherwood
Posts: 36
Joined: Fri Apr 02, 2004 4:46 pm
Location: USGS, Woods Hole, USA

### Re: How to force a 1D vertical case?

Answering my own question, with thanks to John Wilkin (who showed it to me) and John Warner (who found the key issue for bbl simulations: #undef splines).

You can use atmospheric pressure in ana_pair to impose a horizontal pressure gradient. Use it to make the dp/dx you want, recalling that in steady, uniform, unstratified, open-channel flow Taub = rho*g*h*ds/dx = h *dp/dx. You have to turn off the part that applies periodic boundary conditions for pair in ana_pair.h.

In sed_toy.h:

Code: Select all

``````# define atm_press
# define ana_pair
# undef splines /* splines messes with bbl profiles */
``````
In ana_pair.h:

Code: Select all

``````...
#elif defined SED_TOY
! 0.02 Pa/m corresponds to a slope of 1 m / 500 km
! in this case, dx and dy are 100 m, h = 10 m
dx = 100.0_r8
dpdx = 0.2_r8 !Pa (bottom stress will be h*dpdx)
DO j=JstrT-2,JendT+2 ! in older versions: JstrR-2,JendR+2
DO i=IstrT-2,IendT+2 ! IstrR-2,IendR+2 !
! convert Pa to millibars with 0.01
Pair(i,j)=1000.0_r8+0.01_r8*dpdx*FLOAT(i-1)*dx
END DO
END DO
#else
ana_pair.h: no values provided for Pair.
#endif

#if !defined SED_TOY
!  Exchange boundary data, except when SED_TOY is defined, per J. Wilkin
!
IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
CALL exchange_r2d_tile (ng, tile,                               &
&                          LBi, UBi, LBj, UBj,                     &
&                          Pair)
END IF
#endif
``````
Here is a plot of the eddy viscosity profile (modeled with GLS) and the stress profile, looking as they should.
stress.png (7.44 KiB) Viewed 3193 times
Here is what they look like with splines on:
stress_splines.png (7.58 KiB) Viewed 3193 times
Chris Sherwood, USGS
1 508 457 2269

kate
Posts: 3678
Joined: Wed Jul 02, 2003 5:29 pm
Location: IMS/UAF, USA

### Re: How to force a 1D vertical case?

Note that cpp defines are case sensitive and generally all upper case. Try ATM_PRESS and ANA_PAIR.