Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
9afdcdd
commit 2d39258
Showing
31 changed files
with
2,190 additions
and
2 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,109 @@ | ||
undef("PrnOscPat_driver") | ||
function PrnOscPat_driver(eof[*][*][*]:numeric, eof_ts[*][*]:numeric, kPOP[1]:integer) | ||
; ================================================================= | ||
; compute Principal Oscillation Patterns (POPs) | ||
; ================================================================= | ||
local dim_ts, dim_eof, neof, ntim, nlat, mlon, dnam_ts, dnam_eof, neof, j \ | ||
, cov0, cov1, cov0_inverse, A, z, Z, pr, pi, zr, zi, mean, stdev \ | ||
, evlr, eigi, eigr | ||
begin | ||
|
||
dim_ts = dimsizes(eof_ts) ; (neof,ntim) | ||
dim_eof = dimsizes(eof) ; (neof,nlat,mlon) | ||
|
||
ntim = dim_ts(1) | ||
neof = dim_eof(0) | ||
nlat = dim_eof(1) | ||
mlon = dim_eof(2) | ||
|
||
dnam_ts = getvardims(eof_ts) ; dimension names | ||
dnam_eof= getvardims(eof) ; used at end for meta data | ||
|
||
; ================================================================= | ||
; lag-0 and lag-1 matrices | ||
; ================================================================= | ||
|
||
if (get_ncl_version().eq."6.1.2") then ; bug in 6.1.2 | ||
cov0 = covcorm(eof_ts,(/1,0/)) ; lag-0 covariance matrix | ||
else | ||
cov0 = covcorm(eof_ts,(/0,1/)) ; lag-0 covariance matrix (n x n) | ||
end if | ||
; either | ||
cov1 = covcorm_xy(eof_ts, eof_ts, (/0,1,0/)) ; lag-1 | ||
;cov1 = covcorm_xy(eof_ts(:,0:ntim-2) \ ; alternative, brute force | ||
; ,eof_ts(:,1:ntim-1), (/0,0,0/)) | ||
;printVarSummary(cov1) | ||
|
||
; ================================================================= | ||
; matrix A contains information for evolution of the POP system. | ||
; POPs are eigenvectors of A. | ||
; ================================================================= | ||
|
||
cov0_inverse = inverse_matrix(cov0) | ||
A = cov1#inverse_matrix(cov0) ; [*][*] => neof x neof | ||
|
||
; ================================================================= | ||
; NCL 6.1.1 of dgeevx: evlr(2,2,N,N) ; (left(0)/right(1), real(0)/imag(1),:,:) | ||
; Eigenvalues are returned as attributes: eigi = evlr@eigi ; eigr = evlr@eigr | ||
; ================================================================= | ||
|
||
evlr = dgeevx_lapack(A, "B", "V", "V", "B", False) | ||
|
||
; ================================================================= | ||
; POP time series from eigenvalues and right eigenvectors | ||
; ================================================================= | ||
;PR = (/ evlr(1,0,:,:) /) ; right ev (1), real part (0) | ||
;PI = (/ evlr(1,1,:,:) /) ; right ev (1), imag part (1) | ||
; kPOP is what we want; use righteigenvector | ||
pr = (/ evlr(1,0,kPOP-1,:) /) ; right ev (1), real part (0), row 'kPOP-1' | ||
pi = (/ evlr(1,1,kPOP-1,:) /) ; right ev (1), imag part (1), row 'kPOP-1' | ||
|
||
z = inverse_matrix( (/ (/sum(pr*pr), sum(pr*pi)/) \ | ||
, (/sum(pr*pi), sum(pi*pi)/) /))#(/pr,pi/)#eof_ts | ||
|
||
; complex conjugate | ||
z = (/z(0,:), -z(1,:)/) ; real & imag series | ||
z = dim_rmvmean_n(z,1) | ||
mean = dim_avg_n(z,1) ; calculate mean | ||
stdev= dim_stddev_n(z,1) ; calculate stdev | ||
z = dim_standardize_n(z,1,1) ; standardize time series | ||
|
||
z!0 = "nPOP" ; add meta data | ||
z!1 = dnam_ts(1) | ||
z&nPOP = (/0,1/) | ||
z&$dnam_ts(1)$ = eof_ts&$dnam_ts(1)$ | ||
z@stdev = stdev | ||
z@mean = mean | ||
z@long_name = "POP timeseries" | ||
;printVarSummary(z) | ||
|
||
; ================================================================= | ||
; POP spatial patterns | ||
; ================================================================= | ||
|
||
zr = pr(0)*eof(0,:,:) ; construct POP spatial domain | ||
zi = pi(0)*eof(0,:,:) | ||
do j=1,neof-1 | ||
zr = zr + pr(j)*eof(j,:,:) | ||
zi = zi + pi(j)*eof(j,:,:) | ||
end do | ||
|
||
Z = (/zr*stdev(0), -zi*stdev(1)/) ; scale patterns by time series stdev | ||
|
||
Z!0 = "nPOP" ; add meta data | ||
Z!1 = dnam_eof(1) | ||
Z!2 = dnam_eof(2) | ||
|
||
Z&nPOP = (/0,1/) | ||
Z&$dnam_eof(1)$ = eof&$dnam_eof(1)$ | ||
Z&$dnam_eof(2)$ = eof&$dnam_eof(2)$ | ||
Z@long_name = "POP pattern" | ||
;printVarSummary(Z) | ||
|
||
; ================================================================= | ||
; return POP time series and POP spatial patterns as a | ||
; variable of type 'list' which contains 2 variables | ||
; ================================================================= | ||
|
||
return( [/z, Z/] ) ; this is type "list" | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,115 @@ | ||
;************************************************* | ||
; WRF static: panel different variables | ||
;************************************************ | ||
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/wrf/WRF_contributed.ncl" | ||
begin | ||
;************************************************ | ||
; open file and read in data | ||
;************************************************ | ||
f = addfile("static.wrfsi.nc", "r") | ||
;************************************************ | ||
; Read variables | ||
;************************************************ | ||
use = f->use(0,0,:,:) ; land use dominant category | ||
stl = f->stl(0,0,:,:) ; top layer (0-30cm) dom cat soiltype | ||
sbl = f->sbl(0,0,:,:) ; bottom layer (30-90cm) dom cat soiltype | ||
lat2d = f->lat(0,0,:,:) | ||
lon2d = f->lon(0,0,:,:) | ||
lsMask= f->lnd(0,0,:,:) ; land (1) water (0) mas | ||
|
||
;************************************************ | ||
; Use mask function to set all ocean areas to _FillValue | ||
;************************************************ | ||
use = mask(use,lsMask,1) | ||
stl = mask(stl,lsMask,1) | ||
sbl = mask(sbl,lsMask,1) | ||
|
||
;************************************************ | ||
; Associate 2D coordinates with variables for plotting | ||
;************************************************ | ||
use@lat2d = lat2d | ||
use@lon2d = lon2d | ||
stl@lat2d = lat2d | ||
stl@lon2d = lon2d | ||
sbl@lat2d = lat2d | ||
sbl@lon2d = lon2d | ||
|
||
;************************************************ | ||
; The file should be examined via: ncdump -v grid_type static.wrsi | ||
; This will print the print type. then enter below. | ||
;************************************************ | ||
projection = "mercator" | ||
|
||
;************************************************ | ||
; create plots | ||
;************************************************ | ||
wks = gsn_open_wks("ps" ,"WRF_static") ; ps,pdf,x11,ncgm,eps | ||
gsn_define_colormap(wks ,"BlAqGrYeOrReVi200"); choose colormap | ||
|
||
res = True ; plot mods desired | ||
res@gsnSpreadColors = True ; use full range of colormap | ||
res@cnFillOn = True ; color plot desired | ||
res@cnLinesOn = False ; turn off contour lines | ||
res@cnLineLabelsOn = False ; turn off contour labels | ||
res@cnLevelSpacingF = 1 ; manually specify interval | ||
res@cnFillMode = "RasterFill" ; activate raster mode | ||
res@lbLabelAutoStride = True ; let NCL figure lb stride | ||
|
||
;************************************************ | ||
; Turn on lat / lon labeling | ||
;************************************************ | ||
;;res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks | ||
|
||
dimll = dimsizes(lat2d) | ||
nlat = dimll(0) | ||
mlon = dimll(1) | ||
|
||
res@mpProjection = projection | ||
res@mpLimitMode = "Corners" | ||
res@mpLeftCornerLatF = lat2d(0,0) | ||
res@mpLeftCornerLonF = lon2d(0,0) | ||
res@mpRightCornerLatF = lat2d(nlat-1,mlon-1) | ||
res@mpRightCornerLonF = lon2d(nlat-1,mlon-1) | ||
|
||
res@mpCenterLonF = f->LoV ; set center logitude | ||
|
||
if (projection.eq."LambertConformal") then | ||
res@mpLambertParallel1F = f->Latin1 | ||
res@mpLambertParallel2F = f->Latin2 | ||
res@mpLambertMeridianF = f->LoV | ||
end if | ||
|
||
res@mpFillOn = False ; turn off map fill | ||
res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last | ||
res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries | ||
|
||
;;res@tfDoNDCOverlay = True ; True only for 'native' grid | ||
res@gsnAddCyclic = False ; data are not cyclic | ||
|
||
;************************************************ | ||
; allocate array for 3 plots | ||
;************************************************ | ||
plts = new (3,"graphic") | ||
|
||
;************************************************ | ||
; Tell NCL not to draw or advance frame for individual plots | ||
;************************************************ | ||
res@gsnDraw = False ; (a) do not draw | ||
res@gsnFrame = False ; (b) do not advance 'frame' | ||
|
||
plts(0) = gsn_csm_contour_map(wks,use,res) | ||
plts(1) = gsn_csm_contour_map(wks,stl,res) | ||
plts(2) = gsn_csm_contour_map(wks,sbl,res) | ||
;************************************************ | ||
; create panel: panel plots have their own set of resources | ||
;************************************************ | ||
resP = True ; modify the panel plot | ||
resP@txString = "Land Use and Soil Type" | ||
resP@gsnMaximize = True ; maximize panel area | ||
resP@gsnPanelRowSpec = True ; specify 1 top, 2 lower level | ||
gsn_panel(wks,plts,(/1,2/),resP) ; now draw as one plot | ||
|
||
end | ||
|
Oops, something went wrong.