::::::::::::::
pl.ILD.all.run.sh
::::::::::::::
#!/bin/bash
dt=0.3 # Tempearture threshold [degC]
exe=$(basename $0 .all.run.sh ).run.sh
if [ ! -f $exe ]; then echo "Error in $0, No such file, $exe"; exit0; fi
gen_fort_lib=wrapit.run.sh
if [ ! -f $gen_fort_lib ]; then echo "Error in $0, No such file, $exe"; exit0; fi
$gen_fort_lib
m=1
me=12
while [ $m -le $me ]; do
$exe $m $dt
m=$(expr $m + 1)
done
exit 0
::::::::::::::
pl.ILD.run.sh
::::::::::::::
#!/bin/sh
export LANC=C
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
# Default value
dataset=WOA
yyyy="2005-2012"
version="_01v2"
m=1
vlev=1
dt=0.3
m=$1
dt=$2
echo $m
indir=/work4/data/${dataset}
if [ ! -d $indir ]; then
echo
echo Error in $0 : No such directory, $indir
echo
exit 1
fi
outdir=${HOME}/2015.Ibnu.Indian.Ocean/Fig/${dataset}/$(basename $(pwd))/\
$(basename $0 .sh)/${yyyy}
if [ ! -d $outdir ]; then
mkdir -vp $outdir
fi
mm=$(printf %02d $m)
if [ $m -eq 1 ]; then mmm="Jan"; fi
if [ $m -eq 2 ]; then mmm="Feb"; fi
if [ $m -eq 3 ]; then mmm="Mar"; fi
if [ $m -eq 4 ]; then mmm="Apr"; fi
if [ $m -eq 5 ]; then mmm="May"; fi
if [ $m -eq 6 ]; then mmm="Jun"; fi
if [ $m -eq 7 ]; then mmm="Jul"; fi
if [ $m -eq 8 ]; then mmm="Aug"; fi
if [ $m -eq 9 ]; then mmm="Sep"; fi
if [ $m -eq 10 ]; then mmm="Oct"; fi
if [ $m -eq 11 ]; then mmm="Nov"; fi
if [ $m -eq 12 ]; then mmm="Dec"; fi
in[0]=${indir}"/woa13_A5B2_t"${mm}"${version}.nc"
#in[1]=${indir}"/woa13_A5B2_s"${mm}"${version}.nc"
#$exe $ncl ${in[0]} ${in[1]} ${yyyy} ${m} ${mmm} ${dt} ${outdir}
$exe $ncl ${in[0]} ${yyyy} ${m} ${mmm} ${dt} ${outdir}
echo
echo Done $(basename $0).
echo
exit 0
::::::::::::::
wrapit.run.sh
::::::::::::::
#!/bin/bash
f90=find_ild.f90
f90stub=$(basename $f90 .f90).stub
cat <<EOF>$f90
subroutine find_ild(nx,ny,nz,t3d,t2d,dt,z1d,fillvalue,ild)
integer,intent(in)::nx,ny,nz
real,intent(in)::t3d(nx,ny,nz),t2d(nx,ny)
real,intent(in)::dt
real,intent(in)::z1d(nz)
real,intent(in)::fillvalue
real,intent(inout)::ild(nx,ny)
! z1d(k): Model's vertical coordinate
! ild=isothermal layer depth
do j=1,ny
do i=1,nx
ild(i,j)=fillvalue ! default = undefined
xs= t2d(i,j) - dt
do k=1,nz-1
if(t3d(i,j,k+1) <= xs .and. xs < t3d(i,j,k) )then
y1=z1d(k+1)
y2=z1d(k)
x1=t3d(i,j,K+1)
x2=t3d(i,j,k)
!Linear Interpolation
ys=(y2-y1)/(x2-x1)*(xs-x1) + y1
ild(i,j)=ys
end if
end do !k
! Otherwise, the water column should be vertically homogeneous.
! ild(i,j)=fillvalue
end do !i
end do !j
end subroutine
EOF
cat <<EOF>$f90stub
C NCLFORTSTART
subroutine find_ild(nx,ny,nz,t3d,t2d,dt,z1d,fillvalue,ild)
integer nx,ny,nz
real t3d(nx,ny,nz)
real t2d(nx,ny)
real dt
real z1d(nz)
real fillvalue
real ild(nx,ny)
C NCLEND
EOF
export LANG=C
echo
date
echo
echo Fortran Source file: $f90
echo Stub file: $f90stub
echo
WRAPIT $f90stub $f90
if [ $? -ne 0 ]; then
echo "***"
echo "*** ERROR in $0 while uunning WRAPIT."
echo "***"
exit 1
else
so_file=$(basename $f90 .f90).so
echo Shared Library: $so_file
ls -lh $so_file
echo
fi
date
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.ILD.ncl
::::::::::::::
;
; Plot horizontal map of isothermal-layer depth (ILD) of WOA13F V2 dataset
;
; Data site
; https://www.nodc.noaa.gov/cgi-bin/OC5/woa13/woa13.pl
;
; File name format
; woa13_A5B2_t01_01v2.nc
;
; A5 = 2005
; B2 = 2012
;
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"
external FIND_ILD "./find_ild.so"
begin
wallClock1 = systemfunc("date") ; retrieve wall clock time
;
; default value
;
vlev=2
;
; Read command line arguments
;
scname = getenv("NCL_ARG_1")
infle1 = getenv("NCL_ARG_2")
yyyy = getenv("NCL_ARG_3")
m = getenv("NCL_ARG_4")
mmm = getenv("NCL_ARG_5")
dt_in = getenv("NCL_ARG_6")
outdir = getenv("NCL_ARG_7")
dt=stringtofloat(dt_in)
;infle2 = getenv("NCL_ARG_3")
;
; Read input data
;
ft=addfile(infle1,"r")
;fs=addfile(infle2,"r")
tglobal3d=ft->t_an
;printVarSummary(tglobal3d)
;sglobal3d=fs->s_an
z1d=ft->depth
lat=ft->lat
lon=ft->lon
;printVarSummary(tglobal3d)
t3d=tglobal3d(time|0,depth|:,{lat|-30:30},{lon|30:120})
;
; Find ILD
;
dims=dimsizes(t3d)
im=dims(2)
jm=dims(1)
km=dims(0)
print(im+" "+jm+" "+km)
ild=t3d(0,:,:)
ild@coordinates:="lat lon"
delete(ild@depth)
ild!0:="lat"
ild!1:="lon"
printVarSummary(t3d)
printVarSummary(ild)
;Set depth
klev=stringtointeger(vlev)-1
z=floattointeger(z1d(klev))
;print(z)
t2d = t3d(klev,:,:)
t2d@coordinates:="lat lon"
delete(t2d@depth)
printVarSummary(t2d)
FIND_ILD::find_ild(im,jm,km,t3d,t2d,dt,z1d,t3d@_FillValue,ild)
;print(ild)
;
; Set year and month
;
dateprt= mmm + " (" + yyyy +")"
;
; Set output file name
;
filetype="png" ;"ps"
mm=sprinti("%0.2i",stringtointeger(m))
ofle=outdir + "/" + "WOA13_ILD." + yyyy + "." + mm + "." + 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
;Figure size
res@vpHeightF = 0.35
res@vpWidthF = 0.55
res@vpXF = 0.2
;Map
res@mpFillOn = True
res@gsnAddCyclic = False
;Contour
res@cnFillOn = True ; turn on color fill
res@cnLineColor = "white"
res@lbOrientation = "vertical"
res@cnLevelSelectionMode = "ManualLevels" ; manually set the contour levels with the following 3 resources
res@gsnContourLineThicknessesScale=1.5
;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 = 10 ; set min contour level
res1@cnMaxLevelValF = 120 ; set max contour level
res1@cnLevelSpacingF = 10 ; set contour spacing
res1@gsnLeftString = dateprt ; add the gsn titles
res1@gsnCenterString = "ILD [m]"
res1@gsnRightString = "WOA13"
plot1 = gsn_csm_contour_map(wks,ild,res1) ; Draw a contour plot.
; now assign plotA to array
draw(plot1)
delete(res1)
txres = True ; text mods desired
txres@txFontHeightF = 0.015 ; font smaller. default big
txres@txBackgroundFillColor = "white"
txres@txFontColor="black"
gsn_text_ndc(wks,dateprt,0.65,0.59,txres)
;
;;
;; 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@gsnLeftString = dateprt ; add the gsn titles
;res2@gsnCenterString = "S"
;res2@gsnRightString = z + " m"
;
;plot2 = gsn_csm_contour_map(wks,s3d(0,klev,:,:),res2) ; Draw a contour plot.
;
;draw(plot2)
;delete(res2)
;
;gsn_text_ndc(wks,dateprt,0.65,0.14,txres)
;
;
;
; 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("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