::::::::::::::
pl.clim.ts.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=240
day=01
vlev=1
iys=2005
iye=2012
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
yearday_list="\
n10day_01_09 \
n10day_10_18 \
n10day_19_27 \
n10day_28_37 \
"
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} ${iys} ${iye} ${vlev} ${outdir}
done #yearday
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.ts.ncl
::::::::::::::
;
; Plot horizontal map of Monthly Climatology of Temperature and Salinity
;
; 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("export LANG=C; date") ; retrieve wall clock time
; Default values
dataset="dr080g"
dthr=240
iys=2005 ;1993
iye=2006 ;2012
ndy=9
nmo=3 ; Number of months in each file
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")
iys = stringtointeger(getenv("NCL_ARG_5"))
iye = stringtointeger(getenv("NCL_ARG_6"))
vlev = getenv("NCL_ARG_7")
outdir = getenv("NCL_ARG_8")
;iys=(iys_in)
;iye=floattointeger(iye_in)
;Set depth for plotting
klev=stringtointeger(vlev)-1
print("")
print("yearday = "+yearday)
;Set dimension
iy=iys
id=0
indir1=indir+"/"+dataset+"_"+iy+"/"+yearday+"/"
f =addfile(indir1 + "Tave" + datainfo + dthr + ".cdf" ,"r")
t = f->Tave(id,:,:,:)
lat =f->lat
lat_v=f->lat_v
lon =f->lon
lon_u=f->lon_u
depth=f->depth
time_in=f->time
z=floattointeger(depth(klev))
delete(f)
dim=dimsizes(t)
nyr=iye-iys+1
t_in=new((/dim(0),dim(1),dim(2)/),float)
t_in@_FillValue=t@_FillValue
s_in=new((/dim(0),dim(1),dim(2)/),float)
s_in@_FillValue=t@_FillValue
;printVarSummary(t_in)
;printVarSummary(s_in)
t_temp=new((/dim(0),dim(1),dim(2),nyr/),float)
t_temp@_FillValue=t@_FillValue
s_temp=new((/dim(0),dim(1),dim(2),nyr/),float)
s_temp@_FillValue=t@_FillValue
tclim=new((/nmo,dim(0),dim(1),dim(2)/),float)
tclim@_FillValue=t@_FillValue
sclim=new((/nmo,dim(0),dim(1),dim(2)/),float)
sclim@_FillValue=t@_FillValue
tclim@long_name = "Temperature (Monthly climatology)"
tclim@units = "degree C"
tclim!0="time"
tclim!1="depth"
tclim!2="lat"
tclim!3="lon"
tclim&lon=lon
tclim&lat=lat
sclim@long_name = "Salinity (Monthly climatology)"
sclim@units = "PSU"
sclim!0="time"
sclim!1="depth"
sclim!2="lat"
sclim!3="lon"
sclim&lon=lon
sclim&lat=lat
;printVarSummary(tclim)
;printVarSummary(sclim)
; Month
mm=new( 3, integer)
mmm=new( 3, string )
if (yearday .eq. "n10day_01_09") then
ms=1
mmm(0)="Jan"
mmm(1)="Feb"
mmm(2)="Mar"
end if
if (yearday .eq. "n10day_10_18") then
ms=4
mmm(0)="Apr"
mmm(1)="May"
mmm(2)="Jun"
end if
if (yearday .eq. "n10day_19_27") then
ms=7
mmm(0)="Jul"
mmm(1)="Aug"
mmm(2)="Sep"
end if
if (yearday .eq. "n10day_28_37") then
ms=10
mmm(0)="Oct"
mmm(1)="Nov"
mmm(2)="Dec"
end if
mm(0)=ms
mm(1)=ms + 1
mm(2)=ms + 2
;
; Plot map of T and S
;
;
; Set NCL resources
;
res = True
res@gsnDraw = False ; don't draw
res@gsnFrame = False ; don't advance frame
;Figure size
res@vpHeightF = 0.35
res@vpWidthF = 0.55
res@vpXF = 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
;
; Upper panel
;
res1=res
res1@vpYF = 0.9
res1@cnMinLevelValF = 22 ; set min contour level
res1@cnMaxLevelValF = 32 ; set max contour level
res1@cnLevelSpacingF = 1 ; set contour spacing
res1@gsnCenterString = "dr080 T (Climatology)" ; & ~F33~s ~B~q"
res1@gsnRightString = z + " m"
;
; Lower panel
;
res2=res
res2@vpYF = 0.45
res2@cnMinLevelValF = 30 ; set min contour level
res2@cnMaxLevelValF = 38 ; set max contour level
res2@cnLevelSpacingF = 1 ; set contour spacing
res2@gsnCenterString = "dr080 S (Climatology)" ; & ~F33~s ~B~q"
res2@gsnRightString = z + " m"
do im=0,2 ; LOOP FOR MONTH
do iy=0,nyr-1 ; LOOP FOR YEAR
yyyy=iy+iys
print(mmm(im) + " " + yyyy)
;
; Open input files
;
indir1=indir+"/"+dataset+"_"+yyyy+"/"+yearday+"/"
fT=addfile(indir1 + "Tave" + datainfo + dthr + ".cdf" ,"r")
fS=addfile(indir1 + "Save" + datainfo + dthr + ".cdf" ,"r")
do id=0,2 ; LOOP FOR BI-PENTAD DAYS
;
; Read input data
;
idx=im*3 + id
t_in = fT->Tave(idx,:,:,:)
s_in = fS->Save(idx,:,:,:)
;printVarSummary(t_in)
; For climatrogical average
t_temp(:,:,:,iy) = t_in(:,:,:)
s_temp(:,:,:,iy) = s_in(:,:,:)
end do ;id
end do ; iy
tclim(im,:,:,:)= dim_avg_n(t_temp, 3)
sclim(im,:,:,:)= dim_avg_n(s_temp, 3)
;
; Plotting
;
res1@gsnLeftString = mmm(im) + " ("+iys+"-"+iye+")" ; add the gsn titles
res2@gsnLeftString = mmm(im) + " ("+iys+"-"+iye+")"
ofle=outdir + "/" + "ECCO.dr080_TSclim." + sprinti("%0.5i",z) + "m." + sprinti("%0.2i",mm(im)) + "." + figtype
print("Figure: "+ofle)
wks = gsn_open_wks(figtype,ofle)
tcon = gsn_csm_contour_map(wks,tclim(im,klev,:,:),res1) ; Draw a contour plot.
draw(tcon)
scon = gsn_csm_contour_map(wks,sclim(im,klev,:,:),res2) ; Draw a contour plot.
draw(scon)
; Header and footer
;drawNDCGrid(wks)
txres = True ; text mods desired
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)
end do ; im
;
; Output netCDF files
;
time=new(1,float)
time!0="time"
time@units="Hours from Jan 1st"
do im=0,2
time=(mm(im)-0.5)*(30.0*24.0)
ofle1=indir+"/"+"ECCO.dr080g.Tclm."+iys+"-"+iye+"."+sprinti("%0.2i",mm(im))+".nc"
ofle2=indir+"/"+"ECCO.dr080g.Sclm."+iys+"-"+iye+"."+sprinti("%0.2i",mm(im))+".nc"
system("/bin/rm -f " + ofle1)
fout1 =addfile (ofle1, "c")
system("/bin/rm -f " + ofle2)
fout2 =addfile (ofle2, "c")
setfileoption(fout1,"DefineMode",True)
setfileoption(fout1,"MissingToFillValue",True)
setfileoption(fout2,"DefineMode",True)
setfileoption(fout2,"MissingToFillValue",True)
fAtt = True ; assign file attributes
fAtt@title = "ECCO dr080g (Monthly Climatology)"
fAtt@month = mmm(im)
fAtt@description = "Averaging period: "+iys+"-"+iye
fAtt@Conventions = "None"
fAtt@history= "Converted from MITGCM output on 07/24/2015\n"
fAtt@creation_date = systemfunc ("export LANG=C; date")
fAtt@hostname = systemfunc("hostname")
fAtt@cwd = systemfunc("pwd")
fAtt@script = scname
fAtt@input_dir = indir
fileattdef( fout1, fAtt ) ; copy file attributes
fileattdef( fout2, fAtt ) ; copy file attributes
ntim = 1; dimsizes(time) ; get dimension sizes
nlev = dimsizes(depth)
nlat = dimsizes(lat)
nlon = dimsizes(lon)
dimNames = (/"time", "lat", "lon", "depth"/)
dimSizes = (/ ntim, nlat, nlon, nlev /)
dimUnlim = (/False, False, False, False /)
filedimdef(fout1,dimNames,dimSizes,dimUnlim)
filedimdef(fout2,dimNames,dimSizes,dimUnlim)
filevardef(fout1, "time" ,typeof(time),getvardims(time))
filevardef(fout1, "lat" ,typeof(lat), getvardims(lat))
filevardef(fout1, "lon" ,typeof(lon), getvardims(lon))
filevardef(fout1, "depth" ,typeof(depth),getvardims(depth) )
filevardef(fout1, "tclim" ,typeof(tclim),getvardims(tclim(im,:,:,:)))
filevardef(fout2, "time" ,typeof(time),getvardims(time))
filevardef(fout2, "lat" ,typeof(lat), getvardims(lat))
filevardef(fout2, "lon" ,typeof(lon), getvardims(lon))
filevardef(fout2, "depth" ,typeof(depth),getvardims(depth) )
filevardef(fout2, "sclim" ,typeof(sclim),getvardims(sclim(im,:,:,:)))
filevarattdef(fout1,"time" ,time)
filevarattdef(fout1,"lat" ,lat)
filevarattdef(fout1,"lon" ,lon)
filevarattdef(fout1,"depth" ,depth)
filevarattdef(fout1,"tclim" ,tclim)
filevarattdef(fout2,"time" ,time)
filevarattdef(fout2,"lat" ,lat)
filevarattdef(fout2,"lon" ,lon)
filevarattdef(fout2,"depth" ,depth)
filevarattdef(fout2,"sclim" ,sclim)
fout1->time=(/time/)
fout1->lat =(/lat/)
fout1->lon = (/lon/)
fout1->depth=(/depth/)
fout1->tclim = (/tclim(im,:,:,:)/)
fout2->time=(/time/)
fout2->lat =(/lat/)
fout2->lon = (/lon/)
fout2->depth=(/depth/)
fout2->sclim = (/sclim(im,:,:,:)/)
print("Output file:"+ofle1)
print("Output file:"+ofle2)
print("")
end do ;im
end