MLD was estimated using de Boyer's definition.
[2015年 9月 18日 金曜日 12:21:46 JST]
[~/2015.Ibnu.Indian.Ocean/ECCO.dr080g/Map.MLTS.deBoyer]
[am@aofd165]
$ pl.deBoyer.MLTS.run.sh
::::::::::::::
pl.deBoyer.MLTS.ncl
::::::::::::::
;
;
; Main NCL script
;
;
;
; 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 "./mlts.so"
begin
wallClock1 = systemfunc("export LANG=C; date") ; retrieve wall clock time
; Default values
dataset="dr080g"
dthr=240 ;
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")
outdir = getenv("NCL_ARG_6")
figdir = getenv("NCL_ARG_7")
;
; Read input data
;
indir1=indir+"/"+dataset+"_"+iy+"/"+yearday+"/"
fT =addfile(indir1 + "Tave" + datainfo + dthr + ".cdf" ,"r")
fS =addfile(indir1 + "Save" + datainfo + dthr + ".cdf" ,"r")
fD =addfile(indir1 + "mld" + datainfo + dthr + ".cdf" ,"r")
T = fT->Tave(:,:,:,:)
S = fS->Save(:,:,:,:)
mld = fD->mld(:,:,:)
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)
zref=fD->mld@zref
dt=fD->mld@dt
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
mld2d=new((/ny,nx/),float)
mld2d@_FillValue=T@_FillValue
mlt=new((/nt,ny,nx/),float)
mlt@_FillValue=T@_FillValue
mls=new((/nt,ny,nx/),float)
mls@_FillValue=T@_FillValue
mlt@long_name = "Mixed layer temperature"
mlt@data_source = "http://ecco.jpl.nasa.gov/thredds/fileServer/ecco_Smoother"
mlt@units = "degC"
mlt@zref = zref
mlt@dt = dt
mlt!0="time"
mlt!1="lat"
mlt!2="lon"
mlt&lon=lon
mlt&lat=lat
mls@long_name = "Mixed layer salinity"
mls@data_source = "http://ecco.jpl.nasa.gov/thredds/fileServer/ecco_Smoother"
mls@units = "PSU"
mls@zref = zref
mls@dt = dt
mls!0="time"
mls!1="lat"
mls!2="lon"
mls&lon=lon
mls&lat=lat
mlt_tmp=new((/ny,nx/),float)
mlt_tmp@_FillValue=T@_FillValue
mlt_tmp!0="lat"
mlt_tmp!1="lon"
mlt_tmp&lon=lon
mlt_tmp&lat=lat
mls_tmp=new((/ny,nx/),float)
mls_tmp@_FillValue=T@_FillValue
mls_tmp!0="lat"
mls_tmp!1="lon"
mls_tmp&lon=lon
mls_tmp&lat=lat
;
; 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@gsnPaperOrientation = "landscape"
res@mpShapeMode = "FreeAspect"
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
kk_out=new( (/ny,nx/),integer)
do n=0,nt-1
t3d(:,:,:)=T(n,:,:,:)
s3d(:,:,:)=S(n,:,:,:)
mld2d(:,:)=mld(n,:,:)
;
; Vertical average in the mixed layer
;
ML::mlave(nx,ny,nz,t3d,mld2d,depth,T@_FillValue,mlt_tmp)
ML::mlave(nx,ny,nz,s3d,mld2d,depth,T@_FillValue,mls_tmp)
; For output
mlt(n,:,:) = mlt_tmp(:,:)
mls(n,:,:) = mls_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+"/"+"mlts" + datainfo + dthr + "."+ mm +"."+figtype
print("Figure: "+figfle)
wks = gsn_open_wks(figtype,figfle)
;
; 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@gsnLeftString = "MLT"
res1@gsnCenterString = dateprt
tcon = gsn_csm_contour_map(wks,mlt(n,:,:),res1) ; Draw a contour plot.
draw(tcon)
;
; 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 = "MLS"
res2@gsnCenterString = dateprt
scon = gsn_csm_contour_map(wks,mls(n,:,:),res2) ; Draw a contour plot.
draw(scon)
txres = True ; text mods desired
txres@txBackgroundFillColor = "white"
txres@txFontColor="black"
txres@txFontHeightF = 0.012
text=dataset + " MLT and MLS (de Boyer)"
gsn_text_ndc(wks,text,0.475,0.98,txres)
text="Z~B~ref~N~="+zref+" ~F33~D~F21~T="+dt
gsn_text_ndc(wks,text,0.475,0.96,txres)
;
; Header and footer
;
drawNDCGrid(wks)
txres = True ; text mods desired
txres@txFontHeightF = 0.01 ; font smaller. default big
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)
gsn_text_ndc(wks,outdir,0.8,0.03,txres)
frame(wks)
end do ; n
;
; Output netCDF files
;
outdir1=outdir+"/"+dataset+"_"+iy+"/"+yearday
ofle1=outdir1+"/"+"mlt" + datainfo + dthr + ".cdf"
ofle2=outdir1+"/"+"mls" + datainfo + dthr + ".cdf"
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"
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 = 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)
filevardef(fout1, "time",typeof(time),getvardims(time))
filevardef(fout1, "lat" ,typeof(lat), getvardims(lat))
filevardef(fout1, "lon" ,typeof(lon), getvardims(lon))
filevardef(fout1, "mlt" ,typeof(mlt), getvardims(mlt))
filevardef(fout2, "time",typeof(time),getvardims(time))
filevardef(fout2, "lat" ,typeof(lat), getvardims(lat))
filevardef(fout2, "lon" ,typeof(lon), getvardims(lon))
filevardef(fout2, "mls" ,typeof(mls), getvardims(mls))
filevarattdef(fout1,"time" ,time)
filevarattdef(fout1,"lat" ,lat )
filevarattdef(fout1,"lon" ,lon )
filevarattdef(fout1,"mlt" ,mlt )
filevarattdef(fout2,"time" ,time)
filevarattdef(fout2,"lat" ,lat )
filevarattdef(fout2,"lon" ,lon )
filevarattdef(fout2,"mls" ,mls )
fout1->time=(/time/)
fout1->lat =(/lat/)
fout1->lon =(/lon/)
fout1->mlt =(/mlt/)
fout2->time=(/time/)
fout2->lat =(/lat/)
fout2->lon =(/lon/)
fout2->mls =(/mls/)
print("Output file:"+ofle1)
print("Output file:"+ofle2)
print("")
end
::::::::::::::
pl.deBoyer.MLTS.run.sh
::::::::::::::
#!/bin/sh
#
#
# Driver for NCL script
#
#
export LANG=C
#
# Creating shared library using NCL WRAPIT utility
#
fort=mlts.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
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} ${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
#
#
# Driver for NCL WRAPIT utility
#
#
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
::::::::::::::
mlts.f90
::::::::::::::
!
!
! Fortran subroutines for NCL script
!
!
subroutine mlave(nx,ny,nz,v3d,mld,depth,fillvalue,mla)
integer,intent(in)::nx,ny,nz
real,intent(in)::v3d(nx,ny,nz)
real,intent(in)::mld(nx,ny)
real,intent(in)::depth(nz)
real,intent(in)::fillvalue
real,intent(out)::mla(nx,ny)
mla(:,:)=fillvalue ! default = undefined
do j=1,ny
do i=1,nx
if(mld(i,j)==fillvalue)then
cycle
end if
kk=0
do k=1,nz-1
if(depth(k) > mld(i,j))goto 100
kk=k
end do !k
100 continue
if(kk == 0)then
mla(i,j)=fillvalue
cycle
endif
if(kk <= 2)then
mla(i,j)=v3d(i,j,kk)
cycle
endif
tsum=0.0
do k=1,kk-1
x1=depth(k)
x2=depth(k+1)
y1=v3d(i,j,k)
y2=v3d(i,j,k+1)
tsum=tsum + (y2+y1)*(x2-x1)/2.0
enddo !k
x1=depth(kk)
x2=depth(kk+1)
y1=v3d(i,j,kk)
y2=v3d(i,j,kk+1)
v3d_at_zmld=(y2-y1)/(x2-x1)*(mld(i,j)-x1) + y1
tsum2=0.0
x1=depth(kk)
x2=mld(i,j)
y1=v3d(i,j,kk)
y2=v3d_at_zmld
tsum2=(y2+y1)*(x2-x1)/2.0
mla(i,j)=(tsum + tsum2)/(mld(i,j)-depth(1))
end do !i
end do !j
end subroutine mlave
::::::::::::::
mlts.stub
::::::::::::::
C
C
C Interface for Fortran 90 subroutines
C
C
C NCLFORTSTART
subroutine mlave(nx,ny,nz,v3d,mld,depth,fillvalue,mla)
integer nx,ny,nz
real v3d(nx,ny,nz)
real mld(nx,ny)
real depth(nz)
real fillvalue
real mla(nx,ny)
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
CYCLE文
CYCLE文は「ループの最後まで飛ぶ文」である。
論理IFと組み合わせて用いる文で、ループの最後、つまり端末文(END DOやCONTINUE)に向かって飛ぶ。
言い換えれば、CYCLE文から端末文(END DO)までの処理を行なわずに、ループの先頭に戻る、と言える。
CYCLE文の例：1から10までの偶数の総和を求めるプログラム
program sum_even
integer i,sum
sum=0
do i=1,10
if(mod(i,2) /= 0)cycle
sum=sum+i
end do
print *,'Sum=',sum
stop
end
if(mod(i,2) /= 0)cycle
上の文の意味:
2で割った余りが0でない、つまり奇数なら足し算の処理は行なわず、ループの先頭へと戻る。
偶数であればこの条件に引っかからないので、下へと進む。
注意：
このプログラムはあくまでcycle文の使用例であって, DO文のところを「do i=2,100,2」と書いてしまえばCYCLE文がなくても偶数の総和を求めることができる。