Re: Geostrophic velocity

[Mick Spillane <Mick.Spillane@xxxxxxxx>]

Hello again Mohammed,
	I had a chance to look some more at your problem while waiting
for a long job to run.  The result is attached as a script geovelsect.jnl
which computes basically the same variables as does dynamic_height.jnl
but for an arbitrary "section" defined by a set of locations, i=1,2,...
The user provides the latitude and longitude values for each "station" on
the section in addition to identifying the temperature, salinity and
pressure variables.
	I've tested the results by extracting zonal and meridional
sections from the ocean_atlas_annual dataset and running them through
geovelsect.jnl and dynamic_height.jnl . The minor differences in the
results are I suspect due to my substitution of a great circle distance
calculation as a replacement for the @ddc (my 111.11 factor for km/deg
may differ from the equivalent in ferrets representation of global
distances). I haven't looked at the match of geovelz, or tested it with
datafiles containing multiple sections or representing the southern
hemisphere. I did test it with a rudimentary skewed section and it
worked fine.  If I made the change suggested in (3) below I could compare
the result with those of the EPIC oceanographic package we also use at
PMEL.
	A couple of points to note:
1) geovelsect.jnl provides only "geovel" (and geovelz) - the geostrophic
   velocity at right angles to the section - rather than the geo_u and
   geo_v components that dynamic_height.jnl can provide when a grid
   has both x and y extent.                   ^
2) to interpret the sign of geovel           / \   A positive
   refer to the cartoon --->                /_ _\  geovel value
   A positive value of geovel occurs   1-----|2|------3-----4--
   for flow from right to left through       |_|
   the section when looking toward
   increasing values of the station
   index. If the section is from west to east this corresponds to the
   usual positive northward convention.  But if the section is south to
   north, a positive geovel implies a westward flow, different from the
   oceanographic interpretation of U.  I didn't test for southern
   hemisphere sections where I suppose the opposite would apply.
3) To get values for easy comparison with dynamic_height.jnl I followed
   its use of centered differences in computing the geostrophic velocity.
   Inherent in this is an assumption that the grid points are evenly
   spaced! The usual oceanographic practice is to take each adjacent
   station pair and compute the geostrophic velocity mid-way between them.
   If I were planning to use geovelsect myself I would modify it to
   mimic a forward difference rather than the centered difference and
   add lonout,latout variables to hold the positions at which the geovel
   values apply.

Hope this stuff is of some use to you and others,
Regards,
Mick

|!!! Mick.Spillane@noaa.gov !!!|
|__Room 2070 Bldg#3 NOAA/PMEL__|
|____Phone_:_(206)526-6780_____|

! geostropic velocity calculation for section data. Mick Spillane 5-Aug-04     
cancel mode verify
!
! adapted from the original dynamic_height.jnl
!
! For a set of temperature-salinity data on an x-z grid, this routine
! computes the geostrophic velocity at right angles to the "section".
! The x-axis is just an index for the "stations" in the section, so
! it does not need to have any particular units. Position on the
! earth is given by latitude-longitude variables.
!
say " "
say "GEOVELSECT (Version 1.0):"
say " Computes geostrophic velocity at right angles to a section"
say " of temperature/salinity station profiles. A positive geovel"
say " (in the northern hemisphere)        ^"
say " is associated with flow from       / \   A positive"
say " right to left through the         /_ _\  geovel value"
say " section when looking toward  1-----|2|------3-----4--"
say " increasing station number.         |_|"
say " This version uses the equivalent of the centered difference"
say " used in dynamic_height.jnl. This assumes equal spacing of"
say " the stations and a better approach would probably be to"
say " compute geovel midway between each station pair."
say " " 
say "INSTRUCTIONS:"
say " The following variables should be defined by the user"
say "       my_temp : temperature on the x-z grid (degC)"
say "       my_salt : salinity on the x-z grid (psu or ppt)"
say "       my_lat  : latitudes for the section (degrees)"
say "       my_lon  : longitudes for the section (degrees)"
say "       my_pres : pressure in dbar (or depth in meters)"
say " "
say " then the following are computed"
say "       coriolis : the coriolis parameter at each station"
say "       deltax   : distance (in km) for centered difference"
say "       svanom   : the specific volume anomaly"
say "       dynht    : the dynamic height"
say "       dynhtz   :    (see the original dynamic_height.jnl)"
say "       geovel   : geostrophic velocity normal to section"
say "       geovelz  :    (depth dependent version)"
say " "

let SVanom = 1/rho_UN(my_salt,my_temp,my_pres) - 1/rho_UN(35,0,my_pres)
let/title="Dynamic Height(dyn-cm)" DYNHT  = 1E5 * SVanom[z=@din]
let/title="Dynamic Height(dyn-cm)" DYNHTz = 1E5 * (SVanom[z=@din]-SVanom[z=@iin])

let d2r=atan(1.)/45
let coriolis=d2r*SIN(d2r*my_lat)/120
let lat_1=d2r*my_lat[i=@shf:1] ; let lat_2=d2r*my_lat[i=@shf:-1]
let lon_1=d2r*my_lon[i=@shf:1] ; let lon_2=d2r*my_lon[i=@shf:-1]
let deltax=111.11*acos(sin(lat_1)*sin(lat_2)+cos(lat_1)*cos(lat_2)*cos(lon_2-lon_1))/d2r

let/title="Geostrophic Velocity (cm/s)" geovel=(dynHT[i=@shf:1]-dynHT[i=@shf:-1])/(100*coriolis*deltax)
let/title="Geostrophic Velocity (cm/s)" geovelz=(dynHTZ[i=@shf:1]-dynHTZ)/(100*coriolis*deltax)

set mode/last verify