-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathcalc_dir.ncl
88 lines (61 loc) · 2.16 KB
/
calc_dir.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
;;; NCL script to calculate wind direction from u & v components
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"
;; The following variables should be specified as command-line arguments
;; e.g.: ncl file=\"$file\" script.ncl
;ufile = "uas_CRCM_2001010103.nc"
;vfile = "vas_CRCM_2001010103.nc"
;outfile = "dir_CRCM_2001010103.nc"
finu = addfile(ufile, "r")
finv = addfile(vfile, "r")
system("rm -f "+outfile)
fout = addfile(outfile, "c")
filedimdef(fout,"time",-1,True) ;; make time dimension unlimited
; copy/set global attributes
att_names = getvaratts(finu)
do i = 0,dimsizes(att_names)-1
if(att_names(i) .eq. "history_of_appended_files") then
;; skip
else if(att_names(i) .eq. "title") then
fout@title = str_sub_str(finu@title, "10m Eastward Wind Velocity", "Near-Surface Wind Direction")
else
fout@$att_names(i)$ = fin@$att_names(i)$
end if
end if
end do
history = "Created " + systemfunc("date") + " by "+systemfunc("whoami")+"@"+systemfunc("hostname")+" using NCL script from source files "+ufile+" and "+vfile
fout@history = history
;; update unique tracking_id
fout@tracking_id = systemfunc("uuidgen")
; copy variables
var_names = getfilevarnames (finu) ;
do i = 0,dimsizes(var_names)-1
if (var_names(i) .ne. "uas") then
fout->$var_names(i)$ = finu->$var_names(i)$
end if
end do
; calculate wind direction
rad2deg = 45.0/atan(1.0)
dir = finu->uas
printVarSummary(dir)
dir = dir * 0.
printVarSummary(dir)
nt = dimsizes(dir&time)
printVarSummary(dir)
do t = 0, nt-1
tmp = atan2(finu->uas(t,:,:), finv->vas(t,:,:)) * rad2deg + 180
dir(t,:,:) = dir(t,:,:) + tmp
end do
varatts = (/"missing_value", "_FillValue", "coordinates", "grid_mapping"/)
do i = 0,dimsizes(varatts)-1
dir@$varatts(i)$ = finu->uas@$varatts(i)$
end do
dir@units = "degrees"
dir@units_convention = "degrees clockwise from north"
dir@long_name = "Near-Surface Wind Direction"
dir@standard_name = "wind_from_direction"
fout->dir = dir
exit
;; Copyright 2009-2012 Univ. Corp. for Atmos. Research
;; Author: Seth McGinnis, [email protected]