::::::::::::::
pl.MLD.ILD.BLT.all.run.sh
::::::::::::::
#!/bin/bash
dt=0.3 # Tempearture threshold [degC]
fort=mixed_layer.f90
exe=$(basename $0 .all.run.sh ).run.sh
if [ ! -f $exe ]; then echo "Error in $0, No such file, $exe"; exit 1; fi
gen_fort_lib=wrapit.sh
if [ ! -f $gen_fort_lib ]; then echo "Error in $0, No such file, $exe"; exit 1; fi
$gen_fort_lib $fort
if [ $? -ne 0 ]; then
echo "***"
echo "*** ERROR in $0 while running WRAPIT."
echo "***"
exit 1
fi
echo
m=1
me=1 #12
while [ $m -le $me ]; do
$exe $m $dt
m=$(expr $m + 1)
done
exit 0
::::::::::::::
pl.MLD.ILD.BLT.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
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}
echo
echo Done $(basename $0).
echo
exit 0
::::::::::::::
pl.MLD.ILD.BLT.ncl
::::::::::::::
;
; Plot horizontal maps of
; - isothermal-layer depth (ILD)
; - MLD
; - BLT
; 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 ML "./mixed_layer.so"
begin
wallClock1 = systemfunc("export LANG=C; date") ; retrieve wall clock time
;
; default value
;
vlev=2
;
; Read command line arguments
;
scname = getenv("NCL_ARG_1")
infle1 = getenv("NCL_ARG_2")
infle2 = getenv("NCL_ARG_3")
yyyy = getenv("NCL_ARG_4")
m = getenv("NCL_ARG_5")
mmm = getenv("NCL_ARG_6")
dt_in = getenv("NCL_ARG_7")
outdir = getenv("NCL_ARG_8")
dt=stringtofloat(dt_in)
;
; Read input data
;
ft=addfile(infle1,"r")
fs=addfile(infle2,"r")
tglobal3d=ft->t_an
sglobal3d=fs->s_an
;printVarSummary(tglobal3d)
z1d=ft->depth
lat=ft->lat
lon=ft->lon
;printVarSummary(tglobal3d)
t3d=tglobal3d(time|0,depth|:,{lat|-30:30},{lon|30:120})
s3d=sglobal3d(time|0,depth|:,{lat|-30:30},{lon|30:120})
;
; Find ILD
;
dims=dimsizes(t3d)
im=dims(2)
jm=dims(1)
km=dims(0)
ild=t3d(0,:,:)
ild@coordinates:="lat lon"
delete(ild@depth)
ild!0:="lat"
ild!1:="lon"
;printVarSummary(t3d)
;printVarSummary(ild)
;Set depth
klev=vlev-1 ; stringtointeger(vlev)-1
;print(klev)
z=floattointeger(z1d(klev))
;print(z)
t2d = t3d(klev,:,:)
t2d@coordinates:="lat lon"
delete(t2d@depth)
;printVarSummary(t2d)
ML::find_ild(im,jm,km,t3d,t2d,dt,z1d,t3d@_FillValue,ild)
;print(ild)
pd3d=t3d
pd3d@standard_name:="potential density"
delete(pd3d@long_name)
pd3d@coordinates:="depth lat lon"
pd3d@units="kg/m3"
delete(pd3d@time)
ML::pdens3d(im,jm,km,t3d,s3d,z1d,t3d@_FillValue,pd3d)
;printVarSummary(pd3d)
s2d = s3d(klev,:,:)
s2d@coordinates:="lat lon"
delete(s2d@depth)
;printVarSummary(s2d)
ML::pden2d(im,jm,km,t2d,s2d,z1d,t3d@_FillValue,pd2d)
;
; Set year and month
;
dateprt= mmm + " (" + yyyy +")"
;
; Set output file name
;
filetype="png" ;"ps"
mm=sprinti("%0.2i",stringtointeger(m))
ofle=outdir + "/" + "WOA13_MLD.ILD.BLT." + 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 = "dt="+dt_in+"[~S~o~N~C]"
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 = 20 ; set min contour level
res2@cnMaxLevelValF = 40 ; set max contour level
res2@cnLevelSpacingF = 2 ; set contour spacing
res2@gsnLeftString = dateprt ; add the gsn titles
res2@gsnCenterString = "~F33~s~B~q"
res2@gsnRightString = z + " m"
plot2 = gsn_csm_contour_map(wks,pd3d(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
::::::::::::::
wrapit.sh
::::::::::::::
#!/bin/bash
f90=$1
stub=$(basename $f90 .f90).stub
if [ ! -f $f90 ]; then echo Error in $0. No such file, $f90; exit 1; fi
if [ ! -f $fst ]; then echo Error in $0. No such file, $f90; exit 1; fi
export LANG=C
cat <<EOF
Fortran Source file: $f90
Stub file: $stub
EOF
WRAPIT $stub $f90
if [ $? -ne 0 ]; then
echo "***"
echo "*** ERROR in $0 while running WRAPIT."
echo "***"
exit 1
fi
echo
date
echo Shared Libraries:
ls -lh *.so
echo
echo Done $(basename $0).
echo
exit 0
::::::::::::::
mixed_layer.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 find_ild
subroutine pdens3d(im,jm,km,t3d,s3d,z1d,fillvalue,pd3d)
integer,intent(in)::im,jm,km
real,intent(in),dimension(im,jm,km)::t3d,s3d
real,intent(in),dimension(km)::z1d
real,intent(in)::fillvalue
real,intent(inout),dimension(im,jm,km)::pd3d
pd3d=fillvalue
do k=1,km
do j=1,jm
do i=1,im
if(t3d(i,j,k)/=fillvalue .and. s3d(i,j,k)/=fillvalue)then
call sigma_theta(z1d(k),t3d(i,j,k),s3d(i,j,k),pd3d(i,j,k))
end if
end do!i
end do !j
end do !k
end subroutine pdens3d
subroutine sigma_theta(z,t,s,sigth)
implicit none
real,intent(in) :: z,t,s
real,intent(out) :: sigth
! Description:
! sigma_theta = sigma(theta,s,p=0)
! whete, theta is potential temperature, s is salinity, and p is pressure
! Author: aym
real p0 ! pressure in decibars
real pr ! reference pressure in decibars
real th,theta ! potential temperature
real*8 pd,sd,thd,rhod
! Potential temperature
p0=z ! p0 in decibars
pr=0.0
th=theta(s,t,p0,pr)
! Potential density
pd=0.d0 !dble(z)*0.1d0
sd=dble(s)
thd=dble(th)
call eqst2(pd,thd,sd,rhod)
sigth=sngl(rhod-1000.d0)
end subroutine sigma_theta
!***********************************************************************
! eqst2.f90
! TASK
!j 状態方程式(IES80)により海水の密度を計算する。
! REMARK
!
!j All arguments must be DOUBLE PRECISION.
!
!j 入力データは,P [bar], T[deg.], S [psu]
! 圧力の単位がbarであることに注意! (1.0bar = 0.1*dbar)
!
! [USAGE]
!
! Pressure=0.1*z
! call eqst2(pressure,t,s,dens)
!
! where z is depth in meters, t is temperature, s is salinity, and
! dens is density in kg/m3.
!
!
!j 適用範囲
! -2 < T < 40 [deg.]
! 0 < S < 42 [psu]
! 0 < P < 1000 [bar] = 10000 [dbar]
!
!j 水深[m]＝圧力[dbar]という関係を仮定している。これによる両者の誤差は
!j 1%以内(海洋観測指針)。
! It assumes no pressure variation along geopotential surfaces,
! that is, depth and pressure are interchangeable.
!
! REFERENCE
!j 気象庁編():海洋観測指針,428pp.
! UNESCO (1981): Background Papers and Supporting Data on the
! International Equation of State of Seawater. Tech. Pap. Mar.
! Sci., 38, 192pp.
!***********************************************************************
subroutine eqst2(P,T,S,RHO)
implicit none
double precision,intent(in) :: P,T,S
double precision,intent(out):: RHO
double precision b0,b1,b2,b3,b4, c0,c1,c2, d0
double precision a0,a1,a2,a3,a4,a5
double precision f0,f1,f2,f3, g0,g1,g2
double precision i0,i1,i2, j0
double precision m0,m1,m2
double precision e0,e1,e2,e3,e4
double precision h0,h1,h2,h3
double precision k0,k1,k2
parameter(b0= 8.24493d-1, c0=-5.72466d-3 &
& ,b1=-4.0899d-3 , c1= 1.0227d-4 &
& ,b2= 7.6438d-5 , c2=-1.6546d-6 &
& ,b3=-8.2467d-7 &
& ,b4= 5.3875d-9 , d0= 4.8314d-4)
parameter(a0= 999.842594, a1= 6.793952d-2 &
& ,a2=-9.095290d-3,a3= 1.001685d-4 &
& ,a4=-1.120083d-6,a5= 6.536332d-9)
parameter(f0= 54.6746 ,g0= 7.944d-2 &
& ,f1= -0.603459 ,g1= 1.6483d-2 &
& ,f2= 1.09987d-2,g2=-5.3009d-4 &
& ,f3= -6.1670d-5)
parameter(i0= 2.2838d-3 ,j0= 1.91075d-4 &
& ,i1=-1.0981d-5 &
& ,i2=-1.6078d-6)
parameter(m0=-9.9348d-7, m1= 2.0816d-8 &
& ,m2= 9.1697d-10)
parameter(e0= 19652.21, e1= 148.4206 &
& ,e2=-2.327105, e3= 1.360477d-2 &
& ,e4=-5.155288e-5)
parameter(h0= 3.239908, h1= 1.43713d-3 &
& ,h2= 1.16092d-4, h3=-5.77905d-7)
parameter(k0= 8.50935d-5, k1=-6.12293d-6 &
& ,k2= 5.2787d-8)
integer k
double precision rhow,rhop0,Kw,Aw,Bw,Kp0,A,B,Kp
!-----------------------------------------------------------------------
!j 以下の式番号はすべて海洋観測指針(p.92～93)に記載されている式番号に
!j 対応する。
! Eq.(4)
rhow=a0 + a1*t + a2*t**2 + a3*t**3 + a4*t**4 &
& + a5*t**5
! Eq.(3)
rhop0=rhow + (b0 + b1*t + b2*t**2 + b3*t**3 &
& + b4*t**4)*s &
& + (c0 + c1*t + c2*t**2)*s**(1.5) &
& + d0*s**2
! Eq.(9)
Kw=e0 + e1*t + e2*t**2 + e3*t**3 + e4*t**4
! Eq.(10)
Aw=h0 + h1*t + h2*t**2 + h3*t**3
! Eq.(11)
Bw=k0 + k1*t + k2*t**2
! Eq.(6)
Kp0=Kw + (f0 + f1*t + f2*t**2 + f3*t**3)*s &
& + (g0 + g1*t + g2*t**2)*S**(1.5)
! Eq.(7)
A=Aw + (i0 + i1*t + i2*t**2)*s + j0*S**(1.5)
! Eq.(8)
B=Bw + (m0 + m1*t + m2*t**2)*s
! Eq.(5)
Kp=Kp0 + A*p + B*p**2
! Eq.(1)
rho=rhop0/(1.0d0 - p/Kp)
end subroutine eqst2
function theta(s,t0,p0,pr)
!
! to compute local potential temperature at pr
! using bryden 1973 polynomial for adiabatic lapse rate
! and runge-kutta 4-th order integration algorithm.
! ref: bryden,h.,1973,deep-sea res.,20,401-408
! fofonoff,n.,1977,deep-sea res.,24,489-491
! units:
! pressure p0 *** decibars ***
! temperature t0 deg celsius (ipts-68)
! salinity s (ipss-78)
! reference prs pr decibars
! potential tmp. theta deg celsius
! checkvalue: theta= 36.89073 c,s=40 (ipss-78),t0=40 deg c,
! p0=10000 decibars,pr=0 decibars
!
! set-up intermediate temperature and pressure variables
!
! http://www.marine.csiro.au/~jackett/NeutralDensity/
implicit none
real,intent(in)::s,t0,p0,pr
real p,t,h,xk,q,theta,atg
p=p0
t=t0
h = pr - p
xk = h*atg(s,t,p)
t = t + 0.5*xk
q = xk
p = p + 0.5*h
xk = h*atg(s,t,p)
t = t + 0.29289322*(xk-q)
q = 0.58578644*xk + 0.121320344*q
xk = h*atg(s,t,p)
t = t + 1.707106781*(xk-q)
q = 3.414213562*xk - 4.121320344*q
p = p + 0.5*h
xk = h*atg(s,t,p)
theta = t + (xk-2.0*q)/6.0
end function theta
function atg(s,t,p)
!
! adiabatic temperature gradient deg c per decibar
! ref: bryden,h.,1973,deep-sea res.,20,401-408
! units:
! pressure p decibars
! temperature t deg celsius (ipts-68)
! salinity s (ipss-78)
! adiabatic atg deg. c/decibar
! checkvalue: atg=3.255976e-4 c/dbar for s=40 (ipss-78),
! t=40 deg c,p0=10000 decibars
!
! http://www.marine.csiro.au/~jackett/NeutralDensity/
implicit none
real,intent(in)::s,t,p
real atg,ds
ds = s - 35.0
atg = (((-2.1687d-16*t+1.8676d-14)*t-4.6206d-13)*p &
+((2.7759d-12*t-1.1351d-10)*ds+((-5.4481d-14*t &
+8.733d-12)*t-6.7795d-10)*t+1.8741d-8))*p &
+(-4.2393d-8*t+1.8932d-6)*ds &
+((6.6228d-10*t-6.836d-8)*t+8.5258d-6)*t+3.5803d-5
end function atg
::::::::::::::
mixed_layer.stub
::::::::::::::
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
C NCLFORTSTART
subroutine pdens3d(im,jm,km,t3d,s3d,z1d,fillvalue,pd3d)
integer im,jm,km
real t3d(im,jm,km),s3d(im,jm,km)
real z1d(km)
real fillvalue
real pd3d(im,jm,km)
C NCLEND
C NCLFORTSTART
subroutine sigma_theta(z,t,s,sigth)
real z,t,s
real sigth
C NCLEND
C NCLFORTSTART
subroutine eqst2(P,T,S,RHO)
double precision P,T,S
double precision RHO
C NCLEND
C NCLFORTSTART
function theta(s,t0,p0,pr)
real s,t0,p0,pr
C NCLEND
C NCLFORTSTART
function atg(s,t,p)
real s,t,p
C NCLEND