$ jcope2sst.run.sh
$ ls output
jcope2.sst.2011-06-25_00.txt
$ cd plot
$ pl.sst.3.sh ../output/jcope2.sst.2011-06-25_00.txt
Bash script ./pl.sst.3.sh starts.
Input : ../output/jcope2.sst.2011-06-25_00.txt
Output : ../ps/jcope2.sst.2011-06-25_00.2.ps
pstext: Record 0 is incomplete (skipped)
pstext: Record 6 is incomplete (skipped)
::::::::::::::
jcope2sst.run.sh
::::::::::::::
#!/bin/sh
exe=jcope2sst.exe
# time info
syyyy=2012
smm=07
sdd=11
days=1
dh=6 # Output interval in hours
para=para.$(basename $0 .sh).txt
# JCOPE2 data info
indir=/work4/DATA/JCOPE2/201207/data #/work4/DATA/JCOPE2/201305-08
#indir=input
xlons=108.
xlons_offset_denom=24.0
ylats=10.5
ylats_offset_denom=24.0
im=866 #im=626 # 2000-2009
jm=620 #jm=476 # 2000-2009
km=47
tbias=10.0
deltalat=0.0833333
deltalon=0.0833333
outdir=./output
if [ ! -d $outdir ]; then
mkdir -vp $outdir
fi
outdir_check=./out.check
if [ ! -d $outdir_check ]; then
mkdir -vp $outdir_check
fi
d=0
while [ $d -le $days ]; do
d1=$d
d2=$(expr $d1 + 1)
yyyymmdd1=$(date -d"${syyyy}${smm}${sdd} $d1 day" '+%Y%m%d')
yyyymmdd2=$(date -d"${syyyy}${smm}${sdd} $d2 day" '+%Y%m%d')
yyyy1=${yyyymmdd1:0:4}
mm1=${yyyymmdd1:4:2}
dd1=${yyyymmdd1:6:2}
yyyy2=${yyyymmdd2:0:4}
h=0
while [ $h -lt 24 ]; do
ifile1=${indir}/T_${yyyymmdd1}
ifile2=${indir}/T_${yyyymmdd2}
hh=$(printf "%02d" $h)
ofile="${outdir}/jcope2.sst.${yyyy1}-${mm1}-${dd1}_${hh}.txt"
hdate="${yyyy1}-${mm1}-${dd1}-${hh}:00:00"
cxfcst=${hh}
now=$(date)
cwd=$(pwd)
host=$(hostname)
user=$(whoami)
cat <<EOF > $para
!
! namelist for ${exe}
!
! $host
! $cwd
! $user
! $now
¶
ifile1="${ifile1}",
ifile2="${ifile2}",
ofile="${ofile}",
hdate="${hdate}",
cxfcst="${cxfcst}",
yyyy=$yyyy1,
mm=$mm1,
dd=$dd1,
hh=$hh,
im=${im},
jm=${jm},
km=${km},
tbias = ${tbias}
xlons=${xlons}
xlons_offset_denom=${xlons_offset_denom},
ylats=${ylats},
ylats_offset_denom=${ylats_offset_denom},
deltalat=${deltalat},
deltalon=${deltalon},
&end
EOF
h=$( expr $h + $dh )
$exe < $para
done # h
d=$(expr $d + 1)
done #d
exit 0
::::::::::::::
jcope2sst.f90
::::::::::::::
program jcope2sst
! constant output variables (need configure)
integer, parameter :: version = 5 ! Format version (must =5 for WPS format)
integer, parameter :: iproj = 0 ! Code for projection of data in array:
! 0 = cylindrical equidistant (either global or regional)
! 1 = Mercator (regional only)
! 3 = Lambert conformal conic (regional only)
! 4 = Gaussian (global only!)
! 5 = Polar stereographic (regional only)
real, parameter :: nlats = 0. ! Number of latitudes north of equator (for Gaussian grids)
real, parameter :: xlvl = 200100. ! Vertical level of data in 2-d array (200100 for surface, 201300 for SLP)
real, parameter :: dx = 0. ! x-grid spacing, km
real, parameter :: dy = 0. ! y-grid spacing, km
real, parameter :: xlonc = 0. ! Standard longitude of projection
real, parameter :: truelat1 = 0. ! True latitudes of projection
real, parameter :: truelat2 = 0.! True latitudes of projection
real, parameter :: earth_radius = 6367470. * .001 ! Earth radius, km
logical, parameter :: is_wind_grid_rel = .false. ! Flag indicating whether winds are
! relative to source grid (TRUE) or
! relative to earth (FALSE)
character(len=8), parameter :: startloc = 'SWCORNER' ! Which point in array is given by
! startlat/startlon; set either
! to 'SWCORNER' or 'CENTER '
character(len=9), parameter :: field = 'SST' ! Name of the field
character(len=25) :: units = 'K' ! Units of data
character(len=32) :: map_source = 'JCOPE2' !OISST' ! Source model / originating center
character(len=46) :: desc = 'SST' ! Short description of data
real, parameter :: sst_intv = 24. ! input SST interval hours
! dynamic output variables
real :: xfcst ! Forecast hour of data (hours from ifile1)
real, allocatable :: isst1(:,:), isst2(:,:), osst(:,:) ! The 2-d array holding the data
real, allocatable :: lon(:), lat(:)
character(len=24) :: hdate ! Valid date for data YYYY:MM:DD_HH:00:00
! other variables
character(256) :: ifile1, ifile2, ofile
character(80) :: cxfcst ! String of xfcst
integer, parameter :: iunit1 = 10, iunit2 = 20, ounit = 100
character(500)::out_check,out_ctl
integer yyyy,mm,dd,hh
character(3)::cmm(12)
data cmm/'jan','feb','mar','apr','may','jun','jul','aug',&
& 'sep','oct','nov','dec'/
real ,parameter::fillvalue=-999.0
! interplation in time and extrapolation in space
real,allocatable ::itsst(:,:)
integer,parameter :: nw=5
real,allocatable ::tj(:,:,:) !JCOPE2 3D Temperature data
character(4) :: fildscT
integer iyf,imf,idf,mall
!-------------------------------------------
namelist /para/ifile1,ifile2,ofile,hdate,cxfcst,&
& yyyy,mm,dd,hh, &
& im,jm,km,tbias,xlons,xlons_offset_denom,ylats,&
& ylats_offset_denom,deltalat,deltalon
read(*,para)
read(cxfcst, *) xfcst
allocate(tj(im,jm,km),itsst(im,jm),isst1(im,jm),isst2(im,jm),&
& osst(im,jm),lon(im),lat(jm))
startlon =xlons-1.0/xlons_offset_denom
startlat =ylats-1.0/ylats_offset_denom
do i=1,im
lon(i)=startlon + deltalon*float(i-1)
enddo !i
do j=1,jm
lat(j)=startlat + deltalat*float(j-1)
enddo !i
! read JCOPE2 Temp data
tj(:,:,:)=0.0
print '(A,A,A)','Reading ',trim(ifile1),' ...'
open(iunit1, file=ifile1,form='unformatted',action="read")
read(iunit1)&
& fildscT,iyf,imf,idf &
& ,(((tj(i,j,k),i=1,im),j=1,jm),k=1,km-1),mall
close(iunit1)
isst1(:,:)=tj(:,:,1)+tbias
where(isst1 == tbias) isst1 = fillvalue
print '(A)','Done.'
tj(:,:,:)=0.0
print '(A,A,A)','Reading ',trim(ifile2),' ...'
open(iunit2, file=ifile2,form='unformatted',action="read")
read(iunit2)&
& fildscT,iyf,imf,idf &
& ,(((tj(i,j,k),i=1,im),j=1,jm),k=1,km-1),mall
close(iunit2)
isst2(:,:)=tj(:,:,1)+tbias
where(isst2 == tbias) isst2 = fillvalue
print '(A)','Done.'
! interpolate sst in time
itsst(:,:) = (isst1(:,:) )*(sst_intv-xfcst)/sst_intv + (isst2(:,:) )*xfcst/sst_intv
where(itsst <= -10.0) itsst = fillvalue
osst(:,:)=itsst(:,:)
! do extrapolation using spline technique
! call my_near_neighbor(itsst,osst,im,jm,nw,fillvalue)
! write data
print '(A,A)','Output: ',trim(ofile)
open(ounit, file=ofile)
do j=1,jm
do i=1,im
write(ounit,'(f11.5,f10.5,f8.2)' ) lon(i),lat(j),osst(i,j)
enddo !i
enddo !j
close(ounit)
end program jcope2sst
::::::::::::::
my_near_neighbor.f90
::::::::::::::
subroutine my_near_neighbor &
(z,zout,nx,ny,nw,fillvalue)
! Description:
!
! Author: am
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work1/am/WRF.Pre/OISST.Daily.Mean.Takatama
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 16:45 on 02-21-2014.
integer,intent(in)::nx,ny,nw
real,intent(in)::z(nx,ny),fillvalue
real,intent(inout)::zout(nx,ny)
zout(:,:)=z(:,:)
! print *,im,jm,ni
do j=1+nw,ny-nw
do i=1+nw,nx-nw
if(z(i,j) /= fillvalue)then
zout(i,j)=z(i,j)
! print *,bbb
goto 100
end if
do n=1,nw
do jj=-n,n
do ii=-n,n
if(z(i+ii,j+jj) /= fillvalue)then
zout(i,j)=z(i+ii,j+jj)
! print *,aaa
goto 100
end if
enddo !ii
enddo !jj
enddo !n
100 continue
end do ! i
enddo ! j
end subroutine my_near_neighbor
::::::::::::::
makefile
::::::::::::::
#
#
#j Fortran77とFortran90/95の混成プログラムのコンパイル用makefile
#
#j 作成 Sun Feb 15 15:16:02 2009
#
#j 注：ごく簡単な動作チェックしか行っていない
#j サフィックスルールがこれで正しいかどうかチェックする必要あり
#
#j 参考文献
#j makefile の書き方
#j http://www.eis.osakafu-u.ac.jp/~yabu/soft/makefile.html
#j Makefile小技と新バージョンのDelFEM(FEMとUIの日記@New York)
#j http://d.hatena.ne.jp/etopirika5/20091207/1260207955
#
#-----------------------------------------
#j マクロの定義
#-----------------------------------------
#j コンパイラの指定
FC=ifort
#FC=gfortran
#FC=g77
#FC=f77
#j OBJDIR (オブジェクトファイルをおくディレクトリ)
OBJDIR=../obj
#-----------------------------------------
#j コンパイルオプション(ifort)
#-----------------------------------------
#j デバッグ用オプション(はじめて実行するときには必ずこのオプションをつけてコンパイルする)
#FDFLAGS= # -CB -traceback -fpe0
#j コンパイルオプション
FFLAGS=-module ${OBJDIR}
#j -module : モジュールファイル(.mod)を置く場所を指定する
#j 最適化用(他にも何種類かある. ifort -helpかマニュアルを見て調べる)
#FFLAGS = -O3 -module ${OBJDIR}
#j 倍精度
#FFLAGS = -r8 -module ${OBJDIR}
#j 入力データのバイト・スワップ(大型計算機のバイナリデータを読む時などにつかう)
FFLAGS = -O3 -convert big_endian -module ${OBJDIR}
#-----------------------------------------
#j リンク用のオプション
#-----------------------------------------
LDFLAGS=-module ${OBJDIR}
#j ターゲット名(最終的に作りたい実行ファイル名)
TARGET1=jcope2sst.exe
TARGET2=
TARGETS=$(TARGET1) $(TARGET2)
#j MOD (Fortran90以降のモジュールファイル)
MOD=
#j SRC7 (Fortran77のソースファイル)
SRC7=
#j SRC9 (Fortran90のソースファイル)
SRC9=jcope2sst.f90 my_near_neighbor.f90
#j OBJ1 (ターゲットを作るのに必要なオブジェクトファイル名)
OBJM=$(MOD:.f90=.o)
OBJ7=$(SRC7:.f=.o)
OBJ9=$(SRC9:.f90=.o)
OBJTMP=$(OBJM) $(OBJ7) $(OBJ9)
OBJ1=$(OBJTMP:%=${OBJDIR}/%)
MOD2=
#j SRC7 (Fortran77のソースファイル)
SRC72=
#j SRC9 (Fortran90のソースファイル)
SRC92=
#j OBJ (ターゲットを作るのに必要なオブジェクトファイル名)
OBJM2=$(MOD2:.f90=.o)
OBJ72=$(SRC72:.f=.o)
OBJ92=$(SRC92:.f90=.o)
OBJTMP2=$(OBJM2) $(OBJ72) $(OBJ92)
OBJ2=$(OBJTMP2:%=${OBJDIR}/%)
OBJ=$(OBJ1) $(OBJ2)
#-----------------------------------------
#j ここからコンパイルのルールの記述
#-----------------------------------------
all: mkobjd $(TARGETS)
$(TARGET1): $(OBJ1)
$(FC) -o $@ $(OBJ1) ${LDFLAGS}
$(TARGET2): $(OBJ2)
$(FC) -o $@ $(OBJ2) ${LDFLAGS}
#j 暗黙のサフィックスルールを無効にする
.SUFFIXES:
#j サフィックスの追加
.SUFFIXES: .o .f .f90
#j f77のソースファイルのコンパイル(.f用のサフィックスルール)
${OBJDIR}/%.o : %.f
$(FC) -c ${FDFLAGS} ${FFLAGS} -c -o $@ $<
#j f90のソースファイルのコンパイル(.f90用のサフィックスルール)
${OBJDIR}/%.o : %.f90
$(FC) -c ${FDFLAGS} ${FFLAGS} -c -o $@ $<
#j オブジェクトファイルを削除
clean:
rm -rf $(OBJ) ${OBJDIR}/*.mod
#j オブジェクトファイルと実行ファイルを削除
distclean:
rm -rf $(OBJ) ${OBJDIR}/*.mod $(TARGETS) \
rm -rf ${OBJDIR}
#j オブジェクトファイルを置くディレクトリを作成
mkobjd:
if [ ! -d ${OBJDIR} ]; then \
mkdir -p ${OBJDIR} ; \
fi
run:
cd .. ; \
./run.sh 2>&1 > log.txt ; \
cd src
#j tarファイルを作る
tar:
tar cvf $(TARGET1).tar ./*.f90 ./*.f ./*.c ./*.h ./*.txt makefile
::::::::::::::
pl.sst.3.sh
::::::::::::::
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd30
# Directory: /work2/am/12.Work11/21.Climatology/16.jcope_tend_adv/plot
#
# Revision history:
# This file is created by /usr/local/mybin/ngmt.sh at 18:21 on 07-27-2011.
usage()
{
echo Error in $0 : Wrong number of arguments
echo Usage : $0 [-n] input_file run_name
echo -n: Do not paint land
echo
return
}
CMDNAME=$(basename $0)
while getopts n OPT; do
case $OPT in
"n" ) flagn="true" ;;
* ) usage; exit 1 ;;
esac
done
# オプションが指定されなかった場合, flaglに"NOT_AVILABLE"という値を入れる。
flagn=${flagn:-"NOT_AVAILABLE"}
shift `expr $OPTIND - 1`
. ./gmtpar.sh
echo "Bash script $0 starts."
scale=1.
range1=115/150/23/46
size=5
xanot=a5f1
yanot=a5f1
reso=15m/15m
reso_smth=30m/30m
reso_mask=30m/30m
tension=1
anot=${xanot}/${yanot}
tension=1
cpt=$(basename $0 .sh).cpt.txt
cpt_range="10/32/1"
makecpt -Cno_green -T$cpt_range > $cpt
cpt_range="-5/5/.5"
cpt2=$(basename $0 .sh).cpt.2.txt
makecpt -Cjet -T$cpt_range > $cpt2
plcontour=n
CI_tnd=1
CI_RHS=$CI_tnd
scalex=2.5
scaley=-0.5
scalesize=2.5
scale_anot=a10f2
if [ $# -ne 1 ]; then
usage
exit 1
fi
in=$1
#runname=$2
if [ ! -f $in ]; then
echo Error in $0 : No such file, $in
exit 1
fi
psdir="../ps/"
if [ ! -d ${psdir} ];then
mkdir -p $psdir
fi
out=$psdir$(basename $in .txt).2.ps
echo Input : $in
echo Output : $out
grd=$1.grd #$(basename $in .asc).grd
tstamp1=$(ls -l $in |awk '{print $6,$7,$8}')
title=""
col=3
awk -v scl=$scale '{if (($'"$col"' > -100.) && ($'"$col"' < 1000.)) \
print $1,$2,$'"$col"' }' $in | \
blockmean -R$range1 -I$reso_smth | \
surface -R$range1 -I$reso -T0.8 -G$grd
awk -v scl=$scale '{if (($'"$col"' > -100.) && ($'"$col"' < 1000.)) \
print $1,$2,$'"$col"' }' $in | \
psmask -R$range1 -I$reso_mask -JM$size \
-K -X1.5 -Y2 -P > $out
grdimage $grd -R$range1 -JM \
-C$cpt \
-O -K \
>> $out
grdcontour $grd -R$range1 -JM -A${CI_tnd}f10 -G2/2 -L-100/100 \
-W5/${white} -O -K >> $out
psmask -C -O -K >> $out
if [ $flagn = "true" ]; then
pscoast -R$range1 -JM -B${anot}WSne -Dl -W2 -O -K >> $out
else
pscoast -R$range2 -JM -B${anot}WSne -Dl -W2 -G200 -O -K >> $out
fi
text="116 41.7 18 0 1 LM $title"
pstext << end -R -JM -O -K >> $out
$text
end
psscale -D$scalex/$scaley/$scalesize/0.1h -B${scale_anot}:"${deg}C": -C$cpt -O -K >> $out
xoffset=
yoffset=4
comment=
. ./note.sh
rm -rf ${grd} ${cpt}
echo "Done $0.sh"
::::::::::::::
gmtpar.sh
::::::::::::::
#
# GMTのパラメータの設定 (version 4以上対象）
#
# 長さの単位はインチとする
gmtset MEASURE_UNIT INCH
gmtset LABEL_FONT_SIZE 18
gmtset HEADER_FONT_SIZE 20
gmtset ANOT_FONT_SIZE 18
gmtset TICK_PEN 4
# 地図の縦横軸に縞々を使わない
gmtset BASEMAP_TYPE PLAIN
# 紙のサイズはA4
gmtset PAPER_MEDIA A4
# 地図の°の記号
gmtset DEGREE_SYMBOL = degree # Available for ver. 4 or later
# 空白文字の設定
sp="\040" # White space
# =の記号
eq="\075" # equal
# 温度の°の記号
deg="\260" #deg="\312" # degree symbol
# 色のRGB値を設定 (線や点に色をつけるときRGB値を直接指定するより分かりやすい)
red="255/0/0"
orenge="255/165/0"
pink="255/192/203"
green="0/255/0"
blue="0/0/255"
midnightblue="25/25/112"
yellow="255/255/0"
gold="255/215/0"
purple="160/32/240"
magenta="255/0/255"
white="255/255/255"
# 線種を指定
dash="t15_5:15"
dot="t5_5:5"
dotdash="t12_5_5_5:5"
#その他の線種の例
# -W5t20_5:5
# -W5t15_5:5
# -W5t20_5_5_5_5_5:5
# -W5t30_5_5_5:5
# -W5t20_5_10_5:5
# 数字の出力フォーマット(project, grd2xzy等で使う)
gmtset D_FORMAT %lg #デフォルト(規定値)
# 桁数を指定（例：123.45678)
#gmtset D_FORMAT %12.6f
# 状況に応じて
#gmtset D_FORMAT %12.5g
::::::::::::::
note.sh
::::::::::::::
# Print time, current working directory and output filename
export LANG=C
currentdir=`pwd`
if [ ${#currentdir} -gt 90 ]; then
curdir1=${currentdir:1:90}
curdir2=${currentdir:91}
else
curdir1=$currentdir
curdir2="\ "
fi
now=`date`
host=`hostname`
#comment=" "
time1=$(ls -l $in | awk '{print $6, $7, $8}')
pstext -JX6/1.2 -R0/1/0/1.2 -N -X${xoffset:-0} -Y${yoffset:-0} << EOF -O >> $out
0 1.1 9 0 1 LM $0 $@
0 0.95 9 0 1 LM ${now}
0 0.80 9 0 1 LM ${host}
0 0.65 9 0 1 LM ${curdir1}
0 0.50 9 0 1 LM ${curdir2}
0 0.35 9 0 1 LM Input: ${in} (${time1})
0 0.1 9 0 1 LM ${comment:-""}
EOF
# Output: ${out}