中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Introduction（day 1） Language Basics（day 1） Input / Output（day 1）
11/02，2004@NCU 11/02， 11/04，2004@CCU 11/04， 3/14-16，2005@CWB 3/14- 16， 10/01 ，2005@NTU 7/27，2009@CCU 7/27，
Functions and Statistical methods（day 2） Graphics（day 3）
Update：2009.07.23
Function
average、summary、root-mean-square average、 summary、 root- meanstandard deviation、variance、 running average deviation、 variance、 climatology、seasonal mean climatology、
Statistic method
Correlation、Regression、Spectral、FFT、EOF Correlation、 Regression、 Spectral、 FFT、 SVD、Wavelets、Filter、T-Test、Taylor Diagrams SVD、 Wavelets、 Filter、 Test、
Special
vertical integral、interpolation、gradient、dtrend integral、 interpolation、 gradient、 、moisture、Laplacian、9 point smoothing、Vorticity moisture、 Laplacian、 smoothing、 、Divergence、stream function、velocity potential Divergence、 function、
Why Learn NCL
※ Integrated processing environment
pdf png
(Data)
jpg gif
(Figures) (Datas) Datas)
※ freeware: supported public domain software ※ portable: Linux/Unix, Windows, MacOSX ※ excellent graphics
Tu＠ＣＣＵ
1/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Software and Paths
● Binaries … Free （Download） ※ https://www.doczj.com/doc/c112425665.html,/
● NCAR PATH ( C-Shell )
※ setenv NCARG_ROOT /usr/local/lib/ncarg_4.3.1 ※ set path = ( /usr/local/lib/ncarg_4.3.1/bin )（自訂）
Useful NCL related URLs
● CSM and NCL links
https://www.doczj.com/doc/c112425665.html,
https://www.doczj.com/doc/c112425665.html,
● NCL functions and procedures
https://www.doczj.com/doc/c112425665.html,/Document/Functions/list_alpha.shtml https://www.doczj.com/doc/c112425665.html,/Document/Functions/list_alpha.shtml
● NCL Graphics Example Pages
https://www.doczj.com/doc/c112425665.html,/Applications/index.shtml https://www.doczj.com/doc/c112425665.html,/Applications/index.shtml
Tu＠ＣＣＵ
2/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Functions and Procedures
Color Table
NCL commitment to help !
ncl-talk@https://www.doczj.com/doc/c112425665.html, [cc: haley@https://www.doczj.com/doc/c112425665.html, , shea@https://www.doczj.com/doc/c112425665.html,] https://www.doczj.com/doc/c112425665.html,/mailman/listinfo/ncl-talk
Tu＠ＣＣＵ
3/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
NCL
? interpreted language – no optimization ? subscripts are 0 based ? row major ? right index fastest varying ? \ continue next line ? ; comment ? + string concatenation ? built in file handling for nc, grb, ccm, hdf, (hdfeos )
Fortran {C}
? compiled language – high optimization {good} ? subscripts are 1 based {0} ? column major {row} ? left index fastest varying {right} ? col 6 [f77]; ; & [f90] continue ? c [f77]; ![f90] comment { ; } ? // char concatenation {strcat}
Demo Script: Input/Function/Graphics
NCL
dat ( N , M ) dat ( 0 , 0 ) dat ( 0 , 1 ) dat ( 0 , 2 ) dat ( 1 , 0 ) dat ( 1 , 1 ) dat ( 1 , 2 ) <==> <==> <==> <==> <==> <==> <==>
Fortran {C}
dat ( M , N ) dat ( 1 , 1 ) dat ( 2 , 1 ) dat ( 3 , 1 ) dat ( 1 , 2 ) dat ( 2 , 2 ) dat ( 3 , 2 )
dat ( lev , lat , lon ) <=map=> dat ( lon , lat , lev )
load load load load
" $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl " " $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl " " $NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl“ " $NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
input = addfile("/ptmp/shea/TMP_1958-1997.nc", "r") ; open netCDF file dat = input->T ; dat (time,lev,lat,lon) wks = gsn_open_wks("pdf", "demo") gsn_define_colormap(wks,"gui_default") hov = dim_avg(dat ( lev|: , lat|: , lon|: , time|: )) res = True res@cnFillOn = True res@gsnSpreadColors = True res@gsnMaximize = True res@cnLineLabelFont = 21 ; open a graphic file ; color map ; function ; change default plot ; use colors ; use entire color map ; max plot size ; Label font
（屬性設定）
plot = gsn_csm_contour_map_ce(wks, hov(3 , {-60:60} , {100:300} ), res )
end
Executing NCL
? Interactive Mode (Command line)
prompt> ncl
ncl> record (“file_name”) ncl> stop record
Symbols
; @ $ (/.../) ! & {...} : | \ -> 註解符號（等同於fortran77 的 C 或 fortran90 的 ！） 註解符號（等同於fortran77 設定參數屬性或大小（ ex：res@vpXF = 0.12 ） ex： res@ 沿用先前已設定過的字串（ var = $vname+”-”+$vyr） $vname+” +$vyr） 設定陣列大小（ hov = new((/ny,nx/),"float")） new((/ny,nx/),"float")） 設定陣列維度名稱（ex：hov!0 = “lat” , hov!1 = “lon”） 設定陣列維度名稱（ex： ov! lat” hov! lon” 設定座標系統下各個維度相對應的數值。(ex：hov&lat = rlat) 設定座標系統下各個維度相對應的數值。(ex： ov& 在座標系統下設定座標範圍 (ex: hov=hov({-20:20},{90:200})) hov=hov({設定陣列中每個維度的範圍與間隔 ( hov（20:90:2 , ::5) ) hov（ 用於區隔維度名稱與維度 ( hov(:)=dim_avg(dt(lat|:,lon|:)) ) hov(:)=dim_avg(dt(lat|:,lon|:)) 連續符號（相當於fortran第六格的作用，置於最尾端） 連續符號（相當於fortran第六格的作用，置於最尾端） 用於資料的 I/O ( dtmp = input->olr ) (前後不能有空格). input(前後不能有空格).
? Batch Mode [ < and .ncl are optional ] – prompt> ncl < script.ncl – prompt> ncl script.ncl [also acceptable] – prompt> ncl < script.ncl > &! script.out – prompt> ncl < script.ncl > &! script.out & ※ scripts require begin………end statements ※ appending "&" means put in background ※ note: the >&! & are appropriate for csh and tcsh
Tu＠ＣＣＵ
4/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Manual Array Creation
? array
Automatic Array Creation
?
constructor characters (/…/) 【類似 fortran的data宣告 】
– a_integer – a_float
= (/1,2,3/) = (/1.0, 2.0, 3.0/) , a_double = (/1., 2, 3.2d0 /) – a_string = (/"a","b","c"/) – a_logical = (/True, False,True/) – a_2Darray = (/ (/1,2,3/), (/4,5,6/), (/7,8,9/) /) = 3 x 3的矩陣
variable to variable assignment
– –
y=x y => same size, type as x plus meta data no need to preallocate space for y
?
data importation via supported format
– –
? new function [Fortran dimension, allocate]
– x = new (array_size, type, missing value)（missing value is optional）
– a = new(3, float) – b = new(10, double, 1d+20) – c = new( (/ 5, 6, 7 /), integer) – d = new(dimsizes(U), string) – e = new(dimsizes(ndtooned(U)), logical)
u = input->U same for subset of data: u = f->U(:, 3:9:2, :, 10:20) – meta data (coordinate will reflect subset)
?
functions
return array: no need to preallocate space T42 = f2gsh ( gridi, (/ 64,128/), 42) – gridi(10,30,73,144) T42(10,30,64,128)
– –
? new and (/…/) can appear anywhere in script
T42 = f2gsh_Wrap ( gridi, (/ 64,128/), 42) ; contributed.ncl
delete
Deletes variables, attributes, and coordinate variables. variables,
?
Attributes
assign/access with 「@」 character
T@long_name T@wgts – T@x3d
– –
delete ( data )
Example 1 （integer to float） float） x = (/1,2,3,4,5/) delete(x) ; The delete is necessary because delete(x) x = (/1.,2.,3.,4.,5./) ; was originally an integer array. Example 2 （size remodel） remodel） x = (/1,2,3,4,5/) delete(x) ; The delete is necessary because delete(x) x = (/1,2,3,4/) ; was originally an array of 5 elements.
= “temperature” = (/ 0.25, 0.5, 0.25 /) = x3D = = = = = delx dely cmin cmax cint
res@vpWidthF res@vpHeightF res@cnMinLevelValF res@cnMaxLevelValF res@cnLevelSpacingF
?
delete can eliminate an attribute
–
delete(T@long_name)
Dimensions ( GrADS )
XDEF 128 LINEAR 0 2.8125 YDEF 64 LEVELS -87.864 -85.097 -82.313 -79.526 -76.737 -73.948 -71.158 -68.368 -65.578 -62.787 -59.997 -57.207 -54.416 -51.626 -48.835 -46.045 -43.254 -40.464 -37.673 -34.883 -32.092 -29.301 -26.511 -23.720 -20.930 -18.139 -15.348 -12.558 -9.767 -6.977 -4.186 -1.395 1.395 4.186 6.977 9.767 12.558 15.348 18.139 20.930 23.720 26.511 20.930 29.301 32.092 34.883 37.673 40.464 43.254 37.673 46.045 48.835 51.626 54.416 57.207 59.997 54.416 62.787 65.578 68.368 71.158 73.948 76.737 71.158 79.526 82.313 85.097 87.864 87.864 ZDEF 12 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 TDEF 5000 LINEAR 15jan1860 1mon
Name Dimensions ( NCL)
? May be “named”
陣列的維度由左至右從零開始（0、1、2、……… ）依序編號 ex：T (nz,ny,nx) 【 T(nx,ny,nz) → fortran】
? 設定維度名稱需採用符號「!」 {let T(nz,ny,nx)} 定義完各個座標名稱之後，我們就可以將陣列的維度 – T!0 = "lev" 編排方式依照個人需求重新排列，而座標變數也可以 – T!1 = "lat " 用來代表座標軸的index，例如： reordered_T = T（lon|:,lat|:,lev|:） – T!2 = "lon" ? assign dimension name(s) → functions dd= nameDim(dd , dimNames[*] , longName , units)
【T = nameDim (T, (/ "lev" , "lat " , "lon" /), "Temperature", "C") 】 【u = nameDim (u, (/"time","lev","lat","lon"/), "zonal wind","m/s"】
Tu＠ＣＣＵ
5/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Longitude Manipulators
T ( nz , ny , nx ) -> Tnew ( nx , ny , nz )
經度
rlon = lonGlobeF (nx, "lon", "longitude", "degrees_east") nx, "lon", "degrees_east") lon will run from 0 to （360－rint） 360－ rint）
Fortran
【註】如果nx = 144，那麼經度範圍將會介於 0 到 357.5 之間 如果nx 144，那麼經度範圍將會介於
rlon = lonGlobeFo (nx, "lon", "longitude", "degrees_east") nx, "lon", "degrees_east") lon will run from（rint/2）to（360－rint/2） from（ rint/2） to（ 360－ rint/2）
【註】如果nx = 144，那麼經度範圍將會介於1.25 到 358.75 之間 如果nx 144，那麼經度範圍將會介於1.25
NCL
Latitude Manipulators 緯度
rlat = latGlobeF (ny, "lat", "latitude", "degrees_north") ny, "degrees_north") lat will run from 90S to 90N 【註】如果nx = 73，那麼緯度範圍將會介於 90S 到 90N 之間 如果nx 73，那麼緯度範圍將會介於 rlat = latGlobeFo(ny, "lat", "latitude", "degrees_north") latGlobeFo(ny, "degrees_north") lat will run from （－90＋rint/2）to（90－rint/2） （－90＋rint/2） to（ 90－ rint/2）
【註】如果nx = 72，那麼緯度範圍將會介於 88.75S 到 88.75N 之間 如果nx 72，那麼緯度範圍將會介於
rlon = fspan(0 , 360 *(nx-1) /nx , nx) *(nx- /nx nx) rlon = (/ rlon-180. /) rlon; grid goes DateLine eastward
Lon /Lat/Weights
權重
gwts = latGauWgt (ny, "lat", "gaussian weights", "dimension_less") ny, "gaussian "dimension_less") Generate gaussian weights and meta data gwts = NormCosWgtGlobe(rlat) NormCosWgtGlobe(rlat) normalize the cosine wgts gwts = SqrtCosWgtGlobe(rlat) SqrtCosWgtGlobe(rlat) 也可以使用 gau_info = gaus(ny/2) rlat = gau_info(:,0) gau_info(:,0 ; divide by 2 to get "per hemisphere" 或 gau_info = doubletofloat(gaus(nlat/2)); convert double the float doubletofloat(gaus(nlat/2)); ; gaussian latitudes ( 1st dimension of gau_info) gau_info) ; gaussian weights ( 2nd dimension of gau_info) gau_info) gwts = gau_info(:,1) gau_info(:,1
rlat = latGau
(ny, "lat", "latitude", "degrees_north") ny, "degrees_north")
Generates gaussian latitudes and associated meta data. rlat = rlat(::-1) rlat(::; grid goes from North → South
Create and Assign Coordinate Variables
geographical coordinate arrays with named dimensions and units
named dimensions
Latitude LATITUDE Longitude LONGITUDE
?
create vector array
U = new ((/ 12 , 73 , 144 /) , " float " , U@_FillValue) rlon = lonGlobeF (nx, "lon", "longitude", "degrees_east") rlat = latGlobeF (ny, "lat", "latitude", "degrees_north") plev = (/1000,925,850,700,600,500,400,300,250,200,150,100/) plev皆需給定單位） plev@units = "hPa" （rlon、rlat、plev皆需給定單位）
units
degrees_north degrees-north degree-north degrees north degrees_N Degrees_east degrees_east degrees-east degree-east degrees east degrees_E Degrees_east ?
lat Lat LAT hlat lat_u lat_t
lon Lon LON hlon lon_u lon_t
?
assign dimension name(s) → functions
U = nameDim (U , (/ "lev " , "lat " , "lon" /), "zonal wind", "m/s")
assign coordinate variables to U → via『&』
U&lon = rlon U&lat = rlat – U&lev = plev
–
–
Tu＠ＣＣＵ
6/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Create and Assign Coordinate Variables 假設 SST 是一個二維陣列，大小為（nlat, nlon），那麼 是一個二維陣列，大小為（nlat, nlon），那麼 要定義維度的名稱，可以使用符號「!」，如果要將座標 要定義維度的名稱，可以使用符號「! 軸相對應的數值帶入，則可使用符號「&」。 軸相對應的數值帶入，則可使用符號「& SST!0 SST!1 rlon rlat rlon@units rlat@units SST&lon SST&lat = = = = = = = = "lat" "lon" lon" fspan(0.,355.,72) fspan(-90.,90.,37) fspan( "degrees_east" degrees_east" "degrees_north" degrees_north" rlon rlat
Variable Subscripting (1 of 2)
※ Standard subscripting（標準網格註腳， similar to f90 and f77）
每個註腳的維度開始值從0 到 ( N-1 ) ，由右至左 ( nt, ny, nx ) ，維度之間用逗號 ( , )區隔。標準的註腳寫法 → 【開始值 : 結束值 : 間隔大小】亦可簡化成【：】 T T( 0 , : , : : 5 ) T(0 , : : -1 , : 50 ) T( : 1 , 45 , 10 : 20 ) entire array [ don't use T(:,:,:) ] 1st time index, all lat, every 5th lon 1st time index, reverse lat order,1st 51 lon 1st 2 time indices, 46th index of lat, 10-20 indicies of lon
※ Coordinate subscripting（加註座標軸的註腳，similar to GrADS）
資料陣列中，每個維度都可用一個一維陣列加以表示，數值大小可以依序增加， 也可逐漸減小，座標標註可用 『 {...}』加以表示。假設T(nt,ny,nx) dat = T( : , {-30:30} , : ) dat = T( : , {-20:20} , {90:290:2} ) 取所有時間和經度，緯度從30S到30N的資料。 取所有時間, 緯度從20S到20N 經度從90到290，每間隔2取一筆資料。
Variable Subscripting (2 of 2)
※ Named / Coordinate subscripting（加註座標名稱）
indicated by 「 | 」（維度名稱與座標位置間的區隔符號）
Array Dimension Reduction
?
let T ( 12 , 64 , 128 )
– Tjan = T(
0, : , : )
Tjan(64,128)
假設一個陣列 T(nz,ny,nx) T(nz,ny,nx)
T rlon rlat plev = nameDim (T , (/ "lev " , "lat " , "lon" /), “Temperature", "m/s") "lev "lon" "m/s") = lonGlobeF(nx,"lon","longitude","degrees_east") lonGlobeF(nx,"lon","longitude","degrees_east") = latGlobeF(ny,"lat","latitude","degrees_north") latGlobeF(ny,"lat","latitude","degrees_north") = (/1000,925,850,700,600,500,400,300,250,200,150,100/)
– Tjan automatically becomes 2D – all applicable meta data copied
plev@units = "hPa" "hPa" T&lat = rlat , T&lon = rlon , T&lev=plev lev= dd = T(lat|:, lon|:, lev|:) T(lat|:, lon|:, lev|:) dd = T(lev|:,{lon|90:120},{lat|-20:20}) T(lev|:, lon| 90:120} lat| 20:20} dd = T(:,29:45,{lon|90:120}) T(: ,29: 45,{ lon| 90:120} makes dd (lat,lon,lev) lat,lon,lev)
all levels , 90°E～120°E , 20°S～20°N 90° 120° 20° 20° all levels , 90°E～120°E , 20°S～20°N 90° 120° 20° 20°
?
can override dimension reduction
– Tjan – Tjan
= T( 0 : 0 , : , : )
Tjan( 1 , 64 , 128 ) Tjan(
= new( (/1,64,128/), typeof(T), T@_FillValue) new( (/1,64,128/), typeof(T), T@_FillValue)
Tjan(0,:,:) = T(0,:,:) jan(0,:,:)
Standard and Coordinate Subscripting
_FillValue
? most NCL functions ignore _FillValue – “missing_value” has no special status to NCL missing_value” – if “T” has “missing_value” but no “_FillValue ” missing_value” ? T@_FillValue = T@missing_value – T@_FillValue = -99.99
Latitude coordinate variable (1D)
named/coordinate subscripting can be mixed
Standard
T(8:15,1:12)
Coordinate
T({-30:30},{-180:-90})
Combined
Longitude coordinate variable (1D) T({-30:30}, 1:12)
Tu＠ＣＣＵ
7/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Expressions
Algebraic operators
＋ －
* / ^ >, < 加（除了數字加減外，也可用於文字） 減 乘 除（除法的運算速度較慢，所以建議改用乘） 指數（10^2 = 102） 大於, 小於
Data Types
? ? ? ? ? ? ?
?
algebraic operator
5.3 + 7.95 5.0 * 3.0 13.25 15.0
string concatenator
"alpha" + (5.3 + 7) alpha12.3 "A="+3+" B="+5 A=3 B=5
.lt. .le. .gt. .ne. .eq. .and. .xor. .or. .not.
Logical
小於 小於等於 大於 不等於 等於 and exclusive-or or Not
character float (32 bit) integer (32 bit) double (64 bit) long (32 bit) [SGI 64] short (16 bit) byte ( 8 bit, 0-255) complex NOT supported
Do Loops
do n = start , end , skip [statements] end do do while (expression) [statements] end do break: jump to first statement after end do [f90: exit] continue: proceed to next iteration [f90: cycle]
If Statements
if (scalar_logical_expression) then [statements] end if if (scalar_logical_expression) then [statements] else [statements] end if ※ no “else if ” ※ no if (logical_expression) [statements]
Minimize Use of Loops (1 of 2)
? parameters for each example: – nx=1000, ny=1000 – T=new( (/ny,nx/),float)
; 136 seconds
※ “end if ”, “end do” ※ do loop 運算速度非常緩慢，盡量避免使用。
T (i , j ) = 100 ? (8 i ) + (8 j ) do j= 1,ny do i= 1,nx T( j-1 , i-1 ) = 100.0 - sqrt((8.0*i)^2 + (8.0*j)^2) end do ; 88 seconds end do do j= 1 , ny do i= 1 , nx T(j-1,i-1) = i^2 + j^2 end do end do T = 100.0 - sqrt((8.0^2) * T)
2
2
Minimize Use of Loops (2 of 2)
; <1 sec---------ispn = ispan( 1 , nx , 1 )^2 jspn = ispan( 1 , ny , 1 )^2 do i= 0,dimsizes(ispn)-1 T(:,i) = ispn(i) + jspn end do T = 100.0 - sqrt((8.0^2) *T )
?
– –
print (1 of 3)
Prints out all variable information including
meta data, values T(lat,lon): print (T)
Fortran?
; <1 sec
T = onedtond(ispan(1,nx,1),(/ny,nx/)) T!0 = "x" T!1 = "y" T = 100.0 - sqrt((8.0*T(:,:))^2 + (8.0*T(y | :, x | :))^2)
Variable: T Type: float Total Size: 32768 bytes 8192 values Number of Dimensions: 2 Dimensions / Sizes: [lat | 64] x [lon | 128] Coordinates: lat: [-87.86379 .. 87.86379] lon: [ 0. 0 .. 357.1875] Number of Attributes: 2 long_name: Temperature units: C (0,0) -31.7 (0,1) -31.4 (0,2) -32.3 (0,3) -33.4 (0,4) -31.3 etc. [entire T array will be printed]
Tu＠ＣＣＵ
8/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
print (2 of 3)
?
print (3 of 3)
?
can be used to print a subset of array
– –
print with embedded strings
– –
meta data, values T(nlat,nlon): print( T(:,103) ) or print( T(:,{110}) )
Variable: T (subsection) Type: float Total Size: 256 bytes 64 values Number of Dimensions: 1 Dimensions / Sizes: [lat | 64] Coordinates: lat: [-87.86379 .. 87.86379] Number of Attributes: 3 long_name: Temperature units: C lon: 109.6875 [ added ] (0) -40.7 (1) -33.0 (2) -25.1 (3) -20.0 (4) -15.3 etc.
no meta data print ( "min(T)="+min(T)+" max(T)="+max(T) ) (0) min(T)=-53.8125 max(T)=25.9736
? sprintf
and sprinti provide formatting
- often used in graphics - print ( "min(T) = "+ sprintf("%5.2f ", min(T)) )
? sprinti can left fill with zeros
(ex: let n=3)
- fnam = "h" + sprinti ("%0.5i", n) + ".nc" - print("file name = "+fnam)
printVarSummary
? Print out variable (data object) information
– – – – type dimension information coordinate information (if present) attributes (if present)
? printVarSummary (u)
Variable: u Type: double Total Size: 1179648 bytes 147456 values Number of Dimensions: 4 Dimensions / Sizes: [time | 1] x [lev | 18] x [lat | 64] x [lon | 128] Coordinates: time: [4046..4046] lev: [4.809 .. 992.5282] lat: [-87.86379 .. 87.86379] lon: [ 0. 0 .. 357.1875] Number of Attributes: 2 long_name: zonal wind component units: m/s
Supported formats in NCL
● ASCⅡ ASCⅡ ● Fortran Binary ( sequential / direct ) [open(11,form="unformatted", access="sequential")] [open(11,form="unformatted",access=direct",recl=_)] [open(11,form="unformatted",access=direct",recl=_)] ● NetCDF (Network Common Data Format) Network ● GRIB (Grid in Binary; WMO standard) Grid ● HDF (Hierarchical Data Format (Scientific Data Set only) ) ● HDF-EOS (Earth Observing System) HDF● CCMHT (CCM History Tape) History
Tu＠ＣＣＵ
9/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Reading / Writing ASCII files
Reading
remove any previously exist file
head = readAsciiHead (fname:string , opt:numeric) opt:numeric) data = asciiread (fname:string , dimensions , datatype:string) datatype:string) data = readAsciiTable (fname:string , ncol:integer , datatype:string , opt)
system (path+"/rm "+ofn)
Writing
asciiwrite (fname:string , dat) dat) write_matrix → 2D
number of rows and columns
行 列 nrow = numAsciiRow ( ifn : string) ncol = numAsciiCol ( ifn : string)
ASCII file：Time Series Data Assume you have the file "test.dat" "test.dat"
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
DARWIN SEA LEVEL PRESS
YEAR JAN FEB MAR APR MAY JUN JUL AUG SEP SEP OCT NOV DEC
1951 1952 1953 1954 1955 …… [snip]
5.3 5.1 6.7 7.1 0.4 -1.6 1.0 -1.2 -1.1 2.9
7.0 9.4 11.9 11.2 13.3 13.2 12.5 13.2 8.3 10.3 10.3 12.5 13.0 12.6 12.1 -1.4 -0.1 -3.6 -0.5 -0.2 -3.1 -2.4 -0.5 0.6 0.4 -0.5 0.4 1.3 0.3 0.1 -0.8 1.4 1.7 2.7 2.0 2.5
11.4 9.9 -0.3 0.1 2.5
9.5 7.7 -0.7 0.2 2.2
8.9 8.4 -1.1 2.4 1.7
【方法一】 方法一】
soi = new((/nt,nm/),"float") new((/nt,nm/),"float") head = readAsciiHead( ifn , 2 ) readAsciiHead( head = readAsciiHead( ifn , "YEAR" ) readAsciiHead( "YEAR" SOI = asciiread(ifn , (/nt,nm/),"float" ) (/nt,nm/),"float"
【方法二】 方法二】
nm = 13 soi = new((/nt,nm/),"float") new((/nt,nm/),"float") SOI = readAsciiTable(ifn, nm, "float", 2 ) readAsciiTable(ifn, nm, SOI = readAsciiTable(ifn, nm,"float","YEAR") readAsciiTable(ifn, nm,"float","YEAR")
行 列
ncol = numAsciiCol("test.dat") numAsciiCol("test.dat") nrow = numAsciiRow("test.dat") numAsciiRow("test.dat")
由於不需要使用do loop讀取資料，所以必須給定剛好的時間筆數 由於不需要使用do loop讀取資料，所以必須給定剛好的時間筆數
(0,0)
YEAR MON NINO1+2 ANOM NINO3 NINO1+2 NINO3 ANOM NINO4 NINO4 ANOM NINO3.4 NINO3.4 ANOM
1 2 3 4
2004 1 2004 2 2004 3 2004 4 2004 5 2004 6 2004 7 2004 8 2004 9 2004 10 2004 11 2004 12
【方法一】
nn = 10
23.49 24.56 25.59 23.84 23.41 21.70 20.74 20.52 19.82 20.19 20.46 21.92
-1.00 -1.49 -0.89 -1.62 -0.93 -1.32 -1.11 -0.30 -0.65 -0.72 -1.22 -0.92
23.97 24.51 26.65 26.65 25.91 25.64 25.17 24.51 24.12 24.32 23.80 24.44
-1.64 -1.85 -0.43 -0.76 -1.14 -0.74 -0.41 -0.44 -0.71 -0.57 -1.15 -0.64
27.41 26.77 27.00 27.12 27.72 27.85 27.84 27.54 27.12 27.64 27.13 27.27
-0.73 -1.24 -1.09 -1.28 -0.92 -0.79 -0.73 -0.91 -1.37 -0.77 -1.23 -1.01
25.01 24.92 26.41 26.75 26.30 26.79 26.59 26.11 25.56 26.03 25.42 25.54
-1.50 -1.77 -0.73 -0.94 -1.47 -0.70 -0.48 -0.59 -1.09 -0.57 -1.09 -0.93
(0,1) (0,2) (0,3)
【方法一】 SST = new((/nt,nn/),"float") head = readAsciiHead( ifn , 1 ) SST = asciiread(ifn , (/nm,nn/),"float" ) SST = readAsciiTable(ifn, nn, "float", 1)
YEAR MON NINO1+2 ANOM NINO3 ANOM NINO4 ANOM NINO3.4 ANOM
(0,4) 3.4 (0,5) 2004 (0,6) 1 (0,7) 23.49 (0,8) -1 (0,9) 23.97 (1,0) -1.64 (1,1) 27.41 (1,2) -0.73 (1,3) 25.01 (1,4) -1.5
【方法二】
nn = 10 nm = 12 SST = new((/nm,nn/),"float") SST = readAsciiTable(ifn, nn, "float", 1) SST = readAsciiTable(ifn, nn,"float","YEAR")
SST = new((/nm,nn/),"float") head = readAsciiHead( ifn , 1 ) head = readAsciiHead( ifn , "YEAR" ) SST = asciiread(ifn , (/nm,nn/),"float" )
2004 1 2004 2 2004 3 2004 4 2004 5 2004 6 2004 7 2004 8 2004 9 2004 10 2004 11 2004 12
23.49 24.56 25.59 23.84 23.41 21.70 20.74 20.52 19.82 20.19 20.46 21.92
-1.00 -1.49 -0.89 -1.62 -0.93 -1.32 -1.11 -0.30 -0.65 -0.72 -1.22 -0.92
23.97 24.51 26.65 26.65 25.91 25.64 25.17 24.51 24.12 24.32 23.80 24.44
-1.64 -1.85 -0.43 -0.76 -1.14 -0.74 -0.41 -0.44 -0.71 -0.57 -1.15 -0.64
27.41 26.77 27.00 27.12 27.72 27.85 27.84 27.54 27.12 27.64 27.13 27.27
-0.73 -1.24 -1.09 -1.28 -0.92 -0.79 -0.73 -0.91 -1.37 -0.77 -1.23 -1.01
25.01 24.92 26.41 26.75 26.30 26.79 26.59 26.11 25.56 26.03 25.42 25.54
-1.50 -1.77 -0.73 -0.94 -1.47 -0.70 -0.48 -0.59 -1.09 -0.57 -1.09 -0.93
Tu＠ＣＣＵ
10/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Writing ASCII files
ifn ofn nx ny dat = "/glacier2/tu/data/xy.asc" = "new.asc" "new.asc" = 144 = 73 = asciiread ( ifn , (/ny,nx/) , " float " ) (/ny,nx/)
write_matrix(dat[*][*], format, opt) write_matrix(dat[*][*], format,
pretty-print 2D array to standard out prettyinteger, float, double user format control (format) T(7,5): write_matrix (T, " 5f7.2 ", False)
4.35 4.36 9.73 4.91 0.17 -0.63 -4.29 4.39 4.66 -5.84 4.59 3.68 -4.12 0.07 0.27 3.77 0.89 -3.09 5.08 -2.51 5.85 -3.35 -1.66 8.46 7.55 0.14 1.76 0.87
在螢幕上顯示結果
-6.90 4.06 10.39 4.56 -5.63 -1.43 8.65
asciiwrite (ofn, dat) ofn, dat) 0.8400397 0.55219 0.2622927 -0.02875445 ………
? can also create an ASCII file
Example
(STAND TAHITI - STAND DARWIN) SEA LEVEL PRESS YEAR JAN FEB MAR APR MAY JUN JUL AUG 1951 2.7 1.0 -1.3 -1.1 -1.7 -0.5 -2.3 -1.2 1952 -2.0 -1.8 0.0 -0.9 1.0 0.8 0.7 -0.7 1953 0.4 -1.6 -1.4 -0.1 -3.6 -0.5 -0.2 -3.1 1954 1.0 -1.2 -0.5 0.6 0.4 -0.5 0.4 1.3 1955 -1.1 2.9 0.1 -0.8 1.4 1.7 2.7 2.0
; #################################################
load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
Reading/Writing Binary files
Direct
SEP -2.1 -0.4 -2.4 0.3 2.5 OCT -2.3 0.3 -0.3 0.1 2.5 NOV -1.6 -0.3 -0.7 0.2 2.2 DEC -1.6 -2.6 -1.1 2.4 1.7
specified direct record
data = fbindirread (ifn,rec,dim,type) fbindirwrite (ofn,data) data = fbinrecread (ifn,rec,dim,type) fbinrecwrite (ofn,rec,data) data = craybinrecread (ifn,rec,dim,type) fbinwrite (ofn,data) data = cbinread (ifn,dim,type) cbinwrite (ofn,data)
」
Sequential
multiple unformatted sequential records
Cray Sequential
single unformatted sequential record
; ################################################## begin 0.00 -0.90 ifn = "aa.asc" nm = 13 -1.40 -0.10 nt = 5 -0.50 0.60 soi = new((/nt,nm/),"float")
1.00 0.80 0.70 -3.60 -0.50 -0.20 0.40 -0.50 0.40
Cray(C block IO write)
mimics a C block IO write
※ “rec” which start at 「
soi = readAsciiTable(ifn,nm,"float",2)
※ “type” : double、float、integer、byte……………. ※ 文字或資料型態設定採用雙引號。（var=“sst” ， “float” ， ……………） ※ sequential中的“rec”如果設為「 」，表依個人所要求之陣列大小輸出。
end
Reading NetCDF & Grib & HDF files sequential binary format
fbinrecwrite (path+ofn,-1, lat) fbinrecwrite (path+ofn,-1, lon) fbinrecwrite (path+ofn,-1, sst) lat=fbinrecread(ifn,0, ny, "float") lon=fbinrecread(ifn,1, nx, "float") sst=fbinrecread(ifn,2,(/ny,nx/),"float") real rlon(nx), rlat(ny), sst(nx,ny) open (11,file=ofn,form="unformatted" ,access="sequential") read (11) rlat read (11) rlon read (11) sst SST dat ifn ofn = = = =
direct binary format
new ((/nt,nm,ny,nx/),"float" ) new ((/nm,ny,nx/),"float" ) "OIMSST_V2_1982_2003_mon.bin" "OIMSST_V2_1982_2003_cli.bin"
※ 資料的副檔名必須依檔案類型作設定。 「NetCDF」→「nc」 或「cdf」 「GRIB」 →「grb」或「grib」 「HDF 」 →「hdf」或「hdfeos」
【 read / write】 【 read only】 【 read / write】
do l=0,nt-1 do m=0,nm-1
rec = m + l*nm SST( l , m , : , : ) =fbindirread(ifn,rec,(/ny,nx/),"float")
end do end do do m=0,nm-1 fbindirwrite (ofn,dat( m , : , : )) end do open (11,file=ofn, form="unformatted“ ,access="direct", recl=nx*ny) do m=1,nm read (11,rec=m) & ((dat(i,j,m),i=1,nx),j=1,ny) end do
※ 讀取相關類型資料之function input = addfile （fname:string , status:string）→ 讀取單個檔 input = addfiles（fname:string , status:string）→ 讀取多個檔 status = “r” （表示 read， 適用於所有格式之資料） status = “c” （表示 create，適用於NetCDF和HDF4） status = “w” （表示 write， 適用於NetCDF和HDF4）
Fortran
; output subsection (lat,lon) only do m=0,nm-1 fbinrecwrite(ofn,-1,dd(m,:,:)) end do do m=1,nm read(11)((dd(i,j,m),i=1,nx),j=1,ny) enddo
可以利用「print (input)」 ※ 範例 input= addfile("oisst.nc"," r ") 確認變數名稱 dims= filevardimsizes(input,"sst") dd = new((/dims(0),dims(1),dims(2)/),float) dd = input->sst
Tu＠ＣＣＵ
11/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
ncdump
dimensions: lon = 144 ; lat = 73 ; time = UNLIMITED ; // (12 currently) variables: float lon(lon) ; lon:units = "degrees_east" ; lon:long_name = "Longitude" ; lon:actual_range = 0.f, 360.f ; float lat(lat) ; lat:units = "degrees_north" ; lat:actual_range = 90.f, -90.f ; lat:long_name = "Latitude" ; double time(time) ; time:units = "hours since 1-1-1 00:00:0.0" ; time:long_name = "Time" ; time:actual_range = 0., 8016. ; time:delta_t = "0000-01-00 00:00:00" ; time:avg_period = "0000-01-00 00:00:00" ; time:ltm_range = 17338824., 17487096. ; short olr(time, lat, lon) ; olr:long_name = "OLR monthly means" ; olr:valid_range = 0.f, 500.f ; olr:actual_range = 134.9129f, 293.3665f ; olr:units = "W/m^2" ; olr:add_offset = 327.65f ; olr:scale_factor = 0.01f ; olr:missing_value = 32766s ; olr:var_desc="Outgoing Longwave Radiation\n", "O" olr:dataset = "NOAA Interpolated OLR\n",
ncl
input = addfile("olr.mon.ltm.nc","r") 【抓取資料維度（dimensions）】 dims = filevardimsizes(input,"olr") ; nlat = filevardimsizes(input,"lat") mlon = filevardimsizes(input,"lon") print(dims) Variable: dims Type: integer Total Size: 12 bytes 3 values Number of Dimensions: 1 Dimensions and sizes: [3] Coordinates: (0) 12 (1) 73 (2) 144 【讀variables之資料】 data = input->olr(0,:,:) rlon = input->lon rlat = input->lat time = input->time 【讀variables內之參數】 unit = input->olr@units add = input->olr@add_offset scale = input->olr@scale_factor spv = input->olr@missing_value
Reading NetCDF files
begin ; ####### path = "/glacier1/OLR/src/" ifn = "OLR-197406-200411-mon.nc" "OLR- 197406- 200411ofn = "NOAA_OLR_1979_2004.bin" ifn = path+ifn ; #######
unit = input->olr@units add = input->olr@add_offset scale = input->olr@scale_factor spv = input->olr@missing_value
; ####### dd = new((/dims(0),dims(1),dims(2)/),"float") ; #######
input = addfile(ifn,"r")
; #######
dd = input->olr input; #######
dims = filevardimsizes(input,"olr")
nt = dims(0) ny = dims(1) nx = dims(2) print(dims) print(dims)
do m=0,dims(0)-1 fbindirwrite(ofn,dd(m,:,:)) end do
; ####### end
begin ; ############## path = " /glacier1/ECMWF/ERA15/PRS/PRS_1990/ " ifn = "ERA15_T106_PRS_1990_mon.grb " input = addfile(path+ifn,"r") ; ############## printVarSummary(input) end
GRIB HDF
Variable: input (0) filename: ERA15_T106_PRS_1990_mon path: /glacier1/ECMWF/ERA15/PRS/PRS_1990/ERA15_T106_PRS_1990_mon.grb file global attributes: dimensions: initial_time0 = 12 lv_ISBL1 = 17 g4_lat_2 = 160 g4_lon_3 = 320 initial_time4 = 12 variables: float T_GDS4_ISBL_10 ( initial_time0, lv_ISBL1, g4_lat_2, g4_lon_3 ) center : European Center for Medium-Range Weather Forecasts - Reading long_name : Temperature units : K _FillValue : -999 level_indicator : 100 grid_number : 255 parameter_number : 130 forecast_time : 19
begin (reading ; ####### ifn = "ERA15_T106_PRS_1990_mon.grb" ofn = "ERA15_T106_T_1990_mon.bin" input = addfile(ifn,"r") ; ####### system("/bin/rm "+ofn) ; ####### dims = filevardimsizes(input," T_GDS4_ISBL_10 ") print(dims) ; ####### dd = new((/dims(0),dims(1),dims(2)/),"float") ; ####### dd = input->T_GDS4_ISBL_10 ; ####### do m=0,dims(0)-1 fbindirwrite(ofn,dd(m,:,:)) end do ; ####### end
grib)
get dimension sizes
Numerous functions to query contents of supported file
? ? ? ? ? ? ? ? ?
getfilevarnames dpath = "/fs/cgd/data0/shea/ccm/" file = "H200503" getfilevardims getfilevaratts getfilevartypes isfilevar isfilevaratt isfilevardim isfilevarcoord
vname = getfilevarnames (ifn) if (isfilevarcoord(ifn, "U", "lat")) then … end if ext ifn = ".ccm" = addfile(dpath+file+ext , " r ")
? returns the dimension sizes of a variable ? will return 1D array of integers if the array queried is multi-dimensional.
begin input=addfile("OLR.nc", "r") ; ####### dd =input->olr dims = dimsizes(dd) ; ####### dims = filevardimsizes(input,"olr") ; ####### print(dims) end
getfilevardimsizes
dimsizes (針對陣列) filevardimsizes (針對變數檔案)
Tu＠ＣＣＵ
12/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
get dimension sizes
begin input = addfile(“R2.nc”,”r”) t = input->T
Variable: dims 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) rank=4
addfiles
(1 of 2)
? spans multiple files fnames = (/ "reAnal1","reAnal2","reAnal3","reAnal4"/)
?
dims = dimsizes(t) print(dims) rank = dimsizes(dims) print ("rank="+rank) end
input = addfiles (fnames, "r") ? fnames is a 1D array of file names (strings) ? can be used for any supported format
T = input[:]->T ? read T [values only] from each file in list 'input' T will not have any meta data ? T must exist in each file and be same shape [rank] ? a list is used to sequence results of addfiles ? normal file variable selection is used with "[…]"
addfiles
?
? ?
(2 of 2)
begin
addfiles
fnames = (/ "reAnal1","reAnal2","reAnal3","reAnal4"/) input = addfiles (path+fnames, "r") ListSetType (input, "cat") ListSetType (input, "join") T end
Number of Dimensions: Dimensions and sizes: Number of Dimensions:
2 options on variable merging
ListSetType (a, "cat") ListSetType (a, "join”) [default; “cat” => concatenation]
T1(time,https://www.doczj.com/doc/c112425665.html,t,lon)
?
when to use "cat" and "join" [rule of thumb]
? ?
cat : continuous record join : creating ensembles a record dimension will be added
= input[:]->T
?
2 options on variable merging ListSetType (a, "cat") continuous record ListSetType (a, "join”) creating ensembles
cat join
?
suggest using addfiles_GetVar [contributed.ncl]
?attaches
meta data to variable
Dimensions and sizes:
Tu＠ＣＣＵ
13/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
abs( value:integer) fabs(value:numeric) sin(value:numeric) asin(value:numeric) cos(value:numeric) acos(value:numeric) tan(value:numeric) atan(value:numeric) atan2(y:numeric,x:numeric) tanh(value:numeric) exp(value:numeric) log(value:numeric) log10(value:numeric) sqrt(value:numeric) max(arg:numeric) min(arg:numeric)
→ computes absolute value for integer input （ abs(-2) = 2） → computes absolute value for numeric input （fabs(-1.2) = 1.2） → computes sin from radian input （sin(3.1415/180*45)=0.707） → computes inverse sin in radians → computes cosine from radian input → computes inverse cosine in radians →computes tangent from radian input（tan(3.1415/180*45)=1.0） → computes inverse tangent in radians → returns the arctangent of y/x in radians → computes hyperbolic tangent from radian input → computes e to the power of the input（exp(1)=2.718282） → computes natural log （log(10)=2.302585） → computes log base 10 （log10(10)=1.0） → computes square root of input （sqrt(2)=1.414） → returns the maximum value of an array → returns the minimum value of an array
Conversion between data types
function integertobyte function integertocharacter function integertoshort function shorttobyte function shorttocharacter function longtobyte function longtocharacter function longtoshort function longtointeger function floattobyte function floattocharacter function floattoshort function floattointeger function floattolong (int_val : integer) (int_val : integer) (int_val : integer) (short_val : short) (short_val : short) (long_val : long) (long_val : long) (long_val : long) (long_val : long) (float_val : float) (float_val : float) (float_val : float) (float_val : float) (float_val : float)
function doubletobyte function doubletocharacter function doubletoshort function doubletointeger function doubletolong function doubletofloat function charactertodouble function charactertofloat function charactertolong function charactertoshort function charactertostring function charactertointeger function stringtocharacter function stringtodouble function stringtofloat function stringtointeger function stringtolong function stringtoshort
(double_val : double) (double_val : double) (double_val : double) (double_val : double) (double_val : double) (double_val : double) (char_val : character) (char_val : character) (char_val : character) (char_val : character) (char_val : character) (char_val : character) (char_val : string) (char_val : string) (char_val : string) (char_val : string) (char_val : string) (char_val : string)
從長字串中擷取出某個片段字串 假設有一個字串 "2005120800Z" 若欲取出 “1208” 則可按照下面的方式 1208” time = "2005120800Z" word = stringtocharacter(time) stringtocharacter(time) nwod = charactertostring(word(4:7)) nwod 就會等於 1208 求餘數 val = 182 print(val%3) 2
Scalar regridding routines
f2fosh f2fsh f2gsh fo2fsh g2fsh g2gsh
from one fixed grid (including pole points) to a fixed-offset grid using spherical harmonics
假設原本資料之陣列維度為 dd（nlat,nlon），轉換後之新資料維度為 gnew（ny,nx） dd（ nlat,nlon） gnew（ ny,nx）
【Example 1】由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=181 , nx=360） 1】 由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=181 nx=360）
gnew = f2fsh（dd ,（/ny,nx/）） f2fsh（ ny,nx/
【Example 2】由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=72 , nx=144） 2】 由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=72 nx=144）
from one fixed grid to another using spherical harmonics
gnew = f2fosh（dd） f2fosh（ dd）
【Example 3】由固定網格（nlat=72, nlon=144）轉換不同解析度之固定網格（ny=73 , nx=144） 3】 由固定網格（nlat=72, nlon=144）轉換不同解析度之固定網格（ny=73 nx=144）
from one fixed grid to a Gaussian grid using spherical harmonics
from one fixed-offset grid (including pole points) to a fixed grid using spherical harmonics
gnew = fo2fsh（dd） fo2fsh（ dd）
【Example 4】由固定網格（nlat=73, nlon=144）轉換到高斯網格（T42）（ny=64 , nx=128） 4】 由固定網格（nlat=73, nlon=144）轉換到高斯網格（T42）（ny=64 nx=128）
from one Gaussian grid to a fixed grid using spherical harmonics
gnew = f2gsh（dd ,（/ny,nx/）, 42） f2gsh（ ny,nx/ 42）
【Example 5】由高斯網格（T63） (nlat=96, nlon=192) 轉換到固定網格 (ny=73 , nx=144） 5】 由高斯網格（T63） nlat=96, nlon=192) ny=73 nx=144）
from one Gaussian grid to another using spherical harmonics
gnew = g2fsh（dd ,（/ny,nx/）） g2fsh（ ny,nx/
【Example 6】由高斯網格(T63)（nlat=96, nlon=192）轉換到高斯網格(T42)（ny=64 , nx=128） 6】 由高斯網格( T63)（ nlat=96, nlon=192）轉換到高斯網格( T42)（ ny=64 nx=128）
f：包含極點。fo：不包含極點。g：gaussian grid points
gnew = g2gsh（dd ,（/ny,nx/）, 42） g2gsh（ ny,nx/ 42）
Tu＠ＣＣＵ
14/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
linint2
Interpolates from one grid to another grid using bilinear interpolation. interpolation.
Vector regridding routines
f2foshv f2fshv f2gshv fo2fshv g2fshv g2gshv
interpolates a vector pair on a fixed grid (including pole points) to a fixed-offset grid using spherical harmonics interpolates a vector pair from one fixed grid to another using spherical harmonics interpolates a vector pair on a fixed grid to a Gaussian grid using spherical harmonics interpolates a vector pair on a fixed-offset grid (including pole points) to a fixed grid using spherical harmonics interpolates a vector pair on a Gaussian grid to a fixed grid using spherical harmonics interpolates a vector pair from one Gaussian grid to another using spherical harmonics
function linint2 ( xold, yold, data, CyclicX , xnew, ynew, opt0 ) xold, yold, xnew, ynew,
if the rightmost dimension of data is cyclic. Set CyclicX to True
xi = (30,33,35,39,....., 80) [does not have to be equally spaced] spaced] yi = (0,1,2,3,....,28,29) xo = (any coordinates between 30 and 80) [inclusive] yo = (any coordinates between 0 and 29) [inclusive] fo = linint2 (xi,yi,fi, False, xo,yo, 0) (xi,yi,fi, xo,yo, sst (time,lat,lon) 1x2degree grid , @_FillValue is assigned time,lat,lon) @_FillValue lats go from 64N to 74S , lons from 1 to 359E... latnew = ispan (-74,64,2) lonnew = ispan (1,359,2) sstr = linint2_Wrap (sst&lon,sst&lat,sst, True, lonnew,latnew, 0) (sst&lon,sst&lat,sst, lonnew,latnew,
假設原本資料之陣列維度為 → u（nt,nlat,nlon）、 v（nt,nlat,nlon） nt, nlat,nlon）、 nt, nlat,nlon） 轉換後之新資料維度為 → unew（nt,ny,nx）、 vnew（nt,ny,nx） unew（ nt, ny,nx）、 vnew（ nt, ny,nx）
【Example 1】由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=181 , nx=360） 1】 由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=181 nx=360）
Compute Wind field - Global
uv2dvf (u,v,div) uv2dvg (u,v,div) uv2vrf (u,v,vor) uv2vrg (u,v,vor) uv2vrdvf (u,v,vor,div) uv2vrdvg (u,v,vor,div) uv2sfvpf (u,v, sf, vp) uv2sfvpg (u,v, sf, vp) dv2uvf (div, udiv,vdiv) dv2uvg (div, udiv,vdiv) vr2uvf (vor,urot,vrot) vr2uvg (vor,urot,vrot) vrdv2uvf(vor,div,uu,vv) vrdv2uvg(vor,div,uu,vv) sfvp2uvf(sf,vp,uu,vv) sfvp2uvg(sf,vp,uu,vv) ; u,v ==> divergence ; u,v ==> relative vorticity ; u,v ==> divergence and vorticity ; u,v ==> stream function + velocity potential ; div ==> divergent wind components ; vor ==> rotational wind components ; vor,div =>reconstruct original wind ; sf,vp ==> reconstruct original wind
f2fshv（ u , v , unew , vnew） f2fshv（ vnew）
【Example 2】由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=72 , nx=144） 2】 由固定網格（nlat=73, nlon=144）轉換不同解析度之固定網格（ny=72 nx=144）
f2foshv（ u , v , unew , vnew） f2foshv（ vnew）
【Example 3】由固定網格（nlat=72, nlon=144）轉換不同解析度之固定網格（ny=73 , nx=144） 3】 由固定網格（nlat=72, nlon=144）轉換不同解析度之固定網格（ny=73 nx=144）
fo2fshv（ u , v , unew , vnew ） fo2fshv（
【Example 4】由固定網格（nlat=73, nlon=144）轉換到高斯網格（T42）（ny=64 , nx=128） 4】 由固定網格（nlat=73, nlon=144）轉換到高斯網格（T42）（ny=64 nx=128）
f2gshv（ u , v , unew , vnew , 42） f2gshv（ 42）
【Example 5】由高斯網格（T63） (nlat=96, nlon=192) 轉換到固定網格 (ny=73 , nx=144） 5】 由高斯網格（T63） nlat=96, nlon=192) ny=73 nx=144）
g2fshv（ u , v , unew , vnew ） g2fshv（
【Example 6】由高斯網格(T63)（nlat=96, nlon=192）轉換到高斯網格(T42)（ny=64 , nx=128） 6】 由高斯網格( T63)（ nlat=96, nlon=192）轉換到高斯網格( T42)（ ny=64 nx=128）
g2gshv（ u , v , unew , vnew , 42） g2gshv（ 42）
Compute Wind field - Regional
vor = uv2vr_cfd (u,v,lat,lon, opt) ==> vor (nlat,mlon) (u,v,lat,lon, nlat,mlon) div = uv2dv_cfd (u,v,lat,lon, opt) ==> div (nlat,mlon) (u,v,lat,lon, (nlat,mlon)
opt =0 ：boundary point are set to the missing value (=u@_FillValue) (=u @_FillValue) =1 ：the u and v arrays are cyclic in longitude. The arrays longitude. should NOT include the cyclic point. The upper and lower boundaries will be set to missing (=u@_FillValue) (=u @_FillValue) =2 ：boundary points are estimated using one-sided difference oneschemes normal to the boundary. =3 ：the u and v arrays are cyclic in longitude. The arrays should NOT include the cyclic point. The upper and lower boundaries are estimated using a one-sided difference onescheme normal to the boundary. boundary.
Variable Shaping
?
can and should be done without loops
– –
use NCL syntax or functions very fast for variables in memory
Fortran : 3D 1D dd(nx,ny,nz) dd1(nn)
Tu＠ＣＣＵ
15/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
ndtooned ( val ) onedtond ( val , dims : integer )
a = (/ (/1,2,3/) , (/4,5,6/) , (/7,8,9/) /) → 3 x 3 print (ndtooned(a)) (ndtooned(a)) → (/ 1, 2, 3, 4, 5, 6, 7, 8, 9 /) a = (/1, 2, 3, 4, 5, 6, 7, 8/) print (onedtond(a,(/ 2 , 4 /))) (onedtond(a,(/ → a = ( / (/1,2,3,4/) , (/5,6,7,8/) / ) → 4 x 2
? returns a 1D array of integers or floating – beginning with start and ending with finish.
ispan ( start:integer, finish:integer, stride:integer ) start:integer, finish:integer,
1 4 7
2 5 8
3 6 9
建立一組等間隔的整數序列
time = ispan(1990,2001,2) print (time) time=ispan(1990,2001,2)*1.0 time would be Type: float
fspan ( start:float, finish:float, num:integer ) start:float, finish:float,
建立一組等間隔的實數序列 returns a 1D array of evenly spaced floating point numbers ? num is the number of points including start and finish
?
any、all、num
常用於 「if statement」 statement」
any ：Returns True if any of the values of its input evaluate as True. all ：Returns True if all the elements of the input evaluate as True. num ：Counts the number of True values in the input.
Example 1 x = new(5,float,-999) new(5,float,x = (/1.,2.,-999,4.,5./) (/1.,2.,print(any(x.ge.4)) ; should be True print(any(x.lt.0)) ; should be False print(any(x.gt.2.and.x.lt.4)) ; should be False Example 2 x = new(5,float,-999) new(5,float,x = (/1.,2.,-999,4.,5./) (/1.,2.,print(all(x.ge.5) print(all(x.ge.5) ; should be False if(all(ismissing(x))) then if(all(ismissing(x))) print("x contains all missing values, cannot print("x continue.") return end if Example 3 a = (/1,2,3,4,5/) print(num(a.gt.3)) → 2 N = num(.not.ismissing(x)) num(.not.ismissing(x))
begin b = fspan( 2, 6, 11 ) print(b) end
mask
use mask function to mask out a range of the data sets values to _FillValue that DO NOT equal mask array
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" begin ; read in netCDF file in = addfile("/wp/cgd/murphys/Data/ATMOS/atmos.nc","r") ts = in->TS(0,:,:) oro = in->ORO(0,:,:) ; mask out ocean or land data ; ocean=0, land=1, sea_ice=2 land_only = ts ocean_only = ts land_only = mask(ts,oro,1) ocean_only = mask(ts,oro,0) end
mask
range_only = sst range_only = mask( sst , (sst.ge.26) , True ) range_only = mask( sst , (sst.ge.15.and.sst.le.26) , True )
Tu＠ＣＣＵ
16/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Fortran A(nz,ny,nx)、B(nz,ny,nx) C(nz,ny,nx) T(nz,ny,nx)、P(nz,ny,nx)
do I = 1 , nn ….. enddo C = A＋B 或 C＝A－B 或 C=A*B 或 C=2A+B
sum avg stddev variance min max minind maxind
計算所有資料的總和、平均、標準差、變異數。 → 得到一個數值大小
ex : val = avg(dd) 、 val = sum(dd) avg(dd) sum(dd)
【註】如果陣列中有特殊值，可以利用下面的方式避開 設定欲忽略之特殊值 → dd@_FillValue = - 99.99 計算實際運算之點數 → N = num(.not.ismissing(dd)) num(.not.ismissing(dd))
θ(nz,ny,nx)
T(nt,nz,ny,nx)、P(nz) θ(nt,nz,ny,nx)
θ= T * ( 1000/P ) ^0.286 θ= T * (1000/conform(T,P,1))^0.286
function conform (A:numeric , B:numeric , ndim:integer)
【例】dd = conform (x(time|:,lat|:,lon|:,lev|:), dp, (/0,3/)) → x(time,lev,lat,lon)、dp(ntim,klvl)、dd(ntim,nlat,mlon,klvl) 【例】dd = dim_sum(q*conform(q,dz,3)) → q(nt,ny,nx,nz)、dz(nz) → dd(nt,ny,nx)
x = (/3., 2., 5., 1., 5., 2., 5., 1., 3., 2./)
i = maxind(x) maxind(x) j = minind(x) minind(x)
不同維度陣列的對應
let SST be (100,72,144) and SICE = -1.8 (scalar) SST = SST > SICE [f90: where (sst.lt.sice) sst = sice]
Compute the rightmost dimension index
dim_avg dim_sum dim_stddev dim_min dim_max dim_median dim_num dim_rmsd dim_rmvmean dim_rmvmed dim_variance dim_stat4
dim_standardize
compute the rightmost (n-1 th) dimension at all other indices (n- th) 計算最後一個維度的平均、總和、標準差、去除平均值…等。 計算最後一個維度的平均、總和、標準差、去除平均值… 【例】dd1 = new((/nt,nm/),"float") → 定義陣列大小 new((/nt,nm/),"float") dd3 = new((/3,nt/),"float") → 定義陣列大小 dd1!0 = "yy" → 定義第一個維度的名稱為 “yy” "yy" yy” dd1!1 = "mm" → 定義第二個維度的名稱為 “mm” mm” dd3(0,:) = dim_avg(dd1(yy|:,mm|5:7)) → 計算JJA的平均 計算JJA的平均 dd3(1,:) = dim_sum(dd1(yy|:,mm|5:7))→ 計算JJA的總和 dim_sum(dd1(yy|:,mm|5:7))→ 計算JJA的總和 dd3(2,:) = dim_rmvmean(dd3(0,:)) → 去掉平均值
Compute n dimension index
dim_avg_n dim_sum_n dim_stddev_n dim_min_n dim_max_n dim_median_n dim_num_n dim_rmsd_n dim_rmvmean_n dim_rmvmed_n dim_variance_n dim_stat4_n
計算某個維度的平均、總和、標準差、去除平均值…等。 計算某個維度的平均、總和、標準差、去除平均值…
【例】Let q be of size ( nt, ny, nx) and with named nt, ny, nx) dimensions "time", "latitude" and "longitude",
qAvgTime = dim_avg_n( q, 0) )) dim_avg_n( qAvgLon = dim_avg_n( q, 2 ) dim_avg_n(
【例】 Let q be defined q( time , lev , lat , lon ). qSum = dim_sum_n(x, (/ 0 , 1 /)) dim_sum_n(x,
計算平均時，如果採用一般的點數作計算，則裡面所使用的維 度範圍需用整數，如果已用實際經緯度作對應則無此限制。
dim_standardize_n qSum = dim_sum_n(x,(/ 2 , 3 /)) dim_sum_n(x,(/
Available in version 5.1.1 or later.
wgt_runave wgt_areaave wgt_areaave2 wgt_areasum2 wgt_arearmse wgt_arearmse2 wgt_volave wgt_volrmse
wgt_runave(dat, wgt [*], opt) → pass filter wgt_areaave(dat, wgty[*], wgtx[*], opt) wgt_areaave2(dat, wgt [*][*], opt) wgt_areasum2(dat, wgt [*][*], opt) wgt_arearmse(dat1,dat2,wgty[*],wgtx[*],opt) wgt_arearmse2(dat1,dat2, wgt[*][*],opt) wgt_volave(dat1,dat2,wgtz[*],wgty[*],wgtx[*],opt) wgt_volrmse(dat1,dat2,wgtz[*],wgty[*],wgtx[*],opt)
opt = 0, the area average is calculated using available non-missing data. opt = 1, then if any point in dat is missing, the area average is not computed.
; ########################################## rlon = lonGlobeF(nx, "lon", "longitude", "degrees_east") rlat = latGau (ny, "lat", "latitude", "degrees_north") wgty = latGauWgt(ny, "lat", "gaussian weights", "dimension_less") ; ###########################################
dd = nameDim (dd, (/"time","lat","lon"/), “Temp", “K")
dd&lat = rlat dd&lon = rlon wgty!0 = "lat" wgty&lat = rlat ; ############################################# gb = wgt_areaave(dd( : ,{-90:90}, : ) , wgty({-90:90}), 1.0, 0) nh = wgt_areaave(dd( : ,{ 0:90}, : ) , wgty({ 0:90}), 1.0, 0) sh = wgt_areaave(dd( : ,{-90: 0}, : ) , wgty({-90: 0}), 1.0, 0) ; #############################################
zonalAve
zonalAve
Computes a zonal average of the input array.
zonalAve ( x : numeric )
reture_val : typeof (x)
Tu＠ＣＣＵ
17/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Compute gradient、detrend、Laplacian
angmom_atm
計算大氣的相對角動量（relative angular momentum，單位=kg * m^2/s）
angmom_atm(u, dp,lat[*],wgt[*])
( containing gaussian weights if u is on a gaussian grid. Otherwise, set to a scalar (typically 1.0). )
ram_NH = angmom_atm (u(:,:,{0:90},:), dp, lat({0:90}), 1.0) ; ==> ram_NH(time)
gradsf gradsg
dtrend
計算球面座標下x和y方向的梯度（gradient）
gradsf (dat , grad_x , grad_y) dtrend (y:numeric , return_info: logical)
fluxEddy
fluxEddy(x,y) → x‘y’ = ave{x*y} - X*Y（其中x=X+x‘ and y=Y+y’） The rightmost dimension must be "time". Missing values should be indicated by x@_FillValue 【例】upvp = fluxEddy (u(lev|:,lat|:,lon|:,time|:), v(lev|:,lat|:,lon|:,time|:)) upvp will be a 3-dimensional array dimensioned klev x nlat x mlon.
Calculates the pressure of the lifting condensation level（LCL）.
dtrend_msg(x[*],y,remove_mean:logical,return_info: logical) dtrend_msg 【例】yDtrend = dtrend (y(lat|:,lon|:,time|:), False) dtrend_quadratic
lclvl
【例】yDtrend = dtrend_msg (y&time, y(lat|:,lon|:,time|:), True, False)
lclvl(P , T , Td)（units=hPa and K）
【例】plcl = lclvl (1000, 15+273.15, 4+273.15) ; plcl = 848.6 【例】plcl = lclvl (P, T, Td) ; P＆T＆Td(nt,ny,nx)、plcl(nt,ny,nx)
hydro
computing the geopotential height using the hydrostatic equation
hydro( P: numeric, Tv : numeric, Zsfc : numeric) P：pressure (hPa) (last dimension is nlvl and Pressure must start at the surface) Tv：virtual temperature (K) at each P（Tv= (T+273.15)*(1.+q*0.61)） Zsfc：surface geopotential height (gpm) (same dimensions as P)
【例】Tnew = T(time|:,lat|:,lon|:,lev|:)zh=hydro(conform(Tnew,P,3),Tnew, Zsfc) （ntim,nlat,mlon,klev）
prcwater_dp
computing column precipitable water（unit= [kg/m2]）
prcwater_dp(q : numeric, dp : numeric)
q ：specific humidity (kg/kg). The rightmost dimension must be “level". dp ：pressure layer thickness (Pascals). 【例】q = (/ 20.8,19.4,16.5, ... , 1.7 ,1.0 /) ; [g/kg] dp = (/ 10., 30. , 50., ... , 100., 50. /); [hPa] pw = prcwater_dp(q*0.001,dp*100.) ; [kg/m2]
lapsF lapsG lapsf lapsg lapvf lapvg ilapsF ilapsG
Given a scalar z, compute its Laplacian via Spherepack.
lapsF( z : numeric ) lapsG( z : numeric ) lapsf ( z : numeric, zlap : float ; or double ) lapsg ( z : numeric, zlap : float ; or double )
Given vector quantity (u,v), compute the vector Laplacian via Spherepack
lapvf ( u : numeric, v : numeric, ulap : float, vlap : float ) lapvg ( u : numeric, v : numeric, ulap : float, vlap : float) z ：scalar function array (input, two or more dimensions, last two dimensions must be nlat x nlon and input values must be in ascending latitude order). zlap：Laplacian array (output, same dimensions as z, values will be in ascending latitude order)
Running average
function runave( dat : numeric, nave : integer, kopt : integer) runave(
Model
dat ：1或N維資料陣列。可利用 「dat@_FillValue」設定陣列之 dat@_FillValue」設定陣列之
特殊值。若沒設定，於計算過程中將被當作真實值計算。
nave ：Number of points to be included in the running average
(一般使用奇數) 一般使用奇數)
kopt ：邊界條件
Dtrend
kopt < 0 : utilize cyclic conditions kopt = 0 : set unsmoothed beginning and end pts to x@_FillValue kopt > 0 : utilize reflective (symmetric) conditions dd3(3,:) = runave(dd3(0,:),12,0) → 計算12個月的滑動平均 計算12個月的滑動平均
linear interpolation and smoothing
int2p
int2p(prsin , dat , prsout , linlog[1] : integer)
linear interpolation (=1) or linear interpolation in ln (≠1) 【Example】
plev@units = "hPa" plev = (/1000,925,850,700,600,500,400,300,250,200,150,100/) plevn = ispan(100,1000,50) ddt = int2p(plev,tmp,plevn(::-1),2)
vertical average： v2p = dim_avg(int2p(plev({lev1:lev2}),dat(:,:,{lev1:lev2}),plevn(::-1),linlog))
smth9
1------8------7 | | | 2------0------6 | | | 3------4------5
performing nine point local smoothing on a 2D grid
smth9 (x[*][*] : numeric , p : numeric , q : numeric , wrap : logical)
f0 = f0 + (p/4)*(f2+f4+f6+f8-4*f0) + (q/4)*(f1+f3+f5+f7-4*f0)
x ：二維或以上之網格資料陣列。平滑只取最右邊兩個維度 若有特殊值需先作x@_FillValue的設定。 p,q：權重。（通常p設為0.5，q則設為0.25） 【註】 q = -0.25 → results in "light" smoothing q = 0.25 → results in "heavy" smoothing. q = 0.0 → results in a 5-point local smoother. wrap :True when the rightmost dimension is treated as cyclic. False if the rightmost dimension is not to be treated as a cyclic. 【Example】xsmth = smth9 (x , 0.50, 0.25, True) ; heavy local smooth
Tu＠ＣＣＵ
18/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
triple2grid2d
將非網格點資料內差成二維網格點資料
triple2grid ( x [*], y [*], data [*], xgrid [*], ygrid [*], option) triple2grid2d ( x [*], y [*], data [*], xgrid [*][*], ygrid [*][*], option)
option=False, the function will operate under default mode option=False, that is, all grid points will be assigned a value. option= True, this variable may have associated with it option= the attributes mopt and/or distmx. distmx.
grid
= triple2grid2d(xlon,ylat,zVal, lon,lat, False) lon,lat,
opt = True ; setting options opt@mopt = 1 grid = triple2grid2d(xlon,ylat,zVal, lon,lat, opt) lon,lat,
option@mopt - an integer value of 0 or 1
option@mopt=0 (default) option@ mopt=0 - use a quick approximation to determine the distance of an x/y location to a grid point. option@mopt=1 option@ mopt=1 - calculate the distance using the great circle distance formula. ( in slower execution times.)
Setting the option@distmx option will result in any option@ observation greater than distmx to be not used: opt@distmx = 1.2 ; observation more than 1.2 will not be used grid = triple2grid(xlon,ylat,zVal, lon,lat, opt) triple2grid(xlon,ylat,zVal, lon,lat,
option@distmx
Any observation greater than option@distmx from a grid point will not be used. option@ If this option is used, it is possible for some grid points to be returned as missing be
opt = True opt@mopt = 0 opt@distmx = 10.0 res@sfXArray = rlon hov = triple2grid(rlon,rlat,rain,rlon1,rlat1,opt) res@sfYArray = rlat plot = gsn_csm_contour_map(wks,rain,res) hov = smth9(hov,0.5,0.25,False) gsn_csm_contour_map(wks,rain,res) plot = gsn_csm_contour_map(wks,rain,res) gsn_csm_contour_map(wks,rain,res)
Compute climatology and anomaly
計算氣候平均值，時間需為12的倍數。（計算速度慢，所需計算時間大約是 dim_avg的三倍） mean = clmMonLLT (x[*][*][*]:numeric) (lat,lon,time) mean = clmMonLLLT (x[*][*][*][*]:numeric) (lev,lat,lon,time) mean = clmMonTLL (x[*][*][*]:numeric) (time,lat,lon) mean = clmMonTLLL (x[*][*][*][*]:numeric) (time,lev,lat,lon) 【註】執行前需先設定陣列中每個維度的名稱，例如： x!0 = "time" x!1 = "lat" x!2 = "lon" 計算距平值 anom = calcMonAnomLLT(x[*][*][*], xAve[*][*][12]) anom = calcMonAnomLLLT(x[*][*][*][*],xAve[*][*][12]) anom = calcMonAnomTLL(x[*][*][*], xAve[12][*][*]) 計算均方根 x = stdMonLLT (x[*][*][*]:numeric) x = stdMonTLL (x[*][*][*]:numeric) x = stdMonLLLT (x[*][*][*][*]:numeric) x = stdMonTLLL (x[*][*][*][*]:numeric) ; x(lat,lon,time) ; x(time,lat,lon) ; x(lev,lat,lon,time) ; x(time,lev,lat,lon)
clmMon
calcMonAnom
stdMon
rmMonAnnCyc
Removes thae Annual Cycle from monthly (nmos=12) data. x = rmMonAnnCycLLT(x) ; input dim. x(lat,lon,time) x = rmMonAnnCycTLL(x) ; input dim. x(time,lat,lon) x = rmMonAnnCycLLLT(x) ; input dim. x(lev,lat,lon,time)
【結果與下面計算式相同】
xAve = clmMonLLT(x) xAnom = calcMonAnomLLT(x,xAve)
Compute seasonal mean
計算季節平均
[ DJF , JFM , FMA , MAM , AMJ , MJJ , JJA , JAS , ASO , SON , OND , NDJ ] 執行前需先設定陣列中每個維度的名稱，資料比數必須被12整除。資料必須從一月 開始擺放，計算時頭尾兩筆資料實際上是兩個月的平均（DJF=JF ， NDJ=ND）
Vertical integration
Function for doing vertical integration in pressure using beta factors
month_to_season
usage result = month_to_season (xMon:numeric, SEASON:string) example for : xMon(time,lat,lon) , result : xJJA(time/12,lat,lon) xJJA = month_to_season (xMon, "JJA")
function vibeta( p : numeric , x: numeric , linlog : integer , psfc: numeric , pbot: numeric , ptop: numeric)
vibeta
p ：pressure levels（氣壓值從低層（bottom）往上層（top）排） X ：dat（the rightmost dimension of dat must be same length as p） linlog：內差方式。「= 1」表使用線性內差；「= 2」表使用log內差。 Psfc ：地面氣壓（必須與dat相同維度） pbot, ptop：欲積分範圍。pbot表最底層，ptop表最上層。 【Example】：Let dat(nt,nz,ny,nx) , psfc(nt,ny,nx) , p(nz) → vint(nt,ny,nx) vint = vibeta (p ,dat(time|:,lat|:,lon|:,lev|:) , linlog , psfc , pbot , ptop)
計算n個季節平均值
[ DJF , JFM , FMA , MAM , AMJ , MJJ , JJA , JAS , ASO , SON , OND , NDJ ]
month_to_seasonN
usage result = month_to_seasonN (xMon:numeric, SEASONS[*]:string) example for : xMon(time,lat,lon) , result : xJJA(n,time/12,lat,lon) xJJA = month_to_seasonN (xMon, (/"JJA","JFM"/))
ftcurvi simpeq simpne
integrates a sequence datasets using cubic splines cub = ftcurvi(lev2,lev1,plev({lev2:lev1}),dat(:,:,{lev2:lev1})) Integrate a sequence of equally spaced points using Simpson's Rule
function simpseq ( dat: numeric , dx[1] : numeric )
result = simpson ( f (lev|:,lat|:,lon|:,time|:) , dt ) ; result(nz,ny,nx)，dt = 1.0（constant interval between points） Integrate a sequence of unequally spaced points using Simpson's three-point formula.
計算季節的氣候平均值
[ DJF , JFM , FMA , MAM , AMJ , MJJ , JJA , JAS , ASO , SON , OND , NDJ ]
month_to_season12
usage result = month_to_season12 (xMon:numeric) example for : xMon(time,lat,lon) , result : xJJA(12,lat,lon) xJJA = month_to_season (xMon)
function simpne(x: numeric,y: numeric)
【Example】Let p(nz) and f(nt,nz,ny,nx) → result(nt,ny,nx) result = simpne(p,f(time|:,lat|:,lon|:,lev|:))
Tu＠ＣＣＵ
19/41
07/23/2009
中國文化大學大氣科學系暑期電腦訓練課程：ＮＣＬ
Statistic method：correlations
relhum ( Ta: numeric , w: numeric , prs: numeric ) prs:
（給定溫度、混合比、氣壓，計算相對濕度（單位為％））
mixhum_ptrh ( prs: numeric , Ta: numeric , rh: numeric , opt: integer ) prs: rh:
（給定氣壓、溫度以及相對濕度，計算比濕或混合比）
ccr = esacr(dat1,maxlag) ccr = esccr(dat1,dat2,maxlag) ccr = esacv(dat1,maxlag) ccr = esccv(dat1,dat2,maxlag) ccr = escorc(dat1,dat2) ccr = escovc(dat1,dat2)
; auto correlations ; cross correlations（ data1 lead data2） ; auto covariances ; cross covariances ; between two variables at zero lag only
mixhum_ptd( prs: numeric , Td: numeric ,opt : integer ) mixhum_ptd( prs:
（給定氣壓和露點溫度，計算比濕或混合比）
prs Ta Td rh w opt
：pressure (hPa/mb) (any dimensionality) (hPa/mb) ：temperature (K) at each p (same dimensionality as p) ：dew point temperature (K) at each p ：relative humidity (%) (same dimensions as p) : Mixing ratio (kg/kg) ：opt = 1, mixing ratio (kg/kg). 。 ＝ -1（ g/kg） g/kg） opt = 2, specific humidity (kg/kg)。 ＝ -2（ g/kg） (kg/kg)。 g/kg） prs=(/ 1008.,1000.,950.,900.,850.,800.,750.,700.,650.,600/) ; hPA prs=(/ ta=(/ 29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1 /) ; ℃ ta=(/ /) rh=(/ 75.0,60.0, 61.1, 76.7, 90.5, 89.8, 78.3, 76.5, 46.0,55.0 /) ; % rh=(/ w=(/20.38, 19.03,16.14,13.71,11.56,9.80,8.33,6.75,6.06,5.07/) ;g/Kg
【註】假設 q(t) and s(t+k) ，"k" 從0 到 mxlag auto correlations:c(k) = SUM [(q(t)-qAve)*(q(t+k)-qAve)}]/qVar cros correlations:c(k) = SUM [(q(t)-qAve)*(s(t+k)-sAve)}]/(qStd*sStd) auto covariancesc(k) = SUM [(q(t)-qAve)*(q(t+k)-qAve)}]/N cross covariancesc(k) = SUM [(q(t)-qAve)*(s(t+k)-sAve)}]/N cov = SUM [(X(t)-Xave)*(Y(t)-Yave)}]/(NT-1) cor = SUM [(X(t)-Xave)*(Y(t)-Yave)}]/(Xstd*Ystd)
q = mixhum_ptrh (prs , ta+273.15 , rh*0.01 , 1 ) ; mix ratio (kg/kg) rh*0.01 rh = relhum(ta+273.15 , w*0.001 , prs*100.) prs*100.)
Statistic method：correlations
Auto & cross correlation
x(nt) , y(nt) x(nt) y(nt) x(nt) , y(ny,nx,nt) x(nt) y(ny,nx,nt) x(nz,nt) , y(ny,nx,nt) x(nz,nt) y(ny,nx,nt) x(nb,na,nt) , y(nz,ny,nx,nt) x(nb,na,nt) y(nz,ny,nx,nt) x(ny,nx,nt) , y(ny,nx,nt) x(ny,nx,nt) y(ny,nx,nt) → c(mxlag) c(mxlag) → c(ny,nx,mxlag) c(ny,nx,mxlag) → c(nz,ny,nx,mxlag) c(nz,ny,nx,mxlag) → c(nb,na,nz,ny,nx,mxlag) c(nb,na,nz,ny,nx,mxlag) → c(ny,nx,mxlag) c(ny,nx,mxlag)
;************************************************ ; calculate cross correlations ; note, the max lag should not be more than N/4 ;************************************************ x_Lead_y = esccr(x,y,mxlag) maxlag = 25 ccr = esccr(T1,T2,maxlag) (T1,T2,maxlag) xcd = ispan(0,maxlag-1,1) ispan(0,maxlag;************************************************ wks = gsn_open_wks("ps","corel") gsn_open_wks("ps","corel") res = True ; make plot mods res@tiMainString = "37.7N 180E vs 23.72S 149W" res@tiXAxisString = "LAG" ; x-axis label x;************************************************ plot = gsn_xy(wks,xcd,ccr,res) ; plot correlation gsn_xy(wks,xcd,ccr,res)
between two variables
x(nt) , y(nt) x(nt) , y(ny,nx,nt) x(nz,nt) , y(ny,nx,nt) x(nb,na,nt) , y(nz,ny,nx,nt) x(ny,nx,nt) , y(ny,nx,nt) → val (a scalar) → c(ny,nx) → c(nz,ny,nx) → c(nb,na,nz,ny,nx) → c(ny,nx)
ccr2
ccr1
;************************************************ maxlag = 25 totlag = maxlag * 2 ;************************************************ ccr1 = esccr(T1,T2,maxlag) ( ,maxlag) ccr2 = esccr(T2,T1,maxlag) ( ,maxlag) ;************************************************ ccrs = new( (/totlag/),float) (/totlag/),float) ccrs(maxlag:) = ccr1(1:maxlag-1) ccrs(maxlag:) ccr1(1:maxlagccrs(0:maxlag-1) = ccr2(0:maxlag-1:-1) ccrs(0:maxlagccr2(0:maxlag- 1:xcd = ispan(-maxlag+1,maxlag-1,1) ispan(- maxlag+1,maxlag;************************************************ plot = gsn_xy(wks,xcd,ccrs,res) gsn_xy(wks,xcd,ccrs,res)
Statistic method：Regression (1/3)
rc = regline( x: numeric , y: numeric )
Calculates the linear regression coefficient between two series
;************************************************ T1 = new(/ny,nx,nt/) T2 = new(/ny,nx,nt/) ccr = new(/ny,nx,maxlag/) ;************************************************ maxlag = 2 ccr = esccr(T1,T2,maxlag) lag = 0 res@tiMainString = "Correlations at lag "+lag plot = gsn_csm_contour_map_ce(wks,ccr(:,:,lag),res)
regline also returns the following attributes:
xave ：average of x (scalar, float or double) yave ：average of y (scalar, float or double) tva ：t-statistic (assuming null-hypothesis) nullnptxy ：number of points used yintercept：y intercept
Tu＠ＣＣＵ
20/41
07/23/2009
图文详解Windows平台上NCL的安装 NCL在Linux下的安装非常容易，只需下载适当版本的文件，设置好环境变量即可使用。 NCL在Windows下的安装则要麻烦一些，需要先安装一个虚拟Linux环境（Cygwin/X）。 以下内容详细介绍NCL在Windows平台上的安装过程，希望仅具备Windows基本操作技能的用户也能轻松安装NCL。 一、NCL简介 二、准备工作 三、安装Cygwin/X 四、熟悉Cygwin/X环境 五、安装NCL 六、运行NCL范例 七、语法高亮显示（此部分供有兴趣的用户参考） 八、.hluresfile文件（此部分供有兴趣的用户参考） 九、FAQ 十、获取帮助 一、NCL简介 NCL（NCAR Command Language）是由NCAR的“Computational & Information Systems Laboratory”开发的。 NCL是一种编程语言，专门用于分析和可视化数据。主要用于以下三个领域： 文件输入/输出（File input and output）： 资料处理（Data processing）： 图形显示（Graphical display）：可生出出版级别的黑白、灰度或彩色图。 从5.0起，NCL和NCAR Graphics已经打包在一起发行。2009年3月4日，NCL发布了最新的5.1.0版，该版本更新了地图投影，修正了一些bug，增加了更多的函数及资源。下图为新增的含中国省界的地图（见图1-1）。
二、准备工作 2.1 安装环境 安装环境为WinXP Professional SP3，并做如下假定： 计算机名：TEAM 用户名：Grissom 安装目录：D:\download 用户在实际安装中，请根据自己系统的信息替换本教程中的计算机名和用户名。 特别说明：用户名中不能出现空格，否则会在使用中出现一些问题。 2.2 下载Cygwin/X Cygwin/X=Cygwin+X。通俗地说，Cygwin/X可以在Windows平台上实现命令行+图形的Linux模拟环境。 Cygwin/X的下载与安装非常灵活，用户可根据自己的需求定制。为便于大家的安装，我已下载了安装NCL所需的软件包，包括编译器、编辑器、X Server等，用户可直接从以下地址下载，并解压至D:\download\install 目录下。 Cygwin下载：https://www.doczj.com/doc/c112425665.html,/xglm/2009/2/wnx45afnq7.htm 以下关于Cygwin和Cygwin/X的详细介绍供参考：
气象项目NCL开发环境配置手册 1建立和测试源码安装ncl所需要的编译环境 1.1安装ncl所需要的编译器 源码编译安装ncl需要C编译器和Fortran编译器。C编译器使用gcc即可。安装ncl5.2.1版本时，最好使用gfortran或g95作为Fortran编译器，而不要使用g77。 1.2为外部软件配置环境变量 例如设置C编译器环境变量，export CC=gcc；设置Fortran编译器，export FC=gfortran。 2下载安装非可选外部软件 注：安装时最好把所有的外部软件都安装到同一根目录下，这样做便于以后告诉ncl编译系统所有的外部软件的安装位置。本文假设所有的外部软件都安装在/usr/local目录下。官网上说如果源码编译安装ncl，则下面的几款软件都是必须安装的： ●JPEG 支持jpeg图形的软件，我下载的jpeg源码安装文件是jpegsrc.v8c.tar.gz。一旦有了源码，执行以下命令进行安装： ./configure --prefix=/usr/local make all install 如果jpeg版本是v6的，则还要额外执行以下命令： make install-lib make install-headers ●zlib 如果想要支持png图形，或者支持grib2数据，则需要下载安装此软件。我下载的zlib源码安装文件是zlib-1.2.5.tar.gz。一旦有了源码，执行以下命令进行安装： ./configure - -prefix=/usr/local make all install ●NetCDF 支持NetCDF数据格式读取的软件包。（如果不需要NetCDF数据读取的话，应该可以 不用安装） - 1-
Chapter 1: Overview Introduction The Advanced Research WRF (ARW) modeling system has been in development for the past few years. The current release is Version 3, available since April 2008. The ARW is designed to be a flexible, state-of-the-art atmospheric simulation system that is portable and efficient on available parallel computing platforms. The ARW is suitable for use in a broad range of applications across scales ranging from meters to thousands of kilometers, including: ?Idealized simulations (e.g. LES, convection, baroclinic waves) ?Parameterization research ?Data assimilation research ?Forecast research ?Real-time NWP ?Coupled-model applications ?Teaching 简介 Advanced Research WRF (ARW)模式系统在过去的数年中得到了发展。最近公布了第三版，从2008年4月开始可供使用。ARW是灵活的，最先进的大气模拟系统，它易移植，并且有效的应用于各种操作系统。ARW适用于从米到成千上万公里尺度的各种天气系统的模拟，它的功能包括： ?理想化模拟（如，LES，对流，斜压波） ?参数化研究 ?数据同化研究 ?预报研究 ?实时数值天气预报 ?耦合模式应用 ?教学 The Mesoscale and Microscale Meteorology Division of NCAR is currently maintaining and supporting a subset of the overall WRF code (Version 3) that includes: ?WRF Software Framework (WSF) ?Advanced Research WRF (ARW) dynamic solver, including one-way, two-way nesting and moving nest. ?The WRF Preprocessing System (WPS) ?WRF Variational Data Assimilation (WRF-Var) system which currently supports 3DVAR capability ?Numerous physics packages contributed by WRF partners and the research community ?Several graphics programs and conversion programs for other graphics tools And these are the subjects of this document. The WRF modeling system software is in the public domain and is freely available for community use.
叉车电瓶和充电器使用规程 一、电瓶： 当叉车配备了其它种类的电瓶及/或充电器时，遵从其制造商的指导。 （一）远离火种（爆炸性气体） 绝对禁止火焰接近电瓶。电瓶内部会产生爆炸性气体；吸烟、火焰及火花，均会引起电瓶爆炸。 （二）小心触电 严禁引起短路，电瓶带有高电压和能量。当处理电瓶时，要戴护目镜、穿胶鞋和戴橡皮手套。 （三）正确接触 严禁将电瓶的正、负极调乱，否则可导致火花、燃烧及/或爆炸。 （四）远离工具 严禁让工具接近电瓶两极，以免引起火花或短路。 （五）勿过量放电 严禁让叉车的电量耗至叉车不能移动时，才进行充电。（会引致电池寿命缩短）当电瓶负荷显示器显示无电时，请立刻进行充电。 （六）保持清洁 保持电瓶上表面干净。严禁使用干布擦电瓶表面，以免引起静电。清洁电瓶要戴护目镜、穿胶鞋和戴橡皮手套。当清洁电瓶后，才能进行充电。 （七）穿着安全服 为了个人安全，需配戴护目镜、穿胶鞋和戴橡皮手套。 （八）小心电瓶电解液。 电瓶电解液含有硫酸，严禁让皮肤接触电瓶电解液。 （九）急救 电瓶电解液含有硫酸，若与之接触可能造成烧伤。发生意外时，请立即进行急救并请医生治理。 ①溅到皮肤上：用水冲洗10-15分钟 ②溅到眼睛上：用水冲洗10-15分钟 ③误咽：喝大量的水和牛奶 ④溅到衣服上：立即脱下衣服 如不遵守以上各点，可造成严重的伤亡事故。 （十）拧紧电瓶通风盖 拧紧电瓶通风盖，以防泄漏电解液。 （十一）清洗 严禁在叉车上清洗电瓶，以免损坏叉车。 （十二）不正常的电瓶 如电瓶发生下列情况，与和资格的服务机构或电瓶生产商联系。 ①电瓶发臭 ②电解液变浊 ③电解液减少速度过快 ④电解液温度过高 （十三）严禁拆解电瓶 不得让电解液耗尽、拆卸或自行维修电瓶。
WRF-Chem安装及使用说明手册 V1.0 编写者：王彬 2015年6月
目录 1. 引言 (3) 1.1 编写目的 (3) 1.2 背景 (3) 1.3 参考资料 (4) 2. 安装过程 (4) 2.1 设定环境变量 (4) 2.2 WRFV3安装 (5) 2.3 WPS与WRFDA的安装 (6) 3. 使用过程 (6) 3.1 WPS数据准备说明 (7) 3.2 气象初始场的准备说明 (8) 3.3 化学数据前处理程序的使用说明 (8) 3.4 化学数据转化程序的使用说明 (9) 附录1 （prep_chem_sources.inp） (14) 附录2（ncl程序） (18) 附录3（convert.f90） (19)
1. 引言 1.1 编写目的 本手册为指导WRF-Chem数值模拟的用户，方便轻松掌握WRF-Chem数值模式的安装及使用而编写。希望该手册在大家使用和学习过程中起到引导和帮助的作用。 1.2 背景 WRF-Chem模式是由国家大气研究中心（NCAR）等机构联合开发的新一代大气预报模式，真正实现了气象模式与化学传输模式在时空上的耦合。它不仅能够模拟污染气体和气溶胶的排放、传输以及混合过程，还可用于分析空气质量、云与化学之间的相互作用等。研究表明，WRF-Chem对气象场以及污染物的模拟表现出值得信赖的能力。 WRF-Chem在原有WRF模式的基础上耦合了化学反应过程，它通过在WRF 的辐射方案中引入气溶胶光学厚度、单次散射反照率和不对称因子表现气溶胶的直接辐射效应。此外，还增加了云滴数浓度的计算，这会改变原有WRF辐射方案中云粒子有效半径的计算方式，从而影响云的反照率，体现气溶胶的Twomey 效应。另一方面，Liu et al.提出了以云滴数浓度为阈值的新的云水向雨水自动转换参数化方案，取代了原有的Kessler方案，体现气溶胶的Albrecht效应。
小型机UNIX系统操作指南 一、 UNIX基本指令 1、Ls –l/-a/-R 列目录 2、Cd 更改目录 3、ps –ef |grep ora_ 看进程 4、kill -9 pid 杀进程 5、more filename 按屏看文本 6、cat filename 看文本 7、strings filename 看文本 8、df 看文件系统 9、du –k 看空间 10、rm 删文件 11、rm –r dir 删目录 12、find . –name filename 查找文件 13、cp *.* /aaa 拷贝文件 14、mv 移动文件 15、tail –f file1 看文件尾，可临控文件变化 16、makedir 建目录 17、id 看用户 18、mount /cdrom 19、>file1 20、set –o vi 可保存历史命令 21、su – user 切换用户 22、passwd 改密码 23、man 查帮助 24、ifconfig 看网卡状态 25、date 看时间
26、chmod 改权限 27、chown 改组/用户 28、vi 编辑器，常用的操作 二、 AIX （一）、文件系统 1、VG、PV、LV、LP、PP的概念及关系 2、基本文件系统 root(/) :home用户、usr软件、tmp临时、var控制管理文件（二）、常用指令 1、lsdev –C 看设备 2、info 相关书籍 3、lslpp –l 安装的软件 4、clstat 监控群集状态 5、errpt –a|more 错误日志 6、errclear 0 删除错误日志 7、df –k 看文件系统空间 8、ps –aux 列出当前进程 9、lspv –o 看激活的VG 10、smitty mkgroup\group\rmgroup 对组操作 hacmp 对群集操作 user\mkuser\rmuser 对用户操作 chinet 改IP tcpip 网络配置 11、varyonvg sharedvg 挂VG varyoffvg sharedvgb 卸VG 12、diag 故障诊断 13、cfgmgr –v 重新配置硬件 14、netstat –in 看网卡信息 15、lsvg –p 看VG （二）、系统开、关机
用NCL处理卫星数据(一个简单的教程) 2010年03月19日星期五15:13 这个教程主要介绍如何用NCL读取卫星数据，如何将数据存储为二进制格式，如何作图。 本教程所用的数据：https://www.doczj.com/doc/c112425665.html,/file/f19ba404fc 这个数据是AMSR-E的L2A产品，更多信息请查阅相关文档。 在开始这个教程前，请确认已经安装好了NCL。 1. 运行NCL 和GrADS类似，NCL也可以运行在交互式环境中，即输入一条命令执行一条命令。 交互式环境使用： 在命令行下输入 $ ncl 退出请输入命令 ncl 0> exit 交互式环境不方便重复操作，一旦退出程序就需要全部重新来过，所以推荐编写NCL脚本。一个NCL脚本主要包括以下几个部分： load 外部函数 begin ;命令和操作 end load用来载入外部函数。 文件操作、数据操作、数据可视化等代码都放在load 语句之后。 半角分号“;”后的内容是注释，NCL只有行注释功能。 begin 和end 可有可无，但如果出现了begin，就必须有end 与之对应。 NCL脚本扩展名为.ncl。比如有一NCL脚本文件为script.ncl，运行时可以用下面的方法：$ ncl < script.ncl 或$ ncl script.ncl
下面来写一个简单的NCL脚本。打开vi之类的编辑软件，将下面的内容复制进去，然后保存到卫星数据文件所 ;******************************************************** ; Load NCL scripts 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl" ; load语句载入了NCL常用的几个外部程序，建议以后自己写脚本时也保留这6行的内容。 ;******************************************************** begin fi = addfile("./AMSR_E_L2A.hdf","r") print(fi) end 运行这个脚本 $ ncl test.ncl > log 2. 熟悉数据 打开log文件 从第11行相关信息可以知道该文件为HDFEOS_V2.5 第663行开始可以得知该卫星数据的生成时间等信息 第957行开始是卫星数据的维度信息(dimensions) 第983行起是变量信息(variables) 这些对文件内容进行描述的信息被称作元数据(metadata)，元数据对于数据的存储和发布很重要。它记录了数据以及数据集开发的历史等。 3. 读取数据 下面以89V的数据为例 89V的元数据为： short89_0V_Res_1_TB__4 ( DataTrack_lo_Low_Res_Swath, DataXTrack_lo_Low_Res_Swath ) SCALE_FACTOR : 0.01 OFFSET : 327.68 UNIT : Kelvin hdf_name : 89.0V_Res.1_TB hdf_group : Low_Res_Swath/Data Fields hdf_group_id : 4