How to force a 1D vertical case?

Discussion about coupled ecosystem models

Moderators: arango, robertson

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

How to force a 1D vertical case?

#1 Post by csherwood » Sat Nov 02, 2013 12:19 pm

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: 37
Joined: Fri Apr 02, 2004 4:46 pm
Location: USGS, Woods Hole, USA

Re: How to force a 1D vertical case?

#2 Post by csherwood » Fri Dec 27, 2013 3:37 pm

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
stress.png (7.44 KiB) Viewed 3221 times
Here is what they look like with splines on:
stress_splines.png
stress_splines.png (7.58 KiB) Viewed 3221 times
Chris Sherwood, USGS
1 508 457 2269

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

Re: How to force a 1D vertical case?

#3 Post by kate » Sat Dec 28, 2013 1:35 am

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

Post Reply