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/ALLIt 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