[2015年 8月 28日 金曜日 09:34:16 JST]
[~/2015.Ibnu.Indian.Ocean/ECCO.dr080g/Clim.heat.budget]
[am@aofd165]
::::::::::::::
pl.clim.heat.budget.run.sh
::::::::::::::
#!/bin/sh
export LANC=C
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
if [ ! -f $ncl ]; then
echo
echo Error in $0 : No such file, $ncl
echo
exit 1
fi
# Default value
dataset=dr080g
yyyy=1993
yearday=n10day_01_09
dthr=720
day=01
vlev=1
indir=/work4/data/ECCO.dr080g
if [ ! -d $indir ]; then
echo
echo Error in $0 : No such directory, $indir
echo
exit 1
fi
outdir=${HOME}/2015.Ibnu.Indian.Ocean/Fig/$(basename $(pwd))/\
$(basename $0 .sh)/
if [ ! -d $outdir ]; then mkdir -vp $outdir ; fi
vlev_list="\
10 \
"
yearday_list="\
n10day_01_09 \
n10day_10_18 \
n10day_19_27 \
n10day_28_37 \
"
for vlev in $vlev_list; do
for yearday in $yearday_list; do
if [ $yearday = "n10day_01_09" ]; then
datainfo="_08_08.00001_02160_"
elif [ $yearday = "n10day_10_18" ]; then
datainfo="_08_08.02160_04320_"
elif [ $yearday = "n10day_19_27" ]; then
datainfo="_08_08.04320_06480_"
elif [ $yearday = "n10day_28_37" ]; then
datainfo="_08_08.06480_08880_"
fi
$exe $ncl ${indir} ${yearday} ${datainfo} ${vlev} ${outdir}
done #yearday
done #vlev
echo
echo Done $(basename $0).
echo
exit 0
::::::::::::::
runncl.sh
::::::::::::::
#!/bin/bash
#
# Universal wrapper script for ncl.
# Pass arguments from the command line to environment variables
#
# version 0.1, Thierry Corti, C2SM ETH Zurich
#
E_BADARGS=65
if [ ! -n "$1" ]
then
echo "Usage: `basename $0` script.ncl argument1 argument2 etc."
exit $E_BADARGS
fi
# save number of arguments to environment variable NCL_N_ARG
export NCL_N_ARGS=$#
# save command line arguments to environment variable NCL_ARG_#
for ((index=1; index<=$#; index++))
do
eval export NCL_ARG_$index=\$$index
done
# run ncl
ncl -Q -n $1
exit 0
::::::::::::::
pl.clim.heat.budget.ncl
::::::::::::::
;
; Plot horizontal map of heat budget terms
;
; Refer to http://dx.doi.org/10.1002/2014JC010336 for details on dr080
;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
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"
begin
wallClock1 = systemfunc("date") ; retrieve wall clock time
; Default values
dataset="dr080g"
dthr=720
iys=1993
iye=2014
ndy=3
figtype="png" ;"ps"
;
; Read command line arguments
;
scname = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
yearday = getenv("NCL_ARG_3")
datainfo = getenv("NCL_ARG_4")
vlev = getenv("NCL_ARG_5")
outdir = getenv("NCL_ARG_6")
;Set depth
klev=stringtointeger(vlev)-1
print("")
print("yearday = "+yearday)
;Set dimension
iy=iys
id=0
indir1=indir+"/"+dataset+"_"+iy+"/"+yearday+"/"
f =addfile(indir1 + "gtTndave" + datainfo + dthr + ".cdf" ,"r")
tnd = f->gtTndave(id,klev,:,:)
dim=dimsizes(tnd)
nyr=iye-iys+1
tnd_in=new((/ndy,dim(0),dim(1),nyr/),float)
tnd_in@_FillValue=tnd@_FillValue
sfcflx_in=new((/ndy,dim(0),dim(1),nyr/),float)
sfcflx_in@_FillValue=tnd@_FillValue
ext_in=new((/ndy,dim(0),dim(1),nyr/),float)
ext_in@_FillValue=tnd@_FillValue
hadvdif_in=new((/ndy,dim(0),dim(1),nyr/),float)
hadvdif_in@_FillValue=tnd@_FillValue
vdif_in=new((/ndy,dim(0),dim(1),nyr/),float)
vdif_in@_FillValue=tnd@_FillValue
wadv_in=new((/ndy,dim(0),dim(1),nyr/),float)
wadv_in@_FillValue=tnd@_FillValue
rhs_in=new((/ndy,dim(0),dim(1),nyr/),float)
rhs_in@_FillValue=tnd@_FillValue
lat =f->lat
lat_v=f->lat_v
lon =f->lon
lon_u=f->lon_u
depth=f->depth
time=f->time
delete(f)
do iy=0,nyr-1
yyyy=iy+iys
print("Year = " + yyyy)
;
; Open input files
;
indir1=indir+"/"+dataset+"_"+yyyy+"/"+yearday+"/"
fTnd =addfile(indir1 + "gtTndave" + datainfo + dthr + ".cdf" ,"r")
fExt =addfile(indir1 + "gtExtave" + datainfo + dthr + ".cdf" ,"r")
fRelax =addfile(indir1 + "gtRelaxave" + datainfo + dthr + ".cdf" ,"r")
fKF =addfile(indir1 + "gtKFave" + datainfo + dthr + ".cdf" ,"r")
fU =addfile(indir1 + "gtUave" + datainfo + dthr + ".cdf" ,"r")
fV =addfile(indir1 + "gtVave" + datainfo + dthr + ".cdf" ,"r")
fW =addfile(indir1 + "gtWave" + datainfo + dthr + ".cdf" ,"r")
fXdiff =addfile(indir1 + "gtXdiffave" + datainfo + dthr + ".cdf" ,"r")
fYdiff =addfile(indir1 + "gtYdiffave" + datainfo + dthr + ".cdf" ,"r")
fZdiff =addfile(indir1 + "gtZdiffave" + datainfo + dthr + ".cdf" ,"r")
fZdiffGM=addfile(indir1 + "gtZdiffGMave" + datainfo + dthr + ".cdf" ,"r")
fKPP =addfile(indir1 + "gtKPPave" + datainfo + dthr + ".cdf" ,"r")
z1d=fU->depth
;Set depth
z=floattointeger(z1d(klev))
do id=0,2
;
; Read input data
;
tnd = fTnd->gtTndave(id,klev,:,:)
ext = fExt->gtExtave(id,klev,:,:)
relax = fRelax->gtRelaxave(id,:,:)
kf = fKF->gtKFave(id,klev,:,:)
uadv = fU->gtUave(id,klev,:,:)
vadv = fV->gtVave(id,klev,:,:)
wadv = fW->gtWave(id,klev,:,:)
xdif = fXdiff->gtXdiffave(id,klev,:,:)
ydif = fYdiff->gtYdiffave(id,klev,:,:)
zdif = fZdiff->gtZdiffave(id,klev,:,:)
zdifG = fZdiffGM->gtZdiffGMave(id,klev,:,:)
zdifK = fKPP->gtKPPave(id,klev,:,:)
;
; Compute budget terms
;
sfcflx=ext
hadvdif=uadv
vdif=zdif
rhs=ext
if(klev .eq. 0)then
sfcflx=ext + relax
hadvdif=uadv + vadv + xdif + ydif
vdif=zdif + zdifG + zdifK
rhs=sfcflx + hadvdif+ wadv + vdif + kf
else
sfcflx=0.0
hadvdif=uadv + vadv + xdif + ydif
vdif=zdif + zdifG + zdifK
rhs=ext + hadvdif+ wadv + vdif + kf
end if
; For climatrogical averaging
tnd_in(id,:,:,iy) = tnd(:,:)
sfcflx_in(id,:,:,iy) = sfcflx(:,:)
ext_in(id,:,:,iy) = ext(:,:)
hadvdif_in(id,:,:,iy) = hadvdif(:,:)
wadv_in(id,:,:,iy) = wadv(:,:)
vdif_in(id,:,:,iy) = vdif(:,:)
rhs_in(id,:,:,iy) = rhs(:,:)
end do ;id
end do ; iy
tnd_ave = dim_avg_n(tnd_in, 3)
sfcflx_ave = dim_avg_n(sfcflx_in, 3)
ext_ave = dim_avg_n(ext_in, 3)
hadvdif_ave= dim_avg_n(hadvdif_in,3)
wadv_ave = dim_avg_n(wadv_in, 3)
vdif_ave = dim_avg_n(vdif_in, 3)
rhs_ave = dim_avg_n(rhs_in, 3)
tnd_ave@long_name = "Temperature tendency (Climatological average)"
tnd_ave@units = "DEGREES/S"
tnd_ave!0="time"
tnd_ave!1="lat"
tnd_ave!2="lon"
tnd_ave&lon = lon
tnd_ave&lat = lat
sfcflx_ave@long_name = "Temperature tendency due to surface flux + relaxation (Climatological average)"
sfcflx_ave@units = "DEGREES/S"
sfcflx_ave!0="time"
sfcflx_ave!1="lat"
sfcflx_ave!2="lon"
sfcflx_ave&lon = lon
sfcflx_ave&lat = lat
ext_ave@long_name = "Temperature tendency due to external forcing (Climatological average)"
ext_ave@units = "DEGREES/S"
ext_ave!0="time"
ext_ave!1="lat"
ext_ave!2="lon"
ext_ave&lon = lon
ext_ave&lat = lat
hadvdif_ave@long_name = "Temperature tendency due to horizontal advection and diffusion (Climatological average)"
hadvdif_ave@units = "DEGREES/S"
hadvdif_ave!0="time"
hadvdif_ave!1="lat"
hadvdif_ave!2="lon"
hadvdif_ave&lon = lon
hadvdif_ave&lat = lat
wadv_ave@long_name = "Temperature tendency due to vertical advection (Climatological average)"
wadv_ave@units = "DEGREES/S"
wadv_ave!0="time"
wadv_ave!1="lat"
wadv_ave!2="lon"
wadv_ave&lon = lon
wadv_ave&lat = lat
vdif_ave@long_name = "Temperature tendency due to vertical diffusion (Climatological average)"
vdif_ave@units = "DEGREES/S"
vdif_ave!0="time"
vdif_ave!1="lat"
vdif_ave!2="lon"
vdif_ave&lon = lon
vdif_ave&lat = lat
rhs_ave@long_name = "RHS of heat budget equation (Climatological average)"
rhs_ave@units = "DEGREES/S"
rhs_ave!0="time"
rhs_ave!1="lat"
rhs_ave!2="lon"
rhs_ave&lon = lon
rhs_ave&lat = lat
;
; Set NCL resources
;
res = True
res@gsnDraw = False ; don't draw
res@gsnFrame = False ; don't advance frame
;res@gsnMaximize = True
res@gsnPaperOrientation = "landscape"
;Figure size
res@vpHeightF = 0.2; ;0.35
res@vpWidthF = 0.4; 0.55
res@vpXF = 0.05; 0.2
;Map
res@mpFillOn = True
;Contour
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False
res@lbOrientation = "vertical"
res@cnLevelSelectionMode = "ManualLevels" ; manually set the contour levels with the following 3 resources
;Choose subregion
res@mpMinLatF =-20
res@mpMaxLatF = 20
res@mpMinLonF = 40
res@mpMaxLonF =110
res@mpDataBaseVersion = "Ncarg4_1" ; use GMT coastline
res1=res
res1@cnMinLevelValF = -1e-6 ; set min contour level
res1@cnMaxLevelValF = 1e-6 ; set max contour level
res1@cnLevelSpacingF = 2e-7 ; set contour spacing
res1@gsnRightString = z + " m"
res2=res
res2@cnMinLevelValF = -3.0e-5 ; set min contour level
res2@cnMaxLevelValF = 3.0e-5 ; set max contour level
res2@cnLevelSpacingF = 5.0e-6 ; set contour spacing
res2@gsnRightString = z + " m"
do id=0,ndy-1
;
; Set date
;
julian0=greg2jul(1970,01,01,00)
;
; float time(time) ;
; time:units = "hours since 1970-01-01 00:00:0.0" ;
;
julim0=julian0 + time(id)/24.0 - 30
julim1=julian0 + time(id)/24.0
greg0=jul2greg(julim0)
greg1=jul2greg(julim1)
date0=sprinti("%0.2i",greg0(1))+sprinti("%0.2i",greg0(2))
date1=sprinti("%0.2i",greg1(1))+sprinti("%0.2i",greg1(2))
dateprt0=sprinti("%0.2i",greg0(1)) + "/" +sprinti("%0.2i",greg0(2))
dateprt1=sprinti("%0.2i",greg1(1)) + "/" +sprinti("%0.2i",greg1(2))
dateprt=dateprt0 + "-" + dateprt1
print("date = "+ dateprt)
;
; Set output file name
;
ofle=outdir + "/" + "ECCO.dr080.Clim.heat.budget." + sprinti("%0.5i",z) + "m." + date0 + "-" + date1 + "." + figtype
wks = gsn_open_wks(figtype,ofle)
res1@gsnLeftString = dateprt ; add the gsn titles
res2@gsnLeftString = dateprt ; add the gsn titles
;
; Upper left panel
;
res1@vpXF = 0.05
res1@vpYF = 0.85
res1@gsnCenterString = "LHS (Tendency)"
con = gsn_csm_contour_map(wks,tnd_ave(id,:,:),res1)
plot=con
draw(plot)
;
; Upper right panel
;
res1@vpXF = 0.52
res1@vpYF = 0.85
res1@gsnLeftString = dateprt ; add the gsn titles
res1@gsnCenterString = "RHS"
con = gsn_csm_contour_map(wks,rhs_ave(id,:,:),res1)
plot=con
draw(plot)
;
; Middle left panel
;
res2@vpXF = 0.05
res2@vpYF = 0.60
res2@cnMinLevelValF = -3.0e-5 ; set min contour level
res2@cnMaxLevelValF = 3.0e-5 ; set max contour level
res2@cnLevelSpacingF = 5.0e-6 ; set contour spacing
res2@gsnLeftString = dateprt ; add the gsn titles
res2@gsnRightString = z + " m"
if (klev .eq. 0)then
res2@gsnCenterString = "Ext + Relax"
con = gsn_csm_contour_map(wks,sfcflx_ave(id,:,:),res2)
else
res2@gsnCenterString = "Ext"
con = gsn_csm_contour_map(wks,ext_ave(id,:,:),res2)
end if
plot=con
draw(plot)
;
; Middle right panel
;
res2@vpXF = 0.52
res2@gsnCenterString = "H Adv + H Diff"
con = gsn_csm_contour_map(wks,hadvdif_ave(id,:,:),res2)
plot=con
draw(plot)
;
; Lower left panel
;
res2@vpXF = 0.05
res2@vpYF= 0.35
res2@gsnCenterString = "V ADV"
con = gsn_csm_contour_map(wks,wadv_ave(id,:,:),res2) ; Draw a contour plot.
plot=con
draw(plot)
;
; Lower right panel
;
res2@vpXF = 0.52
res2@gsnCenterString = "V DIFF"
con = gsn_csm_contour_map(wks,vdif_ave(id,:,:),res2) ; Draw a contour plot.
plot=con
draw(plot)
;
; Header and footer
;
;drawNDCGrid(wks)
txres = True ; text mods desired
txres@txFontHeightF = 0.025 ; font smaller. default big
gsn_text_ndc(wks,"Climatology (" + dateprt + ")",0.5,0.95,txres)
txres@txFontHeightF = 0.02 ; font smaller. default big
gsn_text_ndc(wks,"Depth = " + z + " m",0.5,0.90,txres)
txres@txFontHeightF = 0.01 ; font smaller. default big
gsn_text_ndc(wks,outdir,0.5,0.99,txres)
text=systemfunc("export LANG=C; date")
gsn_text_ndc(wks,text,0.2,0.010,txres)
text=systemfunc("pwd")+"/"+scname
gsn_text_ndc(wks,text,0.5,0.03,txres)
frame(wks)
print("")
print("Output: "+ ofle)
print("")
end do ;id
end