Plot horizontal map of isothermal-layer depth

::::::::::::::

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