Embedded Files
時間変化項(Tendency)を10日平均のデータから求めている。
dtは20日。
[2015年 8月 18日 火曜日 12:15:51 JST]
[~/2015.Ibnu.Indian.Ocean/ECCO.dr080g/Map.Heat.Budget]
[am@aofd165]
How to run
How to run
$ sh pl.heat.budget.all.run.sh
Source codes
Source codes
::::::::::::::
pl.heat.budget.all.run.sh
::::::::::::::
#!/bin/bash
exe=pl.heat.budget.run.sh
if [ ! -f $exe ]; then
echo
echo Error in $0 : No such file, $exe
echo
exit 1
fi
yyyy=1993
yyyye=2014
ydlist="\
n10day_01_09 \
n10day_10_18 \
n10day_19_27 \
n10day_28_37 \
"
vlev=1
while [ $yyyy -le $yyyye ]; do
for yearday in $ydlist; do
if [ $yearday = "n10day_28_37" ]; then
de=4
else
de=3
fi
day=1
while [ $day -le $de ]; do
$exe $yyyy $yearday $day $vlev
day=$(expr $day + 1)
done #day
done # yearday
yyyy=$(expr $yyyy + 1)
done # yyyy
exit 0
::::::::::::::
pl.heat.budget.run.sh
::::::::::::::
#!/bin/sh
export LANC=C
exe=./runncl.sh
ncl=pl.heat.budget.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
# Read arguments
yyyy=$1
yearday=$2
dthr=720
day=$3
vlev=$4
indir=/work4/data/ECCO.dr080g/${dataset}_${yyyy}/${yearday}
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)/${yyyy}
if [ ! -d $outdir ]; then
mkdir -vp $outdir
fi
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} ${datainfo} ${dthr} ${day} ${vlev} ${outdir}
echo
echo Done $(basename $0).
echo
::::::::::::::
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 $1
::::::::::::::
pl.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
;
; Read command line arguments
;
scname = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
datainfo = getenv("NCL_ARG_3")
dthr = getenv("NCL_ARG_4")
day = getenv("NCL_ARG_5")
vlev = getenv("NCL_ARG_6")
outdir = getenv("NCL_ARG_7")
day_no=stringtointeger(day)
day_no=day_no-1
if ( day_no .eq. 0 ) then
day_no_s1=0
day_no_s2=2
end if
if ( day_no .eq. 1 ) then
day_no_s1=3
day_no_s2=5
end if
if ( day_no .eq. 2 ) then
day_no_s1=6
day_no_s2=8
end if
if ( day_no .eq. 3 ) then
day_no_s1=7
day_no_s2=9
end if
dt_budget=20*(24.0*3600.0)
;Set depth
klev=stringtointeger(vlev)-1
;
; Open input files
;
fExt =addfile(indir + "/gtExtave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fRelax =addfile(indir + "/gtRelaxave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fKF =addfile(indir + "/gtKFave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fU =addfile(indir + "/gtUave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fV =addfile(indir + "/gtVave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fW =addfile(indir + "/gtWave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fXdiff =addfile(indir + "/gtXdiffave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fYdiff =addfile(indir + "/gtYdiffave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fZdiff =addfile(indir + "/gtZdiffave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fZdiffGM=addfile(indir + "/gtZdiffGMave_" + datainfo + "_" + dthr + ".cdf" ,"r")
fKPP =addfile(indir + "/gtKPPave_" + datainfo + "_" + dthr + ".cdf" ,"r")
dthr_s =240
fT =addfile(indir + "/Tave_" + datainfo + "_" + dthr_s + ".cdf", "r")
;
; Read input data
;
ext = fExt->gtExtave(day_no,klev,:,:)
relax = fRelax->gtRelaxave(day_no,:,:)
kf = fKF->gtKFave(day_no,klev,:,:)
uadv = fU->gtUave(day_no,klev,:,:)
vadv = fV->gtVave(day_no,klev,:,:)
wadv = fW->gtWave(day_no,klev,:,:)
xdif = fXdiff->gtXdiffave(day_no,klev,:,:)
ydif = fYdiff->gtYdiffave(day_no,klev,:,:)
zdif = fZdiff->gtZdiffave(day_no,klev,:,:)
zdifG = fZdiffGM->gtZdiffGMave(day_no,klev,:,:)
zdifK = fKPP->gtKPPave(day_no,klev,:,:)
time=fU->time(day_no)
z1d=fU->depth
T1 = fT->Tave(day_no_s1,klev,:,:)
T2 = fT->Tave(day_no_s2,klev,:,:)
time_s1= fT->time(day_no_s1)
time_s2= fT->time(day_no_s2)
;
; 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/24.0 - 30
julim1=julian0 + time/24.0
greg0=jul2greg(julim0)
greg1=jul2greg(julim1)
date0=greg0(0)+sprinti("%0.2i",greg0(1))+sprinti("%0.2i",greg0(2))
date1=greg1(0)+sprinti("%0.2i",greg1(1))+sprinti("%0.2i",greg1(2))
dateprt0=greg0(0) + "/" +sprinti("%0.2i",greg0(1)) + "/" +sprinti("%0.2i",greg0(2))
dateprt1=greg1(0) + "/" +sprinti("%0.2i",greg1(1)) + "/" +sprinti("%0.2i",greg1(2))
dateprt=dateprt0 + "-" + dateprt1
julim_s1=julian0 + time_s1/24.0
julim_s2=julian0 + time_s2/24.0
greg_s1=jul2greg(julim_s1)
greg_s2=jul2greg(julim_s2)
dateprt_s1=greg_s1(0) + "/" +sprinti("%0.2i",greg_s1(1)) + "/" +sprinti("%0.2i",greg_s1(2))
dateprt_s2=greg_s1(0) + "/" +sprinti("%0.2i",greg_s2(1)) + "/" +sprinti("%0.2i",greg_s2(2))
dateprt_s=dateprt_s1 + "-" + dateprt_s2
;Set depth
z=floattointeger(z1d(klev))
;
; Compute budget terms
;
tend=T1
sfcflx=ext
hadvdif=uadv
vdif=zdif
rhs=ext
tend=(T2-T1)/dt_budget
sfcflx=ext+ relax
hadvdif=uadv + vadv + xdif + ydif
vdif=zdif + zdifG + zdifK
rhs=sfcflx + hadvdif+ wadv + vdif
;
; Set output file name
;
filetype="png" ;"ps"
ofle=outdir + "/" + "ECCO.dr080.Heat.Budget." + sprinti("%0.5i",z) + "m." + date0 + "-" + date1 + "." + filetype
wks = gsn_open_wks(filetype,ofle)
;
; 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 = -1.5e-6 ; set min contour level
res1@cnMaxLevelValF = 1.5e-6 ; set max contour level
res1@cnLevelSpacingF = 5.0e-7 ; set contour spacing
res1@gsnLeftString = dateprt_s ; add the gsn titles
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@gsnLeftString = dateprt ; add the gsn titles
res2@gsnRightString = z + " m"
;
; Upper left panel
;
res1@vpXF = 0.05
res1@vpYF = 0.85
res1@gsnCenterString = "LHS (Tendency)"
con = gsn_csm_contour_map(wks,tend,res1) ; Draw a contour plot.
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,res1) ; Draw a contour plot.
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@gsnCenterString = "Sfc Flux + Relax"
res2@gsnRightString = z + " m"
con = gsn_csm_contour_map(wks,sfcflx,res2) ; Draw a contour plot.
plot=con
draw(plot)
;
; Middle right panel
;
res2@vpXF = 0.52
res2@gsnCenterString = "H Adv + H Diff"
con = gsn_csm_contour_map(wks,hadvdif,res2) ; Draw a contour plot.
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,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,res2) ; Draw a contour plot.
plot=con
draw(plot)
;
; Header and footer
;
;drawNDCGrid(wks)
txres = True ; text mods desired
txres@txFontHeightF = 0.03 ; font smaller. default big
gsn_text_ndc(wks,dateprt,0.5,0.95,txres)
txres@txFontHeightF = 0.01 ; font smaller. default big
gsn_text_ndc(wks,outdir,0.5,0.99,txres)
text=systemfunc("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
Google Sites
Report abuse