Hi Dennis. Thank you for the hint. My code is below.
Would the same go for calculating the du/dy term, or would it be
slightly different for the dlat and center_finite_diff loop step?
Sincerely,
Erik
;;;;;;;;;;;;;;;;;;;; NCL Script
; Script calculates dv/dx from Reanalysis II data
; This scripts has a lot of comments so that a user can learn as they go.
; 1) The script reads in a netcdf file that already contains variables
u and v in 4D
; 2) It gets the time, level, lat ,and lon dimensions of each variable.
; 3) It calculates the degree-increment space between each longitude point
; and converts this degree scale to radians.
; 4) It preallocates space to make a new array for dv/dx and makes meta data.
; 5) It calculates dv/dx while looping through each latitude.
; 6) It writes the new data to a new netcdf file.
; The script also contains many print and print summary statements
; so that the user can check data.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; **** User Input
file_in = "ncep2_uv"
file_out = "ncep2_dvdx_test"
; ___ Some Constants
pi = 4.0*atan(1.0)
radian = (pi/180)
print(radian)
radius = 6.37e6 ; redius of Earth, in meters
; **** End User input
; 1) reads in a netcdf file that already contains variables u and v in 4D
a = addfile(file_in+".nc","r")
u = a->U
v = a->V
printVarSummary(u)
printVarSummary(v)
; 2) Get variable dimentions from netcdf file
lat = a->lat
lon = a->lon
lev = a->lev
dt = dimsizes(u(:,0,0,0)) ; dimension of times
nlev = dimsizes(lev) ; dimension of levels
nlat = dimsizes(lat) ; dimention of latitudes
mlon = dimsizes(lon) ; dimention of longitudes
; 3) scale of longitude (space between each point)
dlon = (lon(2)-lon(1))
dlon = dlon*radian ; convert dlon to radians for calculation
print(dlon)
; 4) pre-allocate space for dv/dx
; ___ (i.e. Make new array the same size and dimensions as v)
dVoverdX = new ( dimsizes(v), double, v@_FillValue) ; float ncep2 ;
double merra
printVarSummary(dVoverdX)
; ---- META DATA (for writing out to file or plotting graphics)
dVoverdX!0 = "time"
dVoverdX!1 = "lev"
dVoverdX!2 = "lat"
dVoverdX!3 = "lon"
dVoverdX&lev = lev
dVoverdX&lat = lat
dVoverdX&lon = lon
dVoverdX@long_name = "zonal gradient of meridional wind"
dVoverdX@units = "m^2/s^2"
printVarSummary(dVoverdX)
; 5) Calculate dv/dx by looping over each latitude
do nl=0,nlat-1
dX = radius*cos(radian*lat(nl))*dlon ; constant at this latitude
dVoverdX(:,:,nl,:) = center_finite_diff_n (v(:,:,nl,:), dX , False,0,0)
end do
; result: dV/dX(time,level,lat,lon)
; ---- CHECK YOUR CALCULATIONS TO SEE IF YOU DID THEM CORRECTLY
printMinMax(dVoverdX,True)
printMinMax(dVoverdX(:,:,{5:15},{-30:10}),True) ; all levels
printMinMax(dVoverdX(:,{700},{5:15},{-30:10}),True) ; at 700 mb
; 6) SAVE AS NETCDF FILE
print("***************************************************************")
NCFILE = file_out+".nc"
print( "Writing out netcdf file: "+NCFILE)
system("/bin/rm -f "+NCFILE) ; remove any pre-existing file
ncdf = addfile(NCFILE ,"c") ; open output netCDF file
; make time an UNLIMITED dimension ; recommended for most applications
filedimdef(ncdf,"time",-1,True)
ncdf->dVoverdX=dVoverdX
end
On Sun, Oct 30, 2011 at 9:45 AM, Dennis Shea <shea@ucar.edu> wrote:
> center_finite_diff: See example 6
>
> On 10/28/11 9:19 PM, Erik N wrote:
>>
>> Hi.
>> I used uv2vr_cfd function for calculating relative vorticity.
>> http://www.ncl.ucar.edu/Document/Functions/Built-in/uv2vr_cfd.shtml
>>
>> I only want the dv/dx part values. Is there a simple way of getting
>> these values (this part of equation only) from the resulting values?
>> Thank you for the help.
>> Sincerely,
>> Erik
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Oct 31 01:07:21 2011
This archive was generated by hypermail 2.1.8 : Tue Nov 01 2011 - 13:43:04 MDT