How to force a 1D vertical case?

[csherwood]

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?

[csherwood]

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 10830 times

Here is what they look like with splines on:

stress_splines.png (7.58 KiB) Viewed 10830 times

[kate]

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