Data Processing

1
Data Processing

• Dennis Shea
• National Center for Atmospheric Research

• NCAR is sponsored by the National Science
  Foundation

2

• DP Example multi-formatted data

• ISCCP HDF, binary 20yrs, 3hrly 80GB type
  byte
• calculations, regrid, monthly statistics gt
  netCDF, plot

• NCEP GRIB (same) 20yrs, 6hrly 25GB
• calculations, regrid, monthly statistics gt
  netCDF, plot

• CAM netCDF (same) 20yrs
• ensemble of model runs
• calculations, monthly statistics gt netCDF, plot

• Science datasets, calculations, graphics gt
  paper

3
Data Processing Outline

• Algebraic/logical expression operators
• Manual and automatic array creation
• if statements , do loops
• Built-in and Contributed functions
• User developed NCL functions/procedures
• User developed external procedures
• Sample processing
• Command Line Arguments CLAs
• Fortran external subroutines
• NCL as a scripting tool time permitting
• Global Variables time permitting

4
Algebraic Operators
Algebraic expression operators - Negation

Exponentiation Multiply
/ Divide Modulus integers
only Matrix Multiply Plus
- Minus gt
Greater than selection lt Less than
selection

• Use () to circumvent precedence rules
• All support scalar and array operations like
  f90
• is overloaded operator
• algebraic operator
• 5.3 7.95 ? 13.25
• concatenate string
• "alpha" (5.3 7) ? "alpha12.3

5
Logical Expressions

• Similar to f77

• Logical expressions formed
• by relational operators
• .le. (less-than-or-equal)
• .lt. (less-than)
• .ge. (greater-than-or-equal)
• .gt. (greater-than)
• .ne. (not-equal)
• .eq. (equal)
• .and. (and)
• .xor. (exclusive-or)
• .or. (or)
• .not. (not)

6
Manual Array Creation

• array constructor characters (//)
• a_integer (/1,2,3/)
• a_float (/1.0, 2.0, 3.0/) ,
  a_double (/1., 2, 3.2d /)
• a_string (/"abc",12345",hello, world"/)
  
• a_logical (/True, False,True/)
• a_2darray (/ (/1,2,3/), (/4,5,6/), (/7,8,9/)
  /)

• new function Fortran dimension, allocate C
  malloc
• x new (array_size/shape, type, _FillValue)
• _FillValue is optional assigned default if not
  user specified
• No_FillValue means no missing value assigned
• a new(3, float)
• b new(10, double, 1d20)
• c new( (/5, 6, 7/), integer)
• d new(dimsizes(U), string)
• e new(dimsizes(ndtooned(U)), logical)

• new and (//) can appear anywhere in script
• new is not used that often

7
Automatic Array Creation

• data importation via supported format
• u f-gtU
• same for subset of data u f-gtU(, 392, ,
  1020)
• meta data (coordinate array will reflect subset)

• variable to variable assignment
• y x y gt same size, type as x plus
  meta data
• no need to pre-allocate space for y

• functions
• return array no need to pre-allocate space
  
• grido f2fsh ( gridi, (/ 64,128/))
• gridi(10,30,73,144) ? grido(10,30,64,128)
• grido f2fsh_Wrap ( gridi, (/ 64,128/))
  contributed.ncl

8
Array Dimension Rank Reduction

• subtle point singleton dimensions eliminated
• let T(12,64,128)
• Tjan T(0, , ) ? Tjan(64,128)
  
• Tjan automatically becomes 2D Tjan(64,128)
• array rank reduced considered degenerate
  dimension
• all applicable meta data copied

• can override dimension rank reduction
• Tjan T(00,,) ? Tjan(1,64,128)
• TJAN new( (/1,64,128/), typeof(T),
  T_at__FillValue)
• TJAN(0,,) T(0,,)

• Dimension Reduction is a "feature" really ?

9
Array Syntax/Operators

• similar to f90/f95
• arrays must be same size and shape conform
• let A and B be (10,30,64,128)
• C AB
• D A-B
• E AB
• C, D, E automatically created if they did not
  previously exist
• let T and P be (10,30,64,128)
• theta T(1000/P)0.286 ? theta(10,30,64,128)
• let SST be (100,72,144) and SICE -1.8 (scalar)
• SST SST gt SICE f90 where (sst.lt.sice)
  sst sice
• the operation performed by lt and gt is
  (sometimes) called clipping
• use built-in functions whenever possible
• let T be (10,30,64,128) and P be (30) then
• theta T(1000/conform(T,P,1))0.286
• all array operations automatically ignore
  _FillValue

10
if blocks (1)

• if-then-end if (note end if has space)
• if ( all(a.gt.0.) ) then
• statements
• end if

• if-then-else-end if
• if ( any(ismissing(a)) ) then
• statements
• else
• statements
• end if

• lazy expression evaluation left-to-right
• if ( any(b.lt.0.) .and. all(a.gt.0.) ) then
• statements
• end if

11
if blocks (2)

• Technically, no else if block but .
• str "MAR
• if (str.eq."JAN") then
• print("January")
• else if (str.eq."FEB") then
• print("February")
• else if (str.eq."MAR") then
• print("March")
• else if (str.eq."APR") then
• print("April")
• else
• print("Enough of this!")
• end if must
  group all end if at end not very clear
• end if
• end if
• end if

12
loops

• do loop (traditional structure end do has
  space)
• do iscalar_start_exp, scalar_end_exp ,
  scalar_skip_exp
• do n nStrt, nLast ,stride
• statements
• end do 'end do' has a
  space
• if start gt end
• identifier 'n' is decremented by a positive
  stride
• stride must always be present when startgtend

• do while loop
• do while (x .gt. 100)
• statements
• end do

• break loop to abort
  f90 exit
• continue proceed to next iteration f90 cycle

• minimize loop usage in any interpreted language
• use array syntax, built-in functions, procedures
• use Fortran/C codes when efficient looping is
  required

13
Built-in Functions and Procedures(1 of 2)

• use whenever possible
• learn and use utility functions
• all, any, conform, ind, ind_resolve, dimsizes
• fspan, ispan, ndtooned, onedtond,
• mask, ismissing, where
• system, systemfunc use local system
• functions may require dimension reordering
• must use named dimensions to reorder

compute zonal and time average of variable
T(time,lev,lat,lon) ? T(0,1,2,3)
(zonal average requires rectilinear grid) no
meta data transfered
Tavg dim_avg_n( T, 0)
Tavg(lev,lat,lon)
14
Built-in Functions and Procedures(2 of 2)

• functions NO need to preallocate memory
• y wgt_runave (x, wgt, 0)
• if returning to pre-existing array must conform

• procedures MUST preallocate memory with new
• psi new ( dimsizes(u) , typeof(u) )chi new (
  dimsizes(u) , typeof(u) )uv2sfvpg(u,v,psi,chi)

• functions may be imbedded, procedures can not
• keep code simple avoid 'deep' imbedding

example of a deep imbed x f2gsh( fo2fsh(
fbinrecread(f,6,(/9,18,72,144/),
float)),(/nlat,mlon/),42) without deep
imbedding G fbinrecread(f, 6,
(/1,18,72,144/), "float") g42 f2gsh( fo2fsh(G),
(/nlat,mlon/),42) delete (G)
15
dimsizes(x)

• returns the dimension sizes of a variable
• will return 1D array of integers if the array
  queried is multi-dimensional.

Variable dimt Type integer Total Size 16
bytes 4 values Number of dimensions 1 Dimensions
and sizes(4) (0) 12 (1) 25 (2) 116 (3) 100 (0)
rank4
fin addfile(in.nc,r) t
fin-gtT dimt dimsizes(t) print(dimt)
rank dimsizes(dimt) print ("rank"rank)
16
ispan( startinteger, finishinteger,
strideinteger )

• returns a 1D array of integers
• beginning with start and ending with finish.

time ispan(1990,2001,2) print(time)
Variable time Type integer Number of
Dimensions 1 Dimensions and sizes(6) (0)
1990 (1) 1992 (2) 1994 (3) 1996 (4) 1998 (5) 2000
17
fspan( start, finish, npts )
Variable b Type float Number of Dimensions
1 Dimensions and sizes(10) (0) 1.0 (1) 1.444 (2)
1.888 (3) 2.333 (4) 2.777 (5) 3.222 (6) 3.666 (7)
4.111 (8) 4.555 (9) 5.0

• returns a 1D array of
• evenly spaced float/double
• values
• npts is the integer number
• of points including start
• and finish values

b fspan(1, 5, 10) print(b)
18
ismissing

• must be used to check for _FillValue attribute
• if (x.eq.x_at__FillValue) will not work

• x (/ 1,2, -99, 4, -99, -99, 7/)
  x_at__FillValue -99
• xmsg ismissing(x) print(xmsg)
• (/False, False, True, False, True, True,
  False /)

• often used in combination with array functions
• if (any( ismissing(x) )) then else end if
• nFill num( ismissing(x) )
• nVal num( .not. ismissing(x) )

• if (any(ismissing(xOrig) )) then
• xNew linint2_Wrap(xin,yin,xOrig,True,xout,yo
  ut,0)
• else
• xNew f2fsh_Wrap (x, (/128,256/) )
• end if

19
mask

• sets values to _FillValue that DO NOT equal mask
  array

load "NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_cod
e.ncl" load "NCARG_ROOT/lib/ncarg/nclscripts/cs
m/gsn_csm.ncl" in addfile(atmos.nc","
r") ts in-gtTS(0,,) oro
in-gtORO(0,,) mask ocean ocean0, land1,
sea_ice2 land_only ts land_only
mask(ts,oro,1)

• NCL has 1 degree land-sea mask available
  landsea_mask
• load "NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_ut
  il.ncl
• flags for ocean, land, lake, small island,
  ice shelf

20
where

• performs array assignments based upon a
  conditional array
• function where(conditional_expression \
• , true_value(s)
  \
• , false_value(s) )
• similar to f90 where statement
• components evaluated separately via array
  operations

q is an array qlt0 gt qq256 f90
where(q.lt.0) qq256 q where (q.lt.0,
q256, q)
x where (T.ge.0 .and. ismissing(Z) , a25 ,
1.8b)
salinity where (sst.lt.5 .and. ice.gt.icemax
\ , salinity0.9,
salinity)
can not do y where(y.eq.0, y_at__FillValue,
1./y) instead use y 1.0/where(y.eq.0,
y_at__FillValue, y)
21
dim__n dim_

• perform common operations on an array dimension
• avg, var, sum, sort, median, rmsd, ..
• dim__n functions are newer (added to v5.1.1)
• operate on a user specified dimension
• use less memory
• dim_ functions are original interfaces
• operate on rightmost dimension only
• may require dimension reordering

• consider x(time,lat,lon) gt x(0,1,2)
• function dim_avg_n( x, n )
• xZon dim_avg_n( x, 2 ) gt xZon(ntim,nlat)
• xTim dim_avg_n( x, 0 ) gt xTim(nlat,mlon)
• function dim_avg ( x )
• xZon dim_avg( x ) gt xZon(ntim,nlat)
• xTim dim_avg( x(lat,lon,time) ) gt
  xTim(nlat,mlon)

22
conform, conform_dims
Array operations require that arrays conform

• function conform( x, r, ndim )
• function conform_dims( dims, r, ndim )
• expand array (r) to match (x) or dimensions
  sizes (dims)
• ndim scalar or array indicating which
  dimension(s) of x or dims match the dimensions of
  r
• array r is broadcast to another array size

x(ntim,klev,nlat,mlon), w(nlat) wx conform (x,
w, 2) wx(ntim,klev,nlat,mlon) q xwx
q x conform (x, w, 2) wx
conform_dims ( (/ntim,klev,nlat,mlon/) , w, 2) wx
conform_dims ( dimsizes(x), w, 2)
T(ntim, klev,nlat,mlon,), dp(klev) dpT conform
(T, dp, 1) dpT(ntim,klev,nlat,mlon) T_wg
tAve dim_sum_n (TdpT, 1)/dim_sum(dpT, 1)
23
ind

• ind operates on 1D array only
• returns indices of elements that evaluate to
  True
• generically similar to IDL where and Matlab
  find returns indices

let x, y, z z_at__FillValue create
triplet with only good values iGood
ind (.not. ismissing(z) ) xGood
x(iGood) yGood y(iGood) zGood
z(iGood)
let a, return subscripts can be on lhs
ii ind (a.gt.500 ) a(ii) 3a(ii) 2

• Should check the returned subscript to see if it
  is missing
• if (any(ismissing(ii))) then . end if

24
ind, ndtooned, onedtond

• ind operates on 1D array only
• returns indices of elements that evaluate to
  True
• if nD use with ndtooned reconstruct with
  onedtond, dimsizes
• generically similar to IDL where and Matlab
  find returns indices

let q and x be nD arrays q1D ndtooned
(q) x1D ndtooned (x) ii
ind(q1D.gt.0. .and. q1D.lt.5) jj
ind(q1D.gt.25) kk ind(q1D.lt. -50)
x1D(ii) sqrt( q1D(ii) ) x1D(jj) 72
x1D(kk) -x1D(kk)3.14159 x
onedtond(x1D, dimsizes(x))
25
ind_resolve

• ind_resolve will return the nD indices
  corresponding to ind

• print indices p gt 50 let p(time,lat,lon)
• if ( any(p.gt.50.) ) then
• p1D ndtooned (p)
• i50 ind(p1D.gt.50.)
• ir ind_resolve( i50, dimsizes(p))
  2D npts,3
• print (ind where pgt 50 ir(,0)
  ir(,1) ir(,2))
• print time,lat,lon p gt 50
• print (time(ir(,0)) lat(ir(,1))
  lon(ir(,2)) )
• end if

• CAVEAT should not be use for vector subscripting

u(ir(,0),ir(,1),ir(,2)) 5 wont get what
you might expect
26
str_ string functions

• many new str_ functions
• http//www.ncl.ucar.edu/Document/Functions/string
  .shtml
• greatly enhance ability to handle strings
• can be used to unpack complicated string arrays

• x (/ u_052134_C, q_1234_C,
  temp_72.55_C/)
• var_x str_get_field( x, 1, _)
• result var_x (/u, q, temp/)
  strings
• --------
• col_x str_get_cols( x, 2, 4)
• result col_x (/052, 123, mp_ /)
  strings
• ---------
• N toint( str_get_cols( x(0), 3, 7) )
  N52134 (integer)
• T tofloat( str_get_cols( x(2), 5,9 )
  T72.55 (float)

27
date cd_calendar, cd_inv_calendar

• Date/time functions
• http//www.ncl.ucar.edu/Document/Functions/date.sh
  tml
• cd_calendar, cd_inv_calendar
• deprecated ut_calendar, ut_inv_calendar

• time (/ 17522904, 17522928, 17522952/)
• time_at_units hours since 1-1-1 00000.0
• date cd_calendar(time,0)
• print(date)
• Variable date
• Type float
• Total Size 72 bytes 18 values
• Number of Dimensions 2
• Dimensions and sizes 3 x 6
• (0,05) 2000 1 1 0 0 0
• (1,05) 2000 1 2 0 0 0
• (2,05) 2000 1 3 0 0 0

• date cd_calendar(time,-2)
• print(date)
• Variable date
• Type integer
• Total Size 12 bytes 3 values
• Number of Dimensions 1
• Dimensions and sizes 3
• (0) 20000101
• (1) 20000102
• (2) 20000103

• TIME cd_inv_calendar (iyr, imo, iday, ihr,
  imin, sec \
• ,hours
  since 1-1-1 00000.0 ,0)

28
system, systemfunc (1 of 2)

• system passes to the shell a command to perform
  an action
• NCL executes the Bourne shell (can be changed)

• create a directory if it does not exist
  (Bourne shell syntax)
• DIR /ptmp/shea/SAMPLE
• system (if ! test d DIR
  then mkdir DIR fi)
• same but force the C-shell (csh) to be used
• the single quotes () prevent the Bourne shell
  from interpreting csh syntax
• system ( csh c if (! d DIR)
  then mkdir DIR endif )
• execute some local command
• system (msrcp n mss/SHEA/REANALYSIS/
  /ptmp/shea)
• system (convert foo.eps foo.png
  /bin/rm foo.eps )
• system (ncrcat -v T,Q foo.nc FOO.nc
  )
• system (/bin/rm f file_name)

29
system, systemfunc (1 of 2)

• systemfunc returns to NCL information from the
  system
• NCL executes the Bourne shell (can be changed)

• UTC systemfunc(date)
• Date systemfunc(date a mdy HM )
  single quote
• fils systemfunc (cd /some/directory ls
  foonc) multiple cmds
• city systemfunc ("cut -c100-108 " fname)

30
User-built Functions and Procedures(1 of 4)

• two ways to load existing files w functions/proc
• load "/path/my_script.ncl"
• use environment variable NCL_DEFAULT_SCRIPTS_DIR
  
• must be loaded prior to use
• unlike in compiled language
• avoid function conflict (undef)

undef ("mult") function mult(x1,x2,x3,x4) begin
return ( x1x2x3x4) end
undef ("mult") function mult(x1,x2,x3,x4) begin
return ( x1x2x3x4) end begin x
mult(4.7, 34, 567, 2) end
load /some/path/mult.ncl begin x mult(4.7,
34, 567, 2) end
31
User-Built Functions and Procedures(2 of 4)

• Development process similar to Fortran/C/IDL
• General Structure

undef ("function_name") optional
function function_name (declaration_list) local
local_identifier_list optional
begin
required statements return (return_value)
required end
required
undef ("procedure_name") optional
procedure procedure_name (declaration_list) local
local_identifier_list optional
begin
required statements
end
required
32
User-Built Functions and Procedures (3 of 4)

• arguments are passed by reference fortran
• constrained argument specification
• require specific type, dimensions, and size
• procedure ex(datafloat,reslogical,textstr
  ing)
• generic specification
• type only
• function xy_interp(x1numeric, x2numeric)
• no type, no dimension specification
• procedure whatever (a, b, c)
• combination
• function ex (dfloat, xnumeric, wksgraphic,
  y2, a)
• function prototyping
• built-in functions are prototyped

33
User-Built Functions and Procedures (4 of 4)

• additional (optional) arguments possible
• attributes associated with one or more arguments
• often implemented as a separate argument (not
  required)
• procedure ex(datafloat, textstring,optArg
  logical)

procedure ex(data, text, optlogical) begin
if (opt .and. isatt(opt,scale)) then
d dataopt_at_scale end if if (opt
.and. isatt(opt,wgts)) then end
if if (opt .and. isatt(opt,array)) then
xloc3D opt_at_array_3D nD arrays end if
must be local before use end
optArg True optArg_at_scale
0.01 optArg_at_add 1000 optArg_at_wgts
(/1,2,1/) optArg_at_name sample optArg_at_array
array_3D ex(x2D, Example, optArg)
34
Computations and Meta Data

• computations can cause loss of meta data
• y x variable to variable
  transfer all meta copied
• T T273 T retains all meta data
• T_at_units "K" user responsibility to
  update meta
• z 5x z will have no meta data

• built-in functions cause loss of meta data
• Tavg dim_avg_n(T, 0)
• s sqrt(u2 v2)

• vinth2p is the exception
• retains coordinate variables
• http//www.cgd.ucar.edu/csm/support/Data_P/vert_i
  nterp.shtml
• hybrid to pressure (sigma to pressure) other
  examples

35
Ways to Retain Meta Data(1 of 3)

• use copy functions in contributed.ncl
• copy_VarMeta (coords attributes)
• copy_VarCoords
• copy_VarAtts

load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" begin f addfile("dummy.nc", "r")
x f-gtX
(ntim,nlat,mlon) ----------------
calculations----------------------------
xZon dim_avg _n(x, 2)
xZon(ntim,nlat) ----------------copy meta
data-------------------------- copy_VarMeta(x,
xZon) xZon(time,lat) end
36
Ways to Retain Meta Data(2 of 3)

• f2fosh_Wrap
• g2gshv_Wrap
• g2fshv_Wrap
• f2gshv_Wrap
• f2fshv_Wrap
• f2foshv_Wrap
• linint1_Wrap
• linint2_Wrap
• linint2_points_Wrap
• eof_cov_Wrap
• eof_cov_ts_Wrap
• zonal_mpsi_Wrap
• etc

• use wrapper functions (eg)
• dim_avg_n_Wrap
• dim_variance_n_Wrap
• dim_stddev_n_Wrap
• dim_sum_n_Wrap
• dim_rmsd_n_Wrap
• smth9_Wrap
• g2gsh_Wrap
• g2fsh_Wrap
• f2gsh_Wrap
• f2fsh_Wrap
• natgrid_Wrap

load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl f addfile("dummy.nc", "r") x
f-gtX
time,lev,lat,lon (0,1,2,3) xZon
dim_avg_n_Wrap(x, 3) xZon will have meta
data
37
Ways to Retain Meta Data(3 of 3)

• use variable to variable transfer dimension
  reduction to prefine array before calculation
• requires that user know a priori the array
  structure

load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" f addfile("dummy.nc", "r") x
f-gtX
x(time,lev,lat,lon) -------- var-to-var
transfer dim reduction-------- xZon
x(,,,0)
xZon(time,lev,lat) ---------------------calculat
ions------------------------- xZon
dim_avg_n (x, 0) xZon_at_op "Zonal Avg
"x_at_long_name

• xZon will have all appropriate meta data of x
• NCL will add an attribute here xZon_at_lon
  lon(0)

38
Example read fortran binary, compute psi/chi,
output binary
u fbinrecread ("UV.bin",0,
(/120,18,64,128/),"float") v
fbinrecread ("UV.bin",1, (/120,18,64,128/),"float"
) -----------------------------------------------
---------- calculate psi and chi as inverse
lapacian of divergence --------------------------
------------------------------- psi ilapsG
(uv2vrG(u,v), 0) chi ilapsG (uv2dvG(u,v),
0) ---------------------------------------------
------------ output binary data
------------------------------------------------
---------- fbinrecwrite ("psichi.bin", -1,
psi ) -1 means append fbinrecwrite
("psichi.bin", -1, chi )
39
Ex compute PSI/CHI add meta data
load NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl begin optional f
addfile (UV300.nc, r) u f-gtU v
f-gtV calculate psi and chi

psi ilapsG (uv2vrG(u,v), 0) chi
ilapsG (uv2dvG(u,v), 0) copy coordinate
variables using function in contributed.ncl
copy_VarCoords(u, psi) copy_VarCoords(u,
chi) create unique attributes psi_at_long_name
"PSI" chi_at_long_name "CHI" psi_at_units
"m2/s" chi_at_units
m2/s scale values for plotting scale
1.e06 incorporate this
into units label psi psi/scale chi
chi/scale .. plot .. end only if
begin is present
40
regrid Spherical Harmonics

• some regrid functions use spherical harmonics
• must have global grid to use these functions
• regular functions strip meta data
• wrapper versions preserve attributes
• and create coordinate variables
• no missing data allowed

load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" f addfile
("/fs/cgd/data0/shea/CLASS/T2m.nc", "r") T63
f-gtT T42 g2gsh_Wrap(T63, (/64,128/), 42)
T25 g2fgsh_Wrap(T63, (/73,144/) )
41
regrid bilinear interpolationlinint2

• Cartesian, global or limited area grids
• Must be grids that can be representd by one-dim
  coords
• wrapper versions preserve attributes
• and create coordinate variables
• missing data allowed

load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" f addfile
("/fs/cgd/data0/shea/CLASS/T2m.nc", "r") T
f-gtT TLI linint2_Wrap(Tlon, Tlat,
T, True, LON, LAT, 0 )
42
regrid areal conservative interpolationarea_cons
erve_remap, area_conserve_remap_Wrap

• global rectilinear grids only
• _Wrap preserves attributes creates coordinate
  variables
• missing data (_FillValue) NOT allowed

In particular, use for (say) flux or
precipitation interpolation
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" f addfile (GPCP.nc", "r") p
f-gtPRC P area_conserve_remap_Wrap (plon,
plat, p \
,newlon, newlat,
False)
43
regrid areal average interpolationarea_hi2lores,
area_hi2lores_Wrap

• rectilinear grids can be limited area
• _Wrap preserves attributes creates coordinate
  variables
• missing data allowed
• designed for TRMM data

NOT strictly conservative but close 50S to
50N
Use area_hi2lores for highly structured fields gt
lower res Use linint2 for low resolutiongt high
resolution
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" f addfile (trmm.nc", "r") p
f-gtPRC P area_hi2lores_Wrap (plon, plat, p,
True, LON, LAT, 0 )
44
NCL-ESMF regridding

• Integrated in conjunction with NOAA Cooperative
  Institute for Research in Environmental Sciences
• Available in NCL V6.1.0 (March/April 2012)
• Works with a wide variety of structured and
  unstructured grids
• Multiple interpolation methods available
• Bilinear
• Conservative
• Patch
• Can handle masked points
• Better treatment for values at poles
• Works on global or regional grids
• Can run in parallel or single-threaded mode

45
ESMF versus other regridding software

• Is ESMF regridding better than SCRIP regridding?
• ESMF takes extra care at the pole does
  interpolation in 3D cartesian space.
• ESMF's 'patch' method has no equivalent in SCRIP.
  Weight generation is slow but the results can be
  used to estimate derivatives (i.e. gradients,
  ocean curl)
• The conserve method is fast even on one
  processor
• NCL can be used to regrid data on GRIB1/2 and
  HDF4/5 files, and written to netCDF

46
(No Transcript)
47
(No Transcript)
48
(No Transcript)
49
What is ESMF?

• http//www.earthsystemmodeling.org/

• The Earth System Modeling Framework (ESMF)
  collaboration is building high-performance,
  flexible software infrastructure to increase ease
  of use, performance portability,
  interoperability, and reuse in climate, numerical
  weather prediction, data assimilation, and other
  Earth science applications.The ESMF defines an
  architecture for composing complex, coupled
  modeling systems and includes data structures and
  utilities for developing individual models.

• They have developed some fast parallel regridding
  software that works on structured and
  unstructured grids.

50
Interpolation methods

• "bilinear" - the algorithm used by this
  application to generate the bilinear weights is
  the standard one found in many textbooks. Each
  destination point is mapped to a location in the
  source mesh, the position of the destination
  point relative to the source points surrounding
  it is used to calculate the interpolation
  weights.
• "patch" - this method is the ESMF version of a
  technique called "patch recovery" commonly used
  in finite element modeling. It typically results
  in better approximations to values and
  derivatives when compared to bilinear
  interpolation.
• "conserve" - this method will typically have a
  larger interpolation error than the previous two
  methods, but will do a much better job of
  preserving the value of the integral of data
  between the source and destination grid.

51
poisson_grid_fill

• replaces all _FillValue grid ponts
• Poissons equation solved via relaxation
• values at non-missing locations are used as
  boundary cond.
• works on any grid with spatial dimensions

in addfile (Ocean.nc","r") sst
in-gtSST poisson_grid_fill (sst, True, 1,
1500, 0.02, 0.6, 0)
52
Example Arbitrary Transect
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" begin
open file and read in data diri
/fs/scd/home1/shea/ncldata_input/ fili
01-50.nc f addfile(dirifili , "r")
T f-gtT
T(time,lev,lat,lon)
create arrays of lat and lon
points lonx (/ 295, 291.05, 287.4 , 284,
281, 274.5, 268, 265 /) laty (/ -30,
-25.2, -20.3, -15, -10.3, 0.0, 9.9, 15
/)
must have a regular grid
interpolate data to given
lat/lon points transect linint2_points_Wrap
(Tlon, Tlat, T, True,lonx,laty, 0)
transect(time,lev,
8) end ncl lt PR_ex04.ncl
53
Example Compositing
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl t1 (/ 15, 37, 95, 88,90 /)
indices of desired dates t2 (/ 1,
22, 31, 97, 100, 120/) f
addfile(01-50.nc, "r") T1 f-gtT(t1,,,)
T(time,lev,lat,lon)
T2 f-gtT(t2,,,)
composite
averages T1avg dim_avg_n_Wrap(T1, 0)
(lev,lat,lon) T2avg dim_avg_n_Wrap(T2, 0)
Tdiff T2avg
trick to transfer meta data Tdiff T2avg -
T1avg difference
Tdiff_at_long_name T2_at_long_name composite
difference ------ Also use coordinate
subscripting let time have units yyyymm t1
(/ 190401, 191301, 192001, , 200301/) T1
f-gtT(t1,,,))
54
Ex Climatology, Seasonal Average, EOFs
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
uted.ncl" begin f addfile
(data_in.nc", "r") P f-gtprc
ntime dimsizes(Ptime) calculate monthly
climatology and monthly standard deviation
Pclm clmMonTLL(P) Pstd stdMonTLL
(P) reorder array to calculate running average
Psea runave (P(lat,lon,time), 3, 1
) Psea_at_time_op "3-month running mean
seasonal ave area weight wlat
sqrt(cos(0.01745329Psealat) ) wPsea Psea

transfer meta data wPsea wPseaconform(wPsea,
wlat, 0) wgt at all grid pts calculate
EOFs eof eofunc (wPsea(,,0ntime-112
), 3, opt_eof) eofTs eofunc_ts
(Psea(,,0ntime-112) , eof, opt_ts) end

• eofunc_varimax eofunc_varimax_reorder

55
Empirical Orthogonal Functions (EOFs)

• principal components, eigenvector analysis

• provide efficient representation of variance
• May/may not have dynamical information

• successive eigenvalues should be distinct
• if not, the eigenvalues and associated patterns
  are noise
• 1 from 2, 2 from 1 and 3, 3 from 2 and 4, etc
• North et. al (MWR, July 1982 eq 24-26) provide
  formula

• geophysical variables spatial/temporal
  correlated
• no need sample every grid point
• no extra information gained
• oversampling increases size of covar matrix
  compute time

• patterns are domain dependent

56

• Calculating EOFS, writing a NetCDF file (next
  page)

• load "NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_cod
  e.ncl"
• load "NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm
  .ncl"
• load "NCARG_ROOT/lib/ncarg/nclscripts/csm/contrib
  uted.ncl"
• f addfile("erai_1989-2009.mon.msl_psl.nc","r")
  open file
• p f-gtSLP(12,090,)
  (20,61,240)
• w sqrt(cos(0.01745329platitude) )
  weights(61)
• wp pconform(p, w, 1)
  wp(20,61,240)
• copy_VarCoords(p, wp)
• x wp(latitude,longitude,time)
  reorder data
• neof 3
• eof eofunc_Wrap(x, neof, False)
• eof_ts eofunc_ts_Wrap (x, eof, False)
• printVarSummary( eof ) examine
  EOF variables
• printVarSummary( eof_ts )

57

• Variable eof
• Type float
• Total Size 175680 bytes
• 43920 values
• Number of Dimensions 3
• Dimensions and sizes evn 3 x latitude 61
  x longitude 240
• Coordinates
• evn 1..3
• latitude 0..90
• longitude 0..358.5
• Number Of Attributes 6
• eval_transpose ( 47.2223, 32.42917, 21.44406
  )
• eval ( 34519.5, 23705.72, 15675.61 )
• pcvar ( 26.83549, 18.42885, 12.18624 )
• matrix covariance
• method transpose
• _FillValue 1e20

• printVarSummary output

58

• writing a NetCDF file

• Create netCDF no define mode simple approach
  
• system("/bin/rm -f EOF.nc") rm any
  pre-existing file
• fout addfile("EOF.nc", "c") new
  netCDF file
• fout_at_title "EOFs of SLP 1989-2009"
• fout-gtEOF eof
• fout-gtEOF_TS eof_ts

• Graphics http//www.ncl.ucar.edu/Applications/Scr
  ipts/eof_2.ncl

59
(No Transcript)
60
EOFs via NCL

• eofunc when variances are important
• eofunc (optTrue, opt_at_jopt1) when patterns
  important
• eofunc_ts creates time series of amplitudes
• eofunc_varimax
• - rotate EOFs using Kaiser varimax criterion
• Use _Wrap versions in contributed.ncl

Example http//www.ncl.ucar.edu/Applications/eof
.shtml

• let x(time,lat,lon) and X x(lat , lon ,
  time )
• wgt grid points gt X Xconform(X, wgt_lat, 0)
• eof eofunc_Wrap (X, neof, opteof)
• eofJan eofunc_Wrap (X(, , 0ntim-112),
  neofeof,opt) 12th
• eofJul eofunc_Wrap (X(, , 6ntim-112),
  neof, opt)
• e_ts eofunc_ts_Wrap(x(lat , lon , time
  ) , eof, optts)

• CDO cdo eof,4 anom_file eval_file
  eof_file

61
SVDs via NCL

• svd_lapack general rectangular matrix
• svdcov (svdstd) Bretherton-Smith-Wallace
• svdcov_sv (svdstd_sv) return left/right sing
  vectors

62
Post-processing Tools NCL (WRF-ARW Only)

• WRF Weather Research and Forecast Model

• Cindy Bruyère wrfhelp_at_ucar.edu

63
WRF Functions wrf_

• Special WRF NCL Built-in FunctionsMainly
  functions to calculate diagnosticsSeldom need to
  use these directly slp wrf_slp( z, tk, P,
  QVAPOR )
• Special WRF functionsDeveloped to make it easier
  to generate plotsNCARG_ROOT/lib/ncarg/nclscripts
  /wrf/WRFUserARW.ncl slp wrf_user_getvar(nc_fil
  e,slp,time)

64
WRF Generate Plots A good start - OnLine Tutorial

• http//www.mmm.ucar.edu/wrf/
• OnLineTutorial/
• Graphics/
• NCL/index.html

65
Command Line Arguments CLAs

• CLAs are NCL statements on the command line
• http//www.ncl.ucar.edu/Document/Manuals/Ref_Manua
  l/NclCLO.shtml

ncl tStrt1930 lev(/250, 750/) varT
fNamfoo.nc sample.ncl
if (.not. isvar(fNam") .and. (.not.
isvar(var") ) then print(fNam and/or
variable not specified exit) exit
end if f addfile (fNam, r)
read file x f-gtvar
read variable if (.not.
isvar("tStrt")) then CLA? tStrt
1900 default end
if if (.not. isvar("lev")) then
CLA? lev 500
default end if
66
External Fortran-C Codes
external codes, Fortran, C, or local commerical
libraries may be executed from within NCL

• generic process
• develop a wrapper (interface) to transmit
  arguments
• compile the external code (eg, f77, f90, cc)
• link the external code to create a shared object
• process simplified by operator called WRAPIT

• specifying where shared objects are located
• external statement
• external /dir/code.so
• system environment variable
• LD_LIBRARY_PATH
• NCL environment variable
• NCL_DEF_LIB_DIR
• external functions need not have before use

67
Create/Use a Fortran Shared Object (1 of 2)

• Step 1
• bracket f77 subroutine argument declarations
  with interface deliminators
• only codes actually called from NCL need special
  interface deliminators
• no Fortran argument can be named data (bug)

C NCLFORTSTART subroutine foo (
xin,xout, mlon, nlat, text) integer
mlon, nlat real xin(mlon,nlat),
xout(mlon,nlat) character() text C
NCLEND rest of fortran code may include
many subroutines or other declarations
68
Create/Use a Fortran Shared Object (2 of 2)

• Step 2 create shared object using WRAPIT utility
• WRAPIT foo.f
• WRAPIT specified_f90_compiler foo.stub foo.f90
• WRAPIT pg foo.stub foo.f90

• Step 3 add external statement to NCL script
• external SO_NAME "path_name"
• SO_NAME is arbitrary (capital by
  convention)
• external DEMO "./foo.so" (".so" by
  convention)

• Step 4 invoking the shared object in the script
• SO_NAMEfortran_name(arguments)
• DEMOfoo(x,y,mlon,nlat,label)

69
what WRAPIT does

• automatically creates NCL-fortran interface(s)

• uses wrapit77 to create C interface f77 syntax
• wrapit77 lt foo.f gt! foo_W.c

• only uses code enclosed between delimeters
• input can be code fragment(s) or full
  subroutine(s)

• compiles and creates object modules
• nhlcc -c foo_W.c ? foo_W.o
• nhlf90 -c foo.f ? foo.o

• creates dynamic shared object .so using ld
• SGI ld -64 -o foo.so -shared foo_W.o foo.o
  -fortran
• SUN /usr/ccs/bin/ld -o foo.so foo_W.o
  foo.o -G -lf90 -L /opt/SUNWspro/lib -l sunperf

• removes extraneous intermediate files

• WRAPIT h ltreturngt will show options and
  examples

70
NCL/Fortran Argument Passing

• arrays NO reordering required
• x(time,lev,lat,lon) ltmapgt x(lon,lat,lev,time)

• ncl x(N,M) gt value lt x(M,N) fortran
  M3, N2
• x(0,0) gt 7.23 lt x(1,1)
• x(0,1) gt -12.5 lt x(2,1)
• x(0,2) gt 0.3 lt x(3,1)
• x(1,0) gt 323.1 lt x(1,2)
• x(1,1) gt -234.6 lt x(2,2)
• x(1,2) gt 200.1 lt x(3,2)

• numeric types must match
• integer ltgt integer
• double ltgt double
• float ltgt real

• Character-Strings a nuisance C,Fortran

71
Example Linking to Fortran 77
STEP 1 quad.f C NCLFORTSTART
subroutine cquad(a,b,c,nq,x,quad)
dimension x(nq), quad(nq) C NCLEND do
i1,nq quad(i) ax(i)2 bx(i)
c end do return end C
NCLFORTSTART subroutine prntq (x, q, nq)
integer nq real x(nq), q(nq) C NCLEND
do i1,nq write (,"(i5, 2f10.3)")
i, x(i), q(i) end do return end
STEP 2 quad.so WRAPIT quad.f
STEPS 3-4 external QUPR
"./quad.so" begin a 2.5 b -.5 c
100. nx 10 x fspan(1., 10., 10) q
new (nx, float) QUPRcquad(a,b,c,
nx, x,q) QUPRprntq (x, q, nx) end
72
Example Linking F90 routines
quad.f90 subroutine cquad(a,
b, c, nq, x, quad) implicit none
integer , intent(in) nq real ,
intent(in) a, b, c, x(nq) real ,
intent(out) quad(nq) integer
i ! local quad
ax2 bx c ! array return end
subroutine cquad subroutine prntq(x, q, nq)
implicit none integer , intent(in) nq
real , intent(in) x(nq), q(nq)
integer i !
local do i 1, nq write (,
'(I5, 2F10.3)') i, x(i), q(i) end do
return end
STEP 1 quad90.stub C NCLFORTSTART
subroutine cquad (a,b,c,nq,x,quad)
dimension x(nq), quad(nq) ! ftn default C
NCLEND C NCLFORTSTART subroutine prntq
(x, q, nq) integer nq real x(nq),
q(nq) C NCLEND
prntq_I.f90 module
prntq_I interface subroutine prntq
(x,q,nq) real, dimension(nq) x, q
integer, intent(in) nq end
subroutine end interface end module
cquad_I.f90 module cquad_I
interface subroutine cquad
(a,b,c,nq,x,quad) real, intent(in)
a,b,c integer, intent(in) nq
real,dimension(nq), intent(in) x
real,dimension(nq),intent(out) quad end
subroutine end interface end module
STEP 2 quad90.so WRAPIT pg
quad90.stub printq_I.f90 \
cquad_I.f90 quad.f90 STEP 3-4
same as previous ncl lt PR_quad90.ncl
73
Fortran vs. NCL string - character

• NCL (C) ?? Fortran string/character interchange
  problematical

C NCLFORTSTART subroutine csdemo
(string_in, string_out)
character() string_in ! passed FROM ncl
character8 string_out !
passed TO ncl C NCLEND print ,
string_in string_out "F-to-NCL"
return end
external DEMO "./csdemo.so" begin cstring
new (8, character) DEMOcsdemo("NCL-F" ,
cstring) stringc chartostring( cstring )
print ( stringc ) end

• Can NOT pass arrays of strings or characters

74
Linking Commercial IMSL (NAG,) routines
STEP 1 rcurvWrap.f C NCLFORTSTART
subroutine rcurvwrap (n, x, y, nd, b, s, st,
n1) integer n, nd, n1 real
x(n), y(n), st(10), b(n1), s(n1) C NCLEND
call rcurv(n,x,y,nd,b,s,st) ! IMSL
return end

STEP 2 rcurvWrap.so WRAPIT l mp L
/usr/local/lib64/r4i4 l imsl_mp rcurvWrap.f
external IMSL ./rcurvWrap.so begin x (/
0,0,1,1,2,2,4,4,5,5,6,6,7,7 /) y (/508.1,
498.4, 568.2, 577.3, 651.7, 657.0, 755.3 \
758.9, 787.6, 792.1. 841.4, 831.8, 854.7,
871.4 /) nobs dimsizes(y) nd 2
n1 nd1 b new ( n1,
typeof(y)) s new ( n1, typeof(y))
st new (10, typeof(y))
IMSLrcurvwrap(nobs, x, y, nd, b, s, st, n1) end
75
Accessing LAPACK (1 of 2)

• double precision LAPACK (BLAS) gt distributed
  with NCL
• explicitly link LAPACK lib via fortran
  interface WRAPIT
• eg subroutine dgels solves over/underdetermine
  d real linear systems

C NCLFORTSTART SUBROUTINE DGELSI( M, N,
NRHS, A, B, LWORK, WORK ) IMPLICIT NONE
INTEGER M, N, NRHS, LWORK DOUBLE
PRECISION A( M, N ), B( M, NRHS), WORK(LWORK) C
NCLEND C
declare local variables INTEGER
INFO CHARACTER1 TRANS TRANS
"N CALL DGELS(TRANS, M,N,NRHS,A,LDA,B,LDB,W
ORK,LWORK,INFO) RETURN END
WRAPIT L NCARG_ROOT/lib -l lapack_ncl
dgels_interface.f
76
Accessing LAPACK (2 of 2)
external DGELS "./dgels_interface.so NAG
example http//www.nag.com/lapack-ex/node45.html
These are transposed from the fortran example
A (/ (/ -0.57, -1.93, 2.30, -1.93, 0.15, -0.02
/), \ (4,6) (/ -1.28, 1.08, 0.24,
0.64, 0.30, 1.03 /), \ (/ -0.39,
-0.31, 0.40, -0.66, 0.15, -1.43 /), \
(/ 0.25, -2.14,-0.35, 0.08,-2.13, 0.50 /)
/)1d0 must be double dimA
dimsizes(A) N dimA(0) 4 M
dimA(1) 6 B (/-2.67
,-0.55 ,3.34, -0.77, 0.48, 4.10/)1d0
must be double
LAPACK wants 2D nrhs 1 B2
conform_dims ( (/nrhs,M/), B, 1 ) (1,6)
B2(0,) B lwork 500

allocate space work new ( lwork,
"double", "No_FillValue") must be
double DGELSdgelsiw(M, N, nrhs, A , B2,
lwork, work ) print(B2(0,0N-1))
77
Combining NCL and Fortran in C-shell
!/usr/bin/csh NCL cat
gt! main.ncl ltlt "END_NCL" load
NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl
" load NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_c
sm.ncl load "NCARG_ROOT/lib/ncarg/nclscripts/csm
/contributed.ncl external SUB "./sub.so" begin
... end "END_NCL" FORTRAN
cat gt! sub.f ltlt "END_SUBF" C
NCLFORTSTART ... C NCLEND "END_SUBF"
WRAPIT WRAPIT sub.f
EXECUTE ncl main.ncl gt!
main.out
78
Global Variables and Scope 1 of 2

• Global Variable(s)
• by definition can be accessed from any function
  or procedure
• different from local variables
• NCL does not have explicit global variables
• requires understanding of NCLs variable scope
  identical to Pascal
• http//www.ncl.ucar.edu/Document/Manuals/Ref_Manu
  al/NclStatements.shtmlScoping

load dummy_1.ncl not aware of constants
below GRAVITY 9.8 RGAS
204 load dummy_2.ncl can use GRAVITY and
RGAS REARTH 6371.009 km load
dummy_3.ncl can use GRAVITY, RGAS,
REARTH begin can
use GRAVITY, RGAS, REARTH end
79
Global Variables and Scope 2 of 2

• knowledgeable user can simulate one approach
• create a file GLOBAL.ncl (no
  begin / end )
• populate with desired constants
• best to follow some user defined conventions
  e.g. capital letters

contents of GLOBAL.ncl GRAVITY 9.8
GRAVITY_at_units m/s RDRY
286.9
J/(kg-K) REARTH 63.71.009
km GRAVITY_d 9.81d m/s
(double)
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_cod
e.ncl" load "NCARG_ROOT/lib/ncarg/nclscripts/csm
/gsn_csm.ncl" load "NCARG_ROOT/lib/ncarg/nclscri
pts/csm/contributed.ncl" load "/my/path/GLOBAL.nc
l load "foo_1.ncl begin optional end
only if begin is preent
80
list variable 1 of 2

• container for variables of any type
• can be mixed type (numeric, string, file, )
• three ways to create a variable of list
• f addfiles()
• x // list constructor
• a NewList()

• q addfiles (fil_names, "r") q is of type
  list
• x q-gtX mean all elements
  of list
• s q282-gtX extract x from files
  2,4,6,8

• f addfile(., r) type file
• i ispan(1,10,2) integer
• z (/ example, wkshop/) array constructor
  variable (/.../)
• delete( / f, i, z / ) delete
  multiple variables

81
list variable 2 of 2

• r random_normal(mean, stdev, N)
• s (/ sample , string /)
• lst NewList (fifo) create a list
  variable
• ListPush (lst, r ) add r
• ListPush (lst, r2 ) add r2
• ListPush (lst, s ) add s
• nlst ListCount( lst ) number of items in
  lst
• do n0,nlst-1
• print( lstn ) print
  contents of list item n
• end do

82
NCL as a scripting tool
begin mssi getenv (MSSOCNHIST) get
environment variable diri /ptmp/user/
dir containing input
files fili b20.007.pop.h.0
prefix of input files diro
/ptmp/user/out/ dir
containing output files filo
b20.TEMP. prefix of
output files nyrStrt 300
1st year nyrLast 999
last year do
nyearnyrStrt,nyrLast print (----
nyear ----)
read 12 months for
nyear msscmd msrcp n mss mssi
filinyear -0-10-9.nc diri.
print (msscmdmsscmd) system
(msscmd)
strip off the TEMP variable
ncocmd ncrcat v TEMP
dirifili.nc dirofilonyear.nc
print (ncocmdncocmd) system
(ncocmd)
remove the 12 monthly files
rmcmd /bin/rm dirifilinyear
.nc print (rmcmdrmcmd)
system (rmcmd) end do end

• http//www.ncl.ucar.edu/Applications/system.shtml