geostrophic current

[Namba Takaya <takayanamba7@xxxxxxxxxxx>]

Dear Sir

I have tried to calculate the geostrophic current with the
script below by Dr. David Oxilia
submitted to the mailing list on 15th Jan 1999.
In addition, I have the World Ocean Atlas 1/4 (24 standard
levels) supported by Dr. Jaison Kurian.
However, I could not succeed in this calculation as some
axis errors come. Especially I do not have confidence of
two lines
set reg/i=660:720/j=120:200/k=1:24
and 
let rho = rho_un(ss[G=ggrid],tt[G=ggrid],z)
Could you give me some suggestions?

============= geostrophic current script ==========
! Set your data file (I assume you have a netCDF file of
model output)
!!!use SaltTemp
use woa01_salt.nc
use woa01_temp.nc
let ss=salt[d=woa01_salt.nc,l=1]  !!!WOA 1/4 salinity
let tt=temp[d=woa01_temp.nc,l=1]  !!!WOA 1/4 temperature
! Set parameters
let g = 9.81
let pi = 3.14159
let omega = 7.292e-5
let f = 2*omega*sin(y*pi/180)
define axis/x=160e:180e:0.25/unit=degree xxaxis
define axis/y=60s:40s:0.25/unit=degree yyaxis
define axis/z zaxis = {0, 10, 20, 30, 50, 75, 100, 125,
150,\
200, 250, 300, 400, 500, 600, 700, 800, 900, \
1000, 1100, 1200, 1300, 1400, 1500}
define grid/x=xxaxis/y=yyaxis/z=zaxis ggrid
! Set region to all space (use actual values for ni, nj
and nk)
!set reg/i=1:ni/j=1:nj/k=1:nk
set reg/i=660:720/j=120:200/k=1:24
!!!set reg/x=165e:180e/y=60s:40s/z=0:1500


! Compute density
! I'll assume you have two variables S and T in your data
file and that
! your vertical levels are in dbars (meters would be
pretty close)
let rho = rho_un(ss[G=ggrid],tt[G=ggrid],z)

! Compute vertical shears
let gamma = g/(rho*f)
let u_z = gamma*rho[y=@DDC]
let v_z = gamma*rho[x=@DDC]*(-1)

! Compute geostrophic field
! To compute absolute velocity you either have to know
some reference
! velocity or assume a level of no motion (as is
traditionally done)
! The level of no motion argument is what is used when
calculating
! the surface geostrophic field from "dynamic topography".
You can
! assume that the flow is zero at the bottom for example
and then
! integrate vertically to the level at which you want the
flow field.

! The flow field at the surface
let u = u_z[k=@DIN]
let v = v_z[k=@DIN]
! The flow field at the base of a 100 m mixed layer
assuming the
! depth of no motion is at 5000 m.
let um = u_z[z=5000:100@DIN]
let vm = v_z[z=5000:100@DIN]

Takaya Namba




__________________________________
STOP HIV/AIDS.
Yahoo! JAPAN Redribbon Campaign
http://pr.mail.yahoo.co.jp/redribbon/