Re: [ferret_users] potential vorticity from OGCM

[prescilla kurien <prescy@xxxxxxxxxxxxxxxxx>]

Hi William and Ryo,
I forgot to keep MKS units. My program works fine now.
Thank you very much your comments.
Prescilla
On Wed, 2009-07-29 at 06:57 -0700, William S. Kessler wrote:
> What Ryo said, exactly.
> 
> But if you still think you're not getting the right answer, look at  
> the pieces before assembling them:
> Is the relative vorticity reasonable? (shade it and overlay velocity  
> vectors to check)
> 
> One thing that worries me a little is the use of SET REGION,  
> including for K. Is that causing a problem when you are taking a  
> vertical derivative? Are the vertical grids the same for all quantities?
> If the other pieces look right, then try removing SET REGION/K=10,  
> defining the vertical derivative manually ((density[k=9]-density 
> [k=11])/2*<your interval>), and specifying zeta[k=10].
> 
> And you omitted the most important thing when using SET REGION in a  
> script:
> 
> The final line: CAN REG/ALL
> 
> It is far too easy to get an incorrect result when leftover REGIONS  
> are defined and you have forgotten that they remain. For example, if  
> K=10 remained specfied and you later did a vertical integral of  
> something else, you would get the wrong answer and not know it. I  
> never use SET REGION! It is straightforward to make the specification  
> within the calculation, or within the definition of variables. A  
> further reason to do this is that the JNL files are easier to  
> understand later if the region specifications are done in the  
> calculation, rather than somewhere far above it the listing.
> 
> One man's opinion.
> 
> Billy K
> 
> 
> On 29 Jul 09, at 1:04 AM, Ryo Furue wrote:
> 
> > Hi Prescilla,
> >
> > | The code I prepared is as below
> > |
> > | use 2014_hori_U_AM500_AH250.nc
> > | use 2014_hori_V_AM500_AH250.nc
> > | set region/X=78E:94E/Y=6N:23N/L=8/K=10
> > | use density.nc
> > | LET zeta=(VVEL[d=2,x=@ddc])-(UVEL[d=1,y=@ddc])  ! relative vorticity
> > | let deg2rad=(2*3.14)/360
> > | let omega = 7.292e-5
> > | let f=2*omega*sin(deg2rad*Y[g=UVEL[d=1]]) !planetary V
> > | let rho_zero=1026 ! referece density
> > | let grad=DENSITY[z=@ddc] ! density gradient
> > | fill (grad*(zeta+f))/rho_zero
> >
> > Looks good to me, if DENSITY, VVEL, and UVEL are in MKS units.
> > A potential (minor) problem is that if DENSITY, VVEL, and UVEL
> > are defined on different grids, then it's better to map them on
> > to a common grid before the final multiplication; for example
> >
> >   let pv = grad[g=zeta] * (zeta + f) / rho_zero
> >
> > Density should be a potential density but I don't think
> > that matters much for the upper ocean.
> >
> > By the way, I have a very minor piece of advice, which isn't
> > related to the cause of your problem.  I think it's better
> > not to throw away precision gratuitously: I would write
> > something along the lines of
> >
> >   let pi = 4*atan(1.0)  !! or acos(-1.0)
> >   let deg2rad = pi/180
> >   let omega = 2*pi/(24*60*60)
> >
> > Regards,
> > Ryo
>