::::::::::::::
pl.deBoyer.MLD.run.sh
::::::::::::::
#!/bin/sh
export LANG=C
#
# Creating shared library using NCL WRAPIT utility
#
fort=mixed_layer.f90
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
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
yyyys=1993
yyyye=2014
yearday=n10day_01_09
dthr=240
day=01
dt=0.2
indir=/work4/data/ECCO.dr080g
outdir=/work4/data/ECCO.dr080g
if [ ! -d $indir ]; then
echo
echo Error in $0 : No such directory, $indir
echo
exit 1
fi
figdir=${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 \
"
yyyy=$yyyys
while [ $yyyy -le $yyyye ] ; do
stime=$(date)
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
echo ${yyyy} ${yearday}
$exe $ncl ${indir} ${yyyy} ${yearday} ${datainfo} ${dt} ${outdir} ${figdir}
done #yearday
yyyy=$(expr $yyyy + 1)
etime=$(date)
echo "$yyyy"
echo "Started at $stime"
echo "Finished at $etime"
echo
done #yyyy
echo
echo Done $(basename $0).
echo
exit 0
::::::::::::::
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
::::::::::::::
!
! TTD
! http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Warm_Layer_Depth.php#WLD_DTm02
! MLD
! http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Mixed_Layer_Depth.php#MLD_DReqDTm02
! BLT
! http://www.ifremer.fr/cerweb/deboyer/mld/Subsurface_Barrier_Layer_Thickness.php
!
subroutine find_ttd(nx,ny,nz,t3d,t2d,dt,z1d,fillvalue,ttd)
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)::ttd(nx,ny)
! z1d(k): Model's vertical coordinate
! ttd=Top of Thermocline Depth
do j=1,ny
do i=1,nx
ttd(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
ttd(i,j)=ys
end if
end do !k
! Otherwise, the water column should be vertically homogeneous.
! ttd(i,j)=fillvalue
end do !i
end do !j
end subroutine find_ttd
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 pdens2d(im,jm,t2d,s2d,z,fillvalue,pd2d)
integer,intent(in)::im,jm
real,intent(in),dimension(im,jm)::t2d,s2d
real,intent(in)::z
real,intent(in)::fillvalue
real,intent(inout),dimension(im,jm)::pd2d
pd2d=fillvalue
do j=1,jm
do i=1,im
if(t2d(i,j)/=fillvalue .and. s2d(i,j)/=fillvalue)then
call sigma_theta(z,t2d(i,j),s2d(i,j),pd2d(i,j))
end if
end do!i
end do !j
end subroutine pdens2d
subroutine dsigth(im,jm,pd1,pd2,fillvalue,dpd)
integer,intent(in)::im,jm
real,intent(in),dimension(im,jm)::pd1, pd2
real,intent(in)::fillvalue
real,intent(inout),dimension(im,jm)::dpd
dpd=fillvalue
do j=1,jm
do i=1,im
if(pd1(i,j)/=fillvalue .and. pd2(i,j)/=fillvalue)then
dpd(i,j)=pd1(i,j)-pd2(i,j)
end if
end do!i
end do !j
end subroutine dsigth
subroutine find_mld(nx,ny,nz,pd3d,pd2d,dpd,z1d,fillvalue,mld)
integer,intent(in)::nx,ny,nz
real,intent(in)::pd3d(nx,ny,nz),pd2d(nx,ny),dpd(nx,ny)
real,intent(in)::z1d(nz)
real,intent(in)::fillvalue
real,intent(inout)::mld(nx,ny)
! z1d(k): Model's vertical coordinate
! ttd=Top of Thermocline Depth
do j=1,ny
do i=1,nx
mld(i,j)=fillvalue ! default = undefined
xs= pd2d(i,j) + dpd(i,j)
do k=1,nz-1
if(pd3d(i,j,k) <= xs .and. xs < pd3d(i,j,k+1) )then
y1=z1d(k)
y2=z1d(k+1)
x1=pd3d(i,j,k)
x2=pd3d(i,j,k+1)
!Linear Interpolation
ys=(y2-y1)/(x2-x1)*(xs-x1) + y1
mld(i,j)=ys
end if
end do !k
end do !i
end do !j
end subroutine find_mld
subroutine calc_blt(nx,ny,ttd,mld,fillvalue,blt)
integer,intent(in)::nx,ny
real,intent(in)::ttd(nx,ny),mld(nx,ny)
real,intent(in)::fillvalue
real,intent(inout)::blt(nx,ny)
real,parameter::eps=0.1
! z1d(k): Model's vertical coordinate
! ttd=Top of Thermocline Depth
do j=1,ny
do i=1,nx
blt(i,j)=fillvalue
if(mld(i,j)/=fillvalue .and. ttd(i,j)/=fillvalue )then
if(mld(i,j)>eps .and. ttd(i,j)>eps )then
blt(i,j)=ttd(i,j)-mld(i,j)
end if
end if
end do !i
end do !j
end subroutine calc_blt
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_ttd(nx,ny,nz,t3d,t2d,dt,z1d,fillvalue,ttd)
integer nx,ny,nz
real t3d(nx,ny,nz)
real t2d(nx,ny)
real dt
real z1d(nz)
real fillvalue
real ttd(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 pdens2d(im,jm,t2d,s2d,z,fillvalue,pd2d)
integer im,jm
real t2d(im,jm),s2d(im,jm)
real z
real fillvalue
real pd2d(im,jm)
C NCLEND
C NCLFORTSTART
subroutine dsigth(im,jm,pd1,pd2,fillvalue,dpd)
integer im,jm
real pd1(im,jm), pd2(im,jm)
real fillvalue
real dpd(im,jm)
C NCLEND
C NCLFORTSTART
subroutine find_mld(nx,ny,nz,pd3d,pd2d,dpd,z1d,fillvalue,mld)
integer nx,ny,nz
real pd3d(nx,ny,nz),pd2d(nx,ny),dpd(nx,ny)
real z1d(nz)
real fillvalue
real mld(nx,ny)
C NCLEND
C NCLFORTSTART
subroutine calc_blt(nx,ny,ild,mld,fillvalue,blt)
integer nx,ny
real ild(nx,ny),mld(nx,ny)
real fillvalue
real blt(nx,ny)
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
::::::::::::::
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.deBoyer.MLD.ncl
::::::::::::::
;
; Definitions of variables to be computed
;
; TTD
; http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Warm_Layer_Depth.php#WLD_DTm02
; MLD
; http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Mixed_Layer_Depth.php#MLD_DReqDTm02
; BLT
; http://www.ifremer.fr/cerweb/deboyer/mld/Subsurface_Barrier_Layer_Thickness.php
;
; Dataset
; 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"
external ML "./mixed_layer.so"
begin
wallClock1 = systemfunc("export LANG=C; date") ; retrieve wall clock time
; Default values
dataset="dr080g"
dthr=240 ;
zref=10.0 ;Reference depth for TTD and MLD [m]
zref0=0.0 ;Reference depth for P=0 [db]
dt=0.2 ;delta T for TTD and MLD [degC]
figtype="png" ;"ps"
;
; Read command line arguments
;
scname = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
iy = getenv("NCL_ARG_3")
yearday = getenv("NCL_ARG_4")
datainfo = getenv("NCL_ARG_5")
dt = stringtofloat( getenv("NCL_ARG_6"))
outdir = getenv("NCL_ARG_7")
figdir = getenv("NCL_ARG_8")
;
; Read input data
;
indir1=indir+"/"+dataset+"_"+iy+"/"+yearday+"/"
fT =addfile(indir1 + "Tave" + datainfo + dthr + ".cdf" ,"r")
T = fT->Tave(:,:,:,:)
fS =addfile(indir1 + "Save" + datainfo + dthr + ".cdf" ,"r")
S = fS->Save(:,:,:,:)
lat =fT->lat
lat_v=fT->lat_v
lon =fT->lon
lon_u=fT->lon_u
depth=fT->depth
time=fT->time
nt=dimsizes(time)
dim=dimsizes(T)
nz=dim(1)
ny=dim(2)
nx=dim(3)
t3d=new((/nz,ny,nx/),float)
t3d@_FillValue=T@_FillValue
s3d=new((/nz,ny,nx/),float)
s3d@_FillValue=T@_FillValue
t2d=new((/ny,nx/),float)
t2d@_FillValue=T@_FillValue
s2d=new((/ny,nx/),float)
s2d@_FillValue=T@_FillValue
pd2d_r=new((/ny,nx/),float)
pd2d_r@_FillValue=T@_FillValue
pd2d_b=new((/ny,nx/),float)
pd2d_b@_FillValue=T@_FillValue
dpd2d=new((/ny,nx/),float)
dpd2d@_FillValue=T@_FillValue
pd3d=new((/nz,ny,nx/),float)
pd3d@_FillValue=T@_FillValue
ttd_tmp=new((/ny,nx/),float)
ttd_tmp@_FillValue=T@_FillValue
ttd_tmp!0="lat"
ttd_tmp!1="lon"
ttd_tmp&lon=lon
ttd_tmp&lat=lat
mld_tmp=new((/ny,nx/),float)
mld_tmp@_FillValue=T@_FillValue
mld_tmp!0="lat"
mld_tmp!1="lon"
mld_tmp&lon=lon
mld_tmp&lat=lat
blt_tmp=new((/ny,nx/),float)
blt_tmp@_FillValue=T@_FillValue
blt_tmp!0="lat"
blt_tmp!1="lon"
blt_tmp&lon=lon
blt_tmp&lat=lat
ttd=new((/nt,ny,nx/),float)
ttd@_FillValue=T@_FillValue
mld=new((/nt,ny,nx/),float)
mld@_FillValue=T@_FillValue
blt=new((/nt,ny,nx/),float)
blt@_FillValue=T@_FillValue
ttd@long_name = "Top of thermocline depth"
ttd@data_source = "http://ecco.jpl.nasa.gov/thredds/fileServer/ecco_Smoother"
ttd@units = "m"
ttd@zref = zref + " [m]"
ttd@dt = dt + " [degC]"
ttd!0="time"
ttd!1="lat"
ttd!2="lon"
ttd&lon=lon
ttd&lat=lat
mld@long_name = "Mixed layer depth"
mld@data_source = "http://ecco.jpl.nasa.gov/thredds/fileServer/ecco_Smoother"
mld@units = "m"
mld@zref = zref + " [m]"
mld@dt = dt + " [degC]"
mld!0="time"
mld!1="lat"
mld!2="lon"
mld&lon=lon
mld&lat=lat
blt@long_name = "Barrier layer thickness"
blt@data_source = "http://ecco.jpl.nasa.gov/thredds/fileServer/ecco_Smoother"
blt@units = "m"
blt@zref = zref + " [m]"
blt@dt = dt + " [degC]"
blt!0="time"
blt!1="lat"
blt!2="lon"
blt&lon=lon
blt&lat=lat
zref_check=0.5*(depth(0)+depth(1))
if(zref_check .ne. zref)then
print("Error : Check zref = "+ zref_check + " zref = " + zref)
exit
end if
;
; Set NCL resources
;
res = True
res@gsnDraw = False ; don't draw
res@gsnFrame = False ; don't advance frame
;Figure size
res@vpHeightF = 0.2 ;0.35
res@vpWidthF = 0.25 ;0.55
res@vpYF = 0.6 ;0.2
res@vpXF = 0.1 ;0.2
;Map
res@mpFillOn = True
res@gsnPaperOrientation = "landscape"
res@mpShapeMode = "FreeAspect"
res@gsnAddCyclic = False
;Contour
res@cnFillOn = True ; turn on color fill
res@cnLineColor = "white"
res@lbOrientation = "horizontal"
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
do n=0,nt-1
t3d(:,:,:)=T(n,:,:,:)
s3d(:,:,:)=S(n,:,:,:)
;
; Reference temperature and salinity
;
t2d(:,:)=(t3d(0,:,:)+t3d(1,:,:))*0.5 ;T at zref
t2d_b=t2d - dt
s2d(:,:)=(s3d(0,:,:)+s3d(1,:,:))*0.5 ;S at zref
;
; Compute TTD
;
ML::find_ttd(nx,ny,nz,t3d,t2d,dt,depth,T@_FillValue,ttd_tmp)
;
; Compute MLD
;
ML::pdens2d(nx,ny,t2d,s2d,zref0,T@_FillValue,pd2d_r) ;Pot. Dens (PD)
;at 10m
ML::pdens2d(nx,ny,t2d_b,s2d,zref0,T@_FillValue,pd2d_b)
ML::dsigth(nx,ny,pd2d_b,pd2d_r,T@_FillValue,dpd2d) ;PD diff.
ML::pdens3d(nx,ny,nz,t3d,s3d,depth,T@_FillValue,pd3d)
ML::find_mld(nx,ny,nz,pd3d,pd2d_r,dpd2d,depth,T@_FillValue,mld_tmp)
;
; Compute BLT
;
ML::calc_blt(nx,ny,ttd_tmp,mld_tmp,T@_FillValue,blt_tmp)
; For output
ttd(n,:,:) = ttd_tmp(:,:)
mld(n,:,:) = mld_tmp(:,:)
blt(n,:,:) = blt_tmp(:,:)
;
; 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(n)/24.0 -5
julim1=julian0 + time(n)/24.0 +4
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
;
; Plotting
;
figdir1=figdir+"/"+dataset+"_"+iy+"/"+yearday
foo=systemfunc("mkdir -vp "+figdir1)
mm=sprinti("%0.2i",stringtointeger(n+1))
figfle=figdir1+"/"+"mld" + datainfo + dthr + "."+ mm +"."+figtype
print("Figure: "+figfle)
wks = gsn_open_wks(figtype,figfle)
;
; Left panel
;
res1=res
;res1@vpYF = 0.9
res1@vpXF = 0.05
res1@cnMinLevelValF = 10 ; set min contour level
res1@cnMaxLevelValF = 100 ; set max contour level
res1@cnLevelSpacingF = 10 ; set contour spacing
res1@gsnLeftString = "" ; add the gsn titles
res1@gsnCenterString = "MLD [m]"
res1@gsnRightString = ""
plot1 = gsn_csm_contour_map(wks,mld_tmp(:,:),res1) ; Draw a contour plot.
; now assign plotA to array
draw(plot1)
delete(res1)
;
; Center panel
;
res2=res
;res2@vpYF = 0.45
res2@vpXF = 0.35
res2@cnMinLevelValF = 10 ; set min contour level
res2@cnMaxLevelValF = 100 ; set max contour level
res2@cnLevelSpacingF = 10 ; set contour spacing
res2@gsnLeftString = "" ; add the gsn titles
res2@gsnCenterString = "TTD [m]"
res2@gsnRightString = ""
plot2 = gsn_csm_contour_map(wks,ttd_tmp(:,:),res2) ; Draw a contour plot.
draw(plot2)
delete(res2)
;
; Right panel
;
gsn_define_colormap(wks,"precip3_16lev")
res3=res
;res3@vpYF = 0.45
res3@vpXF = 0.65
res3@cnMinLevelValF = 0 ; set min contour level
res3@cnMaxLevelValF = 40 ; set max contour level
res3@cnLevelSpacingF = 5 ; set contour spacing
res3@gsnLeftString = "" ; add the gsn titles
res3@gsnCenterString = "BLT [m]"
res3@gsnRightString = ""
plot3 = gsn_csm_contour_map(wks,blt_tmp(:,:),res3) ; Draw a contour plot.
draw(plot3)
delete(res3)
txres = True ; text mods desired
txres@txBackgroundFillColor = "white"
txres@txFontColor="black"
txres@txFontHeightF = 0.02
text=dataset
gsn_text_ndc(wks,text,0.475,0.8,txres)
txres@txFontHeightF = 0.015
text="TTD, MLD, and BLT using de Boyer's definition"
gsn_text_ndc(wks,text,0.475,0.75,txres)
text="Z~B~ref~N~="+zref+" [m] ~F33~D~F21~T="+dt+" [~S~o~N~C]"
gsn_text_ndc(wks,text,0.475,0.7,txres)
text=dateprt
gsn_text_ndc(wks,text,0.475,0.65,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)
end do ; n
;
; Output netCDF files
;
outdir1=outdir+"/"+dataset+"_"+iy+"/"+yearday
ofle1=outdir1+"/"+"ttd" + datainfo + dthr + ".cdf"
ofle2=outdir1+"/"+"mld" + datainfo + dthr + ".cdf"
ofle3=outdir1+"/"+"blt" + datainfo + dthr + ".cdf"
system("/bin/rm -f " + ofle1)
fout1 =addfile (ofle1, "c")
system("/bin/rm -f " + ofle2)
fout2 =addfile (ofle2, "c")
system("/bin/rm -f " + ofle3)
fout3 =addfile (ofle3, "c")
setfileoption(fout1,"DefineMode",True)
setfileoption(fout1,"MissingToFillValue",True)
setfileoption(fout2,"DefineMode",True)
setfileoption(fout2,"MissingToFillValue",True)
setfileoption(fout3,"DefineMode",True)
setfileoption(fout3,"MissingToFillValue",True)
fAtt = True ; assign file attributes
fAtt@title = "ECCO dr080g (Monthly Climatology)"
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
fileattdef( fout3, fAtt ) ; copy file attributes
ntim = dimsizes(time) ; get dimension sizes
nlat = dimsizes(lat)
nlon = dimsizes(lon)
dimNames = (/"time", "lat", "lon"/)
dimSizes = (/ ntim , nlat, nlon/)
dimUnlim = (/True , False, False/)
filedimdef(fout1,dimNames,dimSizes,dimUnlim)
filedimdef(fout2,dimNames,dimSizes,dimUnlim)
filedimdef(fout3,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, "ttd" ,typeof(ttd), getvardims(ttd))
filevardef(fout2, "time",typeof(time),getvardims(time))
filevardef(fout2, "lat" ,typeof(lat), getvardims(lat))
filevardef(fout2, "lon" ,typeof(lon), getvardims(lon))
filevardef(fout2, "mld" ,typeof(mld), getvardims(mld))
filevardef(fout3, "time",typeof(time),getvardims(time))
filevardef(fout3, "lat" ,typeof(lat), getvardims(lat))
filevardef(fout3, "lon" ,typeof(lon), getvardims(lon))
filevardef(fout3, "blt" ,typeof(blt), getvardims(blt))
filevarattdef(fout1,"time" ,time)
filevarattdef(fout1,"lat" ,lat )
filevarattdef(fout1,"lon" ,lon )
filevarattdef(fout1,"ttd" ,ttd )
filevarattdef(fout2,"time" ,time)
filevarattdef(fout2,"lat" ,lat )
filevarattdef(fout2,"lon" ,lon )
filevarattdef(fout2,"mld" ,mld )
filevarattdef(fout3,"time" ,time)
filevarattdef(fout3,"lat" ,lat )
filevarattdef(fout3,"lon" ,lon )
filevarattdef(fout3,"blt" ,blt )
fout1->time=(/time/)
fout1->lat =(/lat/)
fout1->lon =(/lon/)
fout1->ttd =(/ttd/)
fout2->time=(/time/)
fout2->lat =(/lat/)
fout2->lon =(/lon/)
fout2->mld =(/mld/)
fout3->time=(/time/)
fout3->lat =(/lat/)
fout3->lon =(/lon/)
fout3->blt =(/blt/)
print("Output file:"+ofle1)
print("Output file:"+ofle2)
print("Output file:"+ofle3)
print("")
end