[MITgcm-support] Problem with OBCs (bug ?)

[Nicolas Grisouard]

David,

You are right, I got a little bit confused in my boundary conditions and 
the first case was kind of useless...

I'll use a different ramp, but have one more question : what does the 
exf package improve when added to the obcs package ? I never used it, 
but it looks like it is designed to apply a forcing at the surface, 
which is not what I want. Am I wrong ?

Thanks,
Nicolas

David Hebert wrote:
> Nicolas,
>
> In your second and third case, are the U and V at the boundary 
> non-zero? If so, it looks like a vertical velocity W is begin 
> generated from continuity, since you can not expressly specify W. This 
> happens because you have zero velocity fields but velocity on the 
> boundary, causing a velocity gradient. I have run into this issue is 
> some of my work.
>
> If you want to start with your U and V fields (and W) at zero, a 
> solution is to use the exf package and start the boundary U and V at 
> zero. Then "ramp" the boundaries to the desired values. Another 
> possibility is to specify U and V fields the same as the boundary 
> values. This will result in zero velocity gradients, and not generate 
> W from continuity.
>
>
> David
>
> Nicolas Grisouard wrote:
>> Hi  everyone,
>>
>> I am running experiments about internal waves forced through open 
>> boundary condition on the western boundary (a little bit like the 
>> internal_wave experiment, but with more complex forcings).
>>
>> I have been getting a weird behaviour for the past months but only 
>> noticed it a few days ago, since this problem now induces unphysical 
>> features. I think I isolated one problem and think it involves the 
>> OBCS package.
>>
>> About the first issue. I tried to implement two forcings that are 
>> about the same. Configuration #1 (C1), in obcs_calc.F :
>>
>>        DO K=1,Nr
>>          DO J=1-Oly,sNy+Oly
>>             OBWu(J,K,bi,bj)=0. _d 0
>>     &            +Uinflow*ramptime2
>>     &            * sin(2*PI*(float(k)/float(Nr)-5.D-1))
>>     &            * sin(2.0 _d 0*PI*futureTime/obTimeScale)
>>             OBWv(J,K,bi,bj)=0. _d 0
>>     &            +Uinflow
>>     &            *f0/(2.0 _d 0*PI/obTimeScale)
>>     &            *cos(2. _d 0*PI*futureTime/obTimeScale )
>>             OBWt(J,K,bi,bj)=tRef(K)
>> #ifdef ALLOW_NONHYDROSTATIC     ! I do allow nonhydrosatic effects
>>             OBWw(J,K,bi,bj)=0.
>> #endif
>>          ENDDO
>>        ENDDO
>>
>> That is, in human-friendly format :
>> u(z,t) = U.f(t) . sin[2.pi.(z/H - 1/2)] . sin(2.pi.t/T)
>> everything_else = 0 as f0 = 0 (nonrotating case)
>> f(t) = function that smoothes the activation of the forcing and equal 
>> to 1.8% at t= 0
>> H = height of the domain
>> kz = 2*pi/H
>> T = period of the forcing so that frequency < N
>>
>> So we have a "stationary" forcing. At time t = 0, the fields are like 
>> they should be : the stratification is as set in gendata and the 
>> velocities fields are = 0.
>>
>> If we now take the exact same configuration except for the forcing 
>> that is now "propagative" (C2) (the one I need) :
>>
>>        DO K=1,Nr
>>          DO J=1-Oly,sNy+Oly
>>             OBWu(J,K,bi,bj)=0.D0
>>     &            +Uinflow*ramptime2
>>     &            * sin(2.D0 *PI*(float(k)/float(Nr)-5. _d -1
>>     &            - futureTime/obTimeScale))
>>             OBWv(J,K,bi,bj)=0. _d 0
>>     &            +Uinflow
>>     &            *f0/(2.0 _d 0*PI/obTimeScale)
>>     &            *cos(2. _d 0*PI*futureTime/obTimeScale )
>>             OBWt(J,K,bi,bj)=tRef(K)
>> #ifdef ALLOW_NONHYDROSTATIC
>>             OBWw(J,K,bi,bj)=0.
>> #endif
>>          ENDDO
>>        ENDDO
>>
>> Again :
>> u(z,t) = U.f(t). sin[2.pi.(z/H - 1/2 - t/T)]
>>
>> Now at time t = 0, the U and T fields are like they should be, but 
>> the W field shows someting odd, as displayed in figure 
>> propagative.jpg. The vertical velocity is non-zero at the column that 
>> is next to the point i = 1. This induces parasite emission of waves 
>> in the best scenario and maybe the blowing of the computation in the 
>> worst scenario.
>>
>>
>> Third configuration (C3), to complete C2 :
>> u(z,t) = U.f(t). sin[2.pi.(z/H - t/T)]
>>
>> I get a W field as shown in propagative2.jpg : same as in 
>> propagative.jpg but opposite values of velocities. I would have 
>> expected a node in the middle...
>>
>>
>>
>> I have no idea how to fix this and am ready to give any additional 
>> information to who could help me of course.

-- 
Nicolas GRISOUARD
PhD Student - ERES division
Laboratoire des Ecoulements Geophysiques et Industriels
BP 53
38041 Grenoble cedex 9 France
tel : +33 (0)476 825 037 - fax : +33 (0)476 825 271